# CIs for random intercepts from gam and stan_gamm4 are similar, but those from gamm4 are not. why?

Question: how can one obtain a good confidence interval for the estimated random effects from gamm4?

Motivation: The example below, using binomial data with a random intercept, shows that estimated random effects are similar with gam, gamm4, and stan_gamm4 but the CI obtained from the fitted gamm4 model by using "postVar" is much narrower than the CI obtained from gam with the corrected covariance matrix, which is similar to the credible interval using stan_gamm4. I assume that the CI from "postVar" is wrong. Why is it wrong and how should it be fixed?

R
pk<-c("mgcv","gamm4","rstanarm")
lapply(pk,library,character.only=T)

set.seed(1)
dat <- gamSim(1,n=1000,scale=2)
dat$$fac <- fac <- as.factor(sample(1:20,1000,replace=TRUE)) rn<-2*rnorm(20) dat$$y <- dat$$y + model.matrix(~fac-1)%*%rn h<-function(u) (1+exp(-u))^-1 p<-sapply(scale(dat$$y),h)
T<-rpois(1000,40)
v<-rbinom(1000,size=T,prob=p)
vmat<-cbind(v,T-v)
datb<-data.frame(vmat,x0=dat$$x0,x1=dat$$x1,x2=dat$$x2,fac=dat$$fac)

br <- gamm4(vmat~s(x0)+x1+s(x2),data=datb,family="binomial",random=~(1|fac))
plot(br$gam,pages=1) bg <- gam(vmat~s(x0)+x1+s(x2)+s(fac,bs="re"), data=datb,family="binomial",method="REML") plot(bg,pages=1) cbr<-ranef(br$$mer,condVar=T)$$fac[,1] ocr<-order(cbr) cbg<-coef(bg)[21:40] plot(cbg,cbr) lines(c(-2,2),c(-2,2),lty=2,col="blue")  #ok sebr<-(attr(ranef(br$$mer,condVar=TRUE)$$fac,"postVar")[1,1,])^0.5 sebg<-(diag(bg$Vc)[21:40])^0.5
plot(sebr,sebg)


#different but related

#stan_gamm4 ##slow, about 30 mins to fit sg4 on my pc

zx0<-with(dat,scale(x0))
zx1<-with(dat,scale(x1))
zx2<-with(dat,scale(x2))
datz<-data.frame(zx0,zx1,zx2,vmat,fac)
pairs(sg4,pars=c("s(zx0).1","s(zx2).1"))
plot_nonlinear(sg4)
mat<-as.matrix(sg4)
dim(mat)
colnames(mat)
matc<-mat[,21:40]

q025<-function(v) quantile(v,0.025)
q975<-function(v) quantile(v,0.975)
mn<-apply(matc,2,mean)
lo<-apply(matc,2,q025)
hi<-apply(matc,2,q975)

plot(cbr,mn)
lines(range(cbr),range(cbr),lty=2,col="seagreen")

plot(qx,cbr[ocr],ylim=range(rn-mean(rn)),ylab="estimated random effect and 95% ci",main="binomial data: gamm4, gam, and stan_gamm4")
for(i in 1:20)
{
lines(rep(qx[i],2),cbr[ocr[i]]+qnorm(c(0.025,0.975))*sebr[ocr[i]])
lines(rep(qx[i]+0.02,2),cbg[ocr[i]]+qnorm(c(0.025,0.975))*sebg[ocr[i]],lty=2,col="blue")
lines(rep(qx[i]+0.04,2),c(lo[ocr[i]],hi[ocr[i]]),lty=3,col="red")
points(qx[i]+0.04,mn[ocr[i]],col="red",pch=16)
}
points(qx,rn[ocr]-mean(rn),pch=16,cex=0.8)
abline(h=0,col="blue",lty=1,lwd=0.5)
lines(c(-2,-2),c(2,2.5))
points(-2,2.25)
text(-2,2,pos=1,cex=0.8,label="gamm4\nest & ci")
lines(c(-1.8,-1.8),c(2,2.5),lty=2,col="blue")
text(-1.8,2,pos=1,label="gam\nci",cex=0.8)
lines(c(-1.6,-1.6),c(2,2.5),lty=3,col="red")
text(-1.6,2,pos=1,label="stan_gamm4\nest & ci",cex=0.8)
points(-1.6,2.25,col="red",pch=16)
points(-1.2,2.1,pch=16,cex=0.8)
text(-1.2,2,pos=1,label="true random effect",cex=0.8)


#good agreement between the ci from gam (dashed blue) and stan_gamm4 (dotted red) #but not gamm4 (black)

• Please register &/or merge your accounts (you can find information on how to do this in the My Account section of our help center), then you will be able to edit & comment on your own question. Jun 13, 2022 at 0:31

here is an answer in the simpler context of lme4, gam, and stan_lmer. I do not know if something like this would work for gamm4 though.

postVar is the variance of the estimated random effects conditional on the estimated fixed effects.

This variance appears in formula 2.15 on p80 of Simon Wood's book "Generalized Additive Models", which gives the conditional distribution.

Further down the page, 2.17 shows the joint distribution of the random and fixed effects.

The joint variance turns out to be identical with the uncorrected posterior variance of the gam model, using bs="re" in the smoother for the random effects.

CIs from the (slightly larger) corrected posterior variance, taking account of smoothing parameter uncertainty, are close to the credible intervals from stan_lmer.

But CIs from the uncorrected posterior variance are not much different, and are identical with those obtained in lmer from the joint distribution. CIs from postVar are much smaller.

All of this can be seen in an example, below. The notation in Wood differs from that in lme4, but iA as defined below is the variance in 2.15, and iM is the variance in 2.17

R
pk<-c("mgcv","lme4","rstanarm")
lapply(pk,library,character.only=T)

set.seed(1)
x1<-runif(400)
x2<-runif(400)
y<-1.8+0.7*x1-0.3*x2+1.2*rnorm(400)
dat <-data.frame(y,x1,x2)
dat$$fac <- fac <- as.factor(sample(1:20,400,replace=TRUE)) rn<-1.5*rnorm(20) dat$$y <- dat$y + model.matrix(~fac-1)%*%rn m1<-lmer(y~x1+x2+(1|fac),data=dat) Z<-getME(m1,"Z") X<-getME(m1,"X") Lambda<-getME(m1,"Lambda") sigma<-getME(m1,"sigma") psi<-sigma^2*Lambda%*%t(Lambda) A<-sigma^-2*t(Z)%*%Z+solve(psi) A<-as.matrix(A) pV<-attr(ranef(m1)$fac,"postVar")[1,1,]
iA<-solve(A)
range(diag(iA)-pV)

B<-sigma^-2*t(X)%*%X
B<-as.matrix(B)
C<-sigma^-2*t(Z)%*%X
C<-as.matrix(C)

M<-matrix(0,23,23)
M[1:20,1:20]<-A
M[21:23,21:23]<-B
M[1:20,21:23]<-C
M[21:23,1:20]<-t(C)
iM<-solve(M)

g1<-gam(y~x1+x2+s(fac,bs="re"),method="REML",data=dat)
Vp<-g1$Vp  to compare with iM, we need to permute the uncorrected posterior variance Vp so the fixed effects come after the smooth term o<-c(4:23,1:3) Vpx<-Vp[o,o] range(iM-Vpx)  so the variance of the joint distribution from lmer (2.17) is, within numerical error, the same as the uncorrected posterior variance in the gam model. now consider the corrected posterior variance and compare the CIs as in the original question est_m1<-ranef(m1)$$fac[,1] est_g1<-g1$$coef[4:23] range(est_m1-est_g1) mean(est_g1) plot(est_g1,rn-mean(rn),main="estimated vs true random effects") lines(range(est_g1),range(est_g1),lty=2,col="blue") se_pV<-pV^0.5 se_Vp<-(diag(Vp)[4:23])^0.5 Vc<-g1$Vc
se_Vc<-(diag(Vc)[4:23])^0.5

sm1<-stan_lmer(y~x1+x2+(1|fac),data=dat)
mat<-as.matrix(sm1)
matc<-mat[,4:23]
q025<-function(v) quantile(v,0.025)
q975<-function(v) quantile(v,0.975)
lo<-apply(matc,2,q025)
hi<-apply(matc,2,q975)
mn<-apply(matc,2,mean)

qx<-qnorm((1:20-0.5)/20)
og<-order(est_g1)
plot(qx,est_g1[og],main="Caterpillar plot with CIs", ylim=c(-4,4),ylab="estimated random effects and CIs",xlab="quantiles")
for (i in 1:20)
{
lines(rep(qx[i],2),est_g1[og[i]]+qnorm(0.975)*c(-1,1)*se_pV[og[i]])
points(qx[i]+0.02,est_g1[og[i]],pch=16,cex=0.8)
lines(rep(qx[i]+0.02,2),est_g1[og[i]]+qnorm(0.975)*c(-1,1)*se_Vp[og[i]],lty=2)
points(qx[i]+0.04,est_g1[og[i]],pch=16,col="blue",cex=0.8)
lines(rep(qx[i]+0.04,2),est_g1[og[i]]+qnorm(0.975)*c(-1,1)*se_Vc[og[i]],lty=3,col="blue")
points(qx[i]+0.06,mn[og[i]],pch=17,cex=0.8,col="red")
lines(rep(qx[i]+0.06,2),c(lo[og[i]],hi[og[i]]),lty=3,col="red")
points(qx[i],(rn-mean(rn))[og[i]],pch=16,col="seagreen")
}

lines(c(-2,-2),c(2,3))
points(-2,2.5)
text(-2,2,pos=1,cex=0.8,label="lmer\nest & ci")

lines(c(-1.7,-1.7),c(2,3),lty=2)
points(-1.7,2.5,pch=16,cex=0.7)
text(-1.7,2,pos=1,label="gam\nest &\n uncorrected ci",cex=0.8)

lines(c(-1.4,-1.4),c(2,3),lty=3,col="blue")
points(-1.4,2.5,pch=16,cex=0.7,col="blue")
text(-1.4,2,pos=1,label="gam\nest &\n corrected ci",cex=0.8)

lines(c(-1.1,-1.1),c(2,3),lty=3,col="red")
points(-1.1,2.5,col="red",pch=17)
text(-1.1,2,pos=1,label="stan_lmer\nest & ci",cex=0.8)

points(-0.7,2.1,pch=16,cex=1,col="seagreen")
text(-0.7,2,pos=1,label="true random effect",cex=0.8)

text(-2,0.8,pos=4,"lmer joint ci = gam uncorrected ci")