r - 更适合线性模型

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;博士 添加 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 summary 不会以任何其他方式表明可能有问题。但是，如果我们尝试获得置信区间，我们会发现可能存在问题:

 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(.))