RCBD土壤入渗数据：R中Philips/Horton模型及代码报错解决

问题 #

处理RCBD（随机完全区组设计）土壤水分入渗实验数据，包含处理组和区组效应，变量有入渗深度（CI）、时间（Time）、干/湿体积含水量（mDry、mWet），无导水率和吸渗率数据。尝试拟合通用非线性模型时报错，希望得到Philips/Horton模型或其他入渗率计算的代码建议。

数据片段 #

Block,Treatment,Time,CI,mDry,mWet
Block1,Control 1,0,0,0.0299,0.564
Block1,Control 1,48,10,0.0299,0.564
Block1,Control 1,217,20,0.0299,0.564
Block1,Control 1,848,30,0.0299,0.564
Block1,Control 1,2288,40,0.0299,0.564
Block1,Control 1,3911,50,0.0299,0.564
Block1,Control 1,5855,60,0.0299,0.564
Block1,Control 1,8624,70,0.0299,0.564
Block1,Control 1,11472,80,0.0299,0.564
Block1,Control 1,13581,90,0.0299,0.564
Block1,Control 1,16724,100,0.0299,0.564
Block1,Biochar25kgPha,0,0,0.027355623,0.299869622
Block1,Biochar25kgPha,90,10,0.027355623,0.299869622
Block1,Biochar25kgPha,200,20,0.027355623,0.299869622
Block1,Biochar25kgPha,341,30,0.027355623,0.299869622
Block1,Biochar25kgPha,604,40,0.027355623,0.299869622
Block1,Biochar25kgPha,1028,50,0.027355623,0.299869622
Block1,Biochar25kgPha,1486,60,0.027355623,0.299869622
Block1,Biochar25kgPha,2168,70,0.027355623,0.299869622
Block1,Biochar25kgPha,3030,80,0.027355623,0.299869622
Block1,Biochar25kgPha,4363,90,0.027355623,0.299869622
Block1,Biochar25kgPha,5461,100,0.027355623,0.299869622

尝试的代码及报错 #

尝试的通用非线性模型代码：

InfilMod <- nlme::gnls(CI ~ (Kf * (t ^ alpha)) / alpha,
data = Infil,
weights = nlme::varIdent(form = ~1 | Treatment),
method = "ML",
control = nlme::glsControl(maxIter = 1000),
correlation = nlme::corAR1(form = ~Time | Block/Treatment),
na.action = na.exclude)

报错信息：

Error in diff[pp != 0] <- pp[pp != 0] * ssize : NAs are not allowed in subscripted assignments 

数据无空值/缺失，仅存在t=0时的零值。

解决方案与代码建议 #

一、修正通用非线性模型报错 #

报错源于t=0时模型表达式的数值不稳定，以及corAR1处理Time=0分组时的异常。修正步骤如下：

1. 移除t=0的观测：入渗模型通常从t>0开始拟合，t=0的CI=0是边界条件，无需参与参数估计；
2. 统一变量名：原代码用t，数据中变量为Time，需保持一致；
3. 添加参数初始值：提升拟合稳定性。

修正后的代码：

# 过滤Time>0的观测
Infil_filtered <- Infil[Infil$Time > 0, ]

# 拟合修正后的通用非线性模型
InfilMod <- nlme::gnls(CI ~ (Kf * (Time ^ alpha)) / alpha,
                       data = Infil_filtered,
                       weights = nlme::varIdent(form = ~1 | Treatment),
                       method = "ML",
                       control = nlme::glsControl(maxIter = 1000),
                       correlation = nlme::corAR1(form = ~Time | Block/Treatment),
                       na.action = na.exclude,
                       start = list(Kf = 0.1, alpha = 0.5))

二、Philips模型（ssPhilips）拟合 #

Philips模型表达式：CI = S*sqrt(Time) + Kf*Time，其中S为吸渗率，Kf为稳定导水率。适配RCBD设计的代码示例：

library(nlme)

# 过滤Time>0的数据
Infil_filtered <- Infil[Infil$Time > 0, ]

# 拟合带随机效应、方差结构和自相关的Philips模型
Philips_mod <- nlme(CI ~ SSphilips(Time, S, Kf),
                    data = Infil_filtered,
                    random = pdDiag(S + Kf ~ 1 | Block/Treatment),
                    weights = varIdent(form = ~1 | Treatment),
                    correlation = corAR1(form = ~Time | Block/Treatment),
                    method = "ML",
                    control = nlmeControl(maxIter = 1000))

# 查看模型结果
summary(Philips_mod)

三、Horton模型（ssHorton）拟合 #

Horton模型表达式：CI = (i0 - Kf)/k * (1 - exp(-k*Time)) + Kf*Time，其中i0为初始入渗率，Kf为稳定导水率，k为衰减系数。代码示例：

library(nlme)

# 过滤Time>0的数据
Infil_filtered <- Infil[Infil$Time > 0, ]

# 拟合Horton模型，设置初始值提升拟合成功率
Horton_mod <- nlme(CI ~ SShorton(Time, i0, Kf, k),
                   data = Infil_filtered,
                   random = pdDiag(i0 + Kf + k ~ 1 | Block/Treatment),
                   weights = varIdent(form = ~1 | Treatment),
                   correlation = corAR1(form = ~Time | Block/Treatment),
                   method = "ML",
                   control = nlmeControl(maxIter = 1000),
                   start = list(i0 = 5, Kf = 0.1, k = 0.001))

# 查看模型结果
summary(Horton_mod)

四、入渗率计算 #

入渗率为CI对Time的导数，可通过模型预测值计算数值导数：

# 以Philips模型为例，生成预测时间序列
pred_time <- seq(10, max(Infil_filtered$Time), by = 10)
pred_data <- expand.grid(Time = pred_time, 
                         Block = unique(Infil_filtered$Block),
                         Treatment = unique(Infil_filtered$Treatment))

# 预测CI值
pred_data$CI_pred <- predict(Philips_mod, newdata = pred_data)

# 计算入渗率（数值导数）
pred_data$infiltration_rate <- c(diff(pred_data$CI_pred)/diff(pred_data$Time), NA)

内容的提问来源于stack exchange，提问作者Segeana