この質問は、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)
結果のプロットは次のようになり、信頼帯は前後に跳ね返ります。
どこかで道に迷ったのではないでしょうか (行列の掛け算でしょうか?)。これが良いアイデアかどうかについての提案を含め、どんな助けも大歓迎です!