更适合线性模型

Question

我正在拟合一些线条，我觉得我在告诉 R 究竟如何拟合它们，但我觉得有一些我不知道的东西（某些因素或影响）阻碍了良好的拟合。

我的实验单元是“地块”，就像田间地块一样，对不起，这令人困惑。

数据可以查到：https://www.dropbox.com/s/a0tplyvs8lxu1d0/rootmeansv2.csv。 与

df$plot.f<-as.factor(df$plot)
dfG<-groupedData(mass ~ year|plot.f, data=df)
dfG30<-dfG[dfG$depth == 30,]

简单地说，随着时间的推移，我有质量，我用模型将它拟合到每个实验单元：

fit <- lme(mass ~ year , random = ~ 1 | plot, data = df)

和plot (augPred(fit)) 我得到这些适合每个实验单元（“情节”）：

我需要怎么做才能让实验单元之间的斜率变化更大？从统计的角度来看，我对此不感兴趣，但从预测的角度来看——所以模型中的任何东西都可以被操纵来让这些线条移动。

Answer 1

这个答案有点像百科全书，因为关于你的问题有几个要点。 tl;dr 添加year*plot 作为随机效果是第一步，但拟合实际上有点问题，虽然一开始看起来不是这样：将year 变量居中问题，但是对数据进行日志转换可能是一个更好的主意。

更新：OP 正在对subset(df,Depth==30) 进行分析。我不小心对整个数据集进行了分析，这可能使事情变得更加困难——下面记录的异方差性对于数据的子集可能并没有那么糟糕——但出于教学目的，我将保持原样（出于懒惰）。

获取数据：

df <- read.csv("rootmeansv2.csv")
library(nlme)
gdf <- groupedData( mass ~ year | plot, data=df)

将逐年的交互作为随机效应添加到模型中：

fit0 <- lme(mass ~ year , random = ~ year | plot, data = gdf)

但是，如果我们查看结果，我们会发现这根本没有太大帮助 - year 项（斜率中的绘图方差）很小。

##             Variance     StdDev       Corr  
## (Intercept) 9.522880e-12 3.085916e-06 (Intr)
## year        7.110616e-08 2.666574e-04 0.32  
## Residual    3.885373e-01 6.233276e-01       

这确实有时会发生（单数适合），但除此之外，lme 总结不会以任何其他方式表明可能有问题。但是，如果我们尝试获得置信区间，我们会发现可能会有问题：

 intervals(fit0)
 ##   Error in intervals.lme(fit0) : 
 ##   cannot get confidence intervals on var-cov components: Non-positive definite approximate variance-covariance

再次检查的另一种方法是在lme4 中拟合相同的模型。

library(lme4)
VarCorr(fit1 <- lmer(mass ~ year +(year | plot), data = gdf))

##  Groups   Name        Std.Dev.   Corr  
##  plot     (Intercept) 0.66159674       
##           year        0.00059471 -1.000
##  Residual             0.62324403       
## Warning messages:
## 1: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
##   Model failed to converge with max|grad| = 0.132739 (tol = 0.002, component 1)
## 2: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
##   Model failed to converge: degenerate  Hessian with 1 negative eigenvalues

我们得到了明显不同的答案，以及关于收敛的警告（这是在开发版本 1.1-7 中，它不受 1.1-6 广泛报道的误报警告的影响）。看起来lme4的合身稍微好一点：

c(logLik(fit1)-logLik(fit0))
## 0.1775161

复杂拟合（例如混合模型）的潜在数据问题之一是预测变量的中心化和缩放不良。在这种情况下，因为year 是从 2008 年开始的连续预测变量，所以截距和斜率非常高度相关（截距是第 0 年的预测值！）。在这种情况下，我们可以通过居中来解决问题——减去最小值也是合理的（即从 0 开始而不是 2008 年）

c. <- function(x) scale(x,center=TRUE,scale=FALSE)
VarCorr(fit2 <- update(fit1,.~ c.(year) +(c.(year) | plot)))

##  Groups   Name        Std.Dev. Corr 
##  plot     (Intercept) 0.53798       
##           c.(year)    0.10515  1.000
##  Residual             0.59634       

我们得到了更合理的答案（并且没有警告），尽管我们仍然有完全（现在正）相关的截距和斜率——这只是说模型稍微过度拟合了数据，但我们对此无能为力（我们可以通过拟合模型~c.(year)+(c.(year)||plot)) 将相关性强制为零，但这有其自身的问题）。

现在去lme试试吧：

VarCorr(fit3 <- update(fit0,
                       fixed.=~c.(year),
                       random=~c.(year)|plot,
         control=lmeControl(opt="optim")))
## plot = pdLogChol(c.(year)) 
##             Variance   StdDev    Corr  
## (Intercept) 0.28899909 0.5375864 (Intr)
## c.(year)    0.01122497 0.1059479 0.991 
## Residual    0.35559015 0.5963138       

结果相似（尽管相关性仅为 0.991 而不是 1.0）：对数似然的差异实际上略大，但仍然很小。 （拟合仍然有些问题——我必须更改优化器，如lmeControl 参数所示才能到达这里。）

现在看起来好多了：

plot(augPred(fit3))

但是，残差与拟合图向您展示了一个您应该担心的问题：

plot(fit3)

漏斗形状表明您存在异方差问题。比例位置图显示得更清楚：

plot(fit3,sqrt(abs(resid(.)))~fitted(.),type=c("p","smooth"))

最明显的解决方法是对数据进行日志转换：

df$logmass <- log10(df$mass)  ## log10() for interpretability
gdfL <- groupedData( logmass ~ year | plot, data=df)
VarCorr(fitL1 <- lme(logmass ~ c.(year) , random = ~ c.(year) | plot,
             data = gdfL))
## plot = pdLogChol(c.(year)) 
##             Variance     StdDev     Corr  
## (Intercept) 0.1842905872 0.42929080 (Intr)
## c.(year)    0.0009702893 0.03114947 -0.607
## Residual    0.1052946948 0.32449144       

年际差异再次变小，但这一次可能是因为不需要。当我们拟合等效模型时，lmer 没有警告，我们得到相同的结果；拟合是非奇异的（没有零方差或完全相关）；和intervals(fitL1) 工作。

预测看起来很合理：

plot(augPred(fitL1))

诊断图也是如此：

plot(fitL1)

虽然 2008 年似乎确实有些有趣（轴被翻转是因为 lme 更喜欢水平绘制箱线图 -- factor(year) 告诉 lme 使用箱线图而不是散点图）。

plot(fitL1,factor(year)~resid(.))