[R] gam plots and seWithMean

[Greg Dropkin]

hadn't realised the answer would be in the source code!

anyway, this appears to work. The only difference is in the last section.

greg

--

library(mgcv)

#simulate some data
x1<-runif(500)
x2<-rnorm(500)
x3<-rpois(500,3)
d<-runif(500)
t<-runif(500,20,50)
linp<--6.5+x1+2*x2-x3+2*exp(-2*d)*sin(2*pi*d)
lam<-t*exp(linp)
y<-rpois(500,lam)
sum(y)
table(y)

#fit the data without d by glm and with d by gam
f1<-glm(y~offset(log(t))+x1+x2+x3,poisson)
f2<-gam(update.formula(as.formula(f1),~.+s(d)),poisson)
anova(f1,f2)
summary(f2)
plot(f2)

#the solid line s(d)
dat<-data.frame(t,d,x1,x2,x3)
datn<-transform(dat,d=0)
dif<-predict(f2)-predict(f2,datn)
cdif<-dif-mean(dif)
points(d,cdif,cex=0.5,col=rgb(0,1,0,0.2))

#another approach to the solid line s(d)
devAskNewPage(ask=T)
plot(f2)
premat2<-PredictMat(f2$smooth[[1]],data=dat)
dim(premat2)
pars<-f2$coef
pars2<-pars[5:13]
pars2<-as.matrix(pars2,9,1)
pars2
points(d,premat2%*%pars2,cex=0.5,col=rgb(0,0.6,0.3,0.2))
#premat2%*%pars2 = cdif

#confidence intervals when seWithMean = FALSE
devAskNewPage(ask=T)
plot(f2)
points(d,cdif,cex=0.5,col=rgb(0,1,0,0.2))
Vp2<-f2$Vp[5:13,5:13]
se2<-sqrt(diag(premat2%*%Vp2%*%t(premat2)))
points(d,cdif+qnorm(0.975)*se2,cex=0.5,col=rgb(1,0,0,0.2))
points(d,cdif-qnorm(0.975)*se2,cex=0.5,col=rgb(0,0,1,0.2))
#numerical output for the confidence bands is given by
#cdif+qnorm(0.975)*se2
#cdif-qnorm(0.975)*se2

#confidence intervals when seWithMean = TRUE
devAskNewPage(ask=T)
plot(f2,seWithMean=T)
points(d,cdif,cex=0.5,col=rgb(0,1,0,0.2))
premat<-cbind(rep(1,500),x1,x2,x3,premat2)
pars<-as.matrix(pars,13,1)
range(predict.gam(f2)-log(t)-as.numeric(premat%*%pars))
#so predict.gam(f2) = log(t) + as.numeric(premat%*%pars)
sew<-sqrt(diag(premat%*%f2$Vp%*%t(premat)))
range(sew-predict(f2,se.fit=T)$se.fit)
#sew = predict(f2,se.fit=T)$se.fit
points(d,cdif+qnorm(0.975)*sew,cex=0.5,col=rgb(1,0,0,0.2))
points(d,cdif-qnorm(0.975)*sew,cex=0.5,col=rgb(0,0,1,0.2))
#so this is not what plot.gam is doing

#but this is
devAskNewPage(ask=T)
plot(f2,seWithMean=T)
points(d,cdif,cex=0.5,col=rgb(0,1,0,0.2))
premat<-cbind(rep(1,500),x1,x2,x3,premat2)
pars<-as.matrix(pars,13,1)
premat1<-premat
premat1[,1:4]<-rep(f2$cmX[1:4],each=500)
sew1<-sqrt(diag(premat1%*%f2$Vp%*%t(premat1)))
points(d,cdif+qnorm(0.975)*sew1,cex=0.5,col=rgb(1,0,0,0.2))
points(d,cdif-qnorm(0.975)*sew1,cex=0.5,col=rgb(0,0,1,0.2))