1

这个问题建立在从 lme fit 中提取预测带的基础上，但使用非线性混合模型。

我有一个按“条目”分组的响应值数据集。我使用 AIC 模型选择程序来测试哪种类型的模型（线性、对数、指数等）最能代表响应和预测变量之间的关系。现在我想绘制每个条目中数据的拟合值，然后跨条目绘制。我还想在整体趋势上绘制随附的置信带——请参阅 Ben Bolker 博客上提供的代码和上面的帖子（我知道这是对解释的一些技巧，但那是另一篇帖子）。后者是我遇到问题的地方 - 请参阅此示例代码：

#Load nlme
library(nlme)

#Create data frame
set.seed(6)
df=data.frame(y=c(1:5+runif(5,0,1),21:25+runif(5,1,5)),
              x=rep(1:5,2),entry=rep(letters[1:2],each=5))

#Group data by entry
df=groupedData(y~x|entry,data=df)

#Build model
mod=nlme(y~a+b*x,fixed=a+b~1,random=~a+b~1,start=c(a=1,b=1),data=df)

#Create data frame for predictions
pred=do.call(rbind,lapply(rownames(coef(mod)),function(i) {
  data.frame(x=seq(0.1,max(df[df$entry==i,"x"]),0.1),entry=i) } ) )
#Grab coefficients from model for a and b
a=coef(mod)[match(pred$entry,rownames(coef(mod))),1]
b=coef(mod)[match(pred$entry,rownames(coef(mod))),2]
#Grab fixed coefficients for overall trend and fitted values for individual entries using a and b above
pred$fixed=fixef(mod)[1]+fixef(mod)[2]*pred$x 
pred$fitted=a+pred$x*b

#Get SEs on predictions using code from Ben Bolker's website
predvar=diag(as.matrix(data.frame(a,b)) %*% vcov(mod)%*% t(as.matrix(data.frame(a,b))))
pred$fixed.SE=sqrt(predvar)
pred$fixed.SE2=sqrt(predvar+mod$sigma^2)

#Plot
ggplot()+geom_point(data=df,aes(x=x,y=y),col="grey50",alpha=0.5)+
  geom_line(data=subset(pred,max(df$x)>x & x>1),aes(x=x,y=fitted,group=entry),col="grey50")+
  geom_ribbon(data=subset(pred,max(df$x)>x & x>1),aes(x=x,ymin=fixed-2*fixed.SE,ymax=fixed+2*fixed.SE),alpha=0.5,fill="grey10")+
  geom_ribbon(data=subset(pred,max(df$x)>x & x>1),aes(x=x,ymin=fixed-2*fixed.SE2,ymax=fixed+2*fixed.SE2),alpha=0.3,fill="grey50")+
  geom_line(data=subset(pred,max(df$x)>x & x>1),aes(x=x,y=fixed),col="black",lwd=1.1)

结果图如下所示，置信带来回跳跃：

在此处输入图像描述

我怀疑我在某个地方出了差错（也许在矩阵乘法中？）。感谢您提供任何帮助，包括有关这是否是个好主意的建议！

4

0 回答 0