[R-sig-ME] question about confidence intervals with glmmPQL

[mintoc at mathstat.dal.ca]

Hi Sally,

> For two models the CIs are not symmetrical around the mean proportion
from the model

The confidence intervals, as implemented in the code you sent, are
symmetric on the linear predictor (logit) scale and not necessarily
symmetric on the response scale.

Also, see Bolker et al. (2009) for a discussion of the performance of PQL
with low numbers of successes and failures; and Zeger et al. (1988) for
extracting marginal parameters in a logistic mixed model such as yours.

Best regards,
Coilin

Bolker, B.M., Brooks, M.E., Clark, C.J., Geange, S.W., Poulsen, J.R.,
Stevens, M.H.H., and White, J.S. (2009). Generalized linear mixed models:
a practical guide for ecology and evolution. Trends in Ecology and
Evolution, 24(3), 127-135.

Zeger, S.L., Liang, K.-Y. and Albert, P.S. (1988). Models for longitudinal
data: a generalized estimating equation approach. Biometrics, 44,
1019-1031.

> Hi - I am using R version 2.13.0. I have run several GLMMs using the
> glmmPQL function in the MASS package to model the proportion of fish
> caught in one net to the total caught in both nets by length for 5
> different species. I started with a polynomial regression full model with
> three length terms: l, l^2, and l^3 (l=length). The length terms and
> intercept were the fixed effects and the random effect was a paired haul
> (n=18).
>
> For the majority of the models, I ended up with a constant model with no
> length effect. I used a step wise backward selection for length terms and
> terms were removed if they were not significant.
>
> The issue I am having is with the confidence intervals that were
> calculated. For two models the CIs are not symmetrical around the mean
> proportion from the model. The CIs for the other constant models are
> symmetrical around the mean. I was wondering if anyone has an idea why
> this would be or if anyone has any suggestions.
>
> I provided a small sample of data for one species along with the CI,
> model, and plotting codes. I used some code that was passed to me from a
> colleague who got it from someone else. I have gone through the code to
> understand it, but got stuck on a couple of areas. I did not hear back
> from the original person who wrote the code.
>
> #subsample of data
> Pair Species Length Experimental Control Proportion second third
> 2 Monkfish 23 0 1 0 529 12167
> 2 Monkfish 27 0 1 0 729 19683
> 7 Monkfish 27 0 2 0 729 19683
> 8 Monkfish 27 0 1 0 729 19683
> 14 Monkfish 28 1 0 1 784 21952
> 13 Monkfish 29 0 1 0 841 24389
> 14 Monkfish 29 0 1 0 841 24389
> 1 Monkfish 30 0 1 0 900 27000
> 5 Monkfish 30 0 2 0 900 27000
> 18 Monkfish 30 1 0 1 900 27000
> 4 Monkfish 31 1 0 1 961 29791
> 5 Monkfish 31 0 1 0 961 29791
> 11 Monkfish 31 0 1 0 961 29791
> 2 Monkfish 32 0 1 0 1024 32768
> 11 Monkfish 32 0 1 0 1024 32768
> 12 Monkfish 32 0 1 0 1024 32768
> 18 Monkfish 32 1 1 0.5 1024 32768
> 4 Monkfish 34 0 1 0 1156 39304
> 17 Monkfish 34 0 1 0 1156 39304
> 18 Monkfish 34 0 1 0 1156 39304
> 18 Monkfish 35 0 1 0 1225 42875
> 1 Monkfish 37 0 1 0 1369 50653
> 2 Monkfish 40 1 0 1 1600 64000
> 14 Monkfish 40 0 1 0 1600 64000
> 18 Monkfish 40 1 0 1 1600 64000
>
>
> This is my model code
> #Models
> #Select trip
> t<-t1
> #select species
> species<-subset(t,Species=="Monkfish")
> #define polynomials
> #2nd order poly
> second<-species$Length^2
> #3rd order poly
> third<-species$Length^3
> #combine species,second & third
> species<-cbind(species,second,third)
>
> #backward selection
> #full model
> m1<-glmmPQL(fixed=Proportion~1+Length+second+third,random=~1|Pair,family=binomial,data=species,verbose=T,niter=2,weight=(Experimental+Control))
> summary(m1)
> #2nd order
> m2<-update(m1, . ~. -third)
> summary(m2)
> #linear
> m3<-update(m2, . ~. -second)
> summary(m3)
> #constant
> m4<-update(m3, . ~. -Length)
> summary(m4)
>
>
> Here is the code for a constant model -
>
> #CIs
> get.sel.and.conf.band.catch=function(parm,varm,l.min,l.max,n=200){
> lgt<-seq(l.min,l.max,length=n)
> X=sapply(1:length(parm),function(i){lgt^(i-1)})
> reg.line=X%*%parm
> se=apply(X,1,function(x,varm){sqrt(t(x)%*%varm%*%x)},varm)
> cbind(lengtht=lgt,min.L50=reg.line-2*se,mean.L50=reg.line,max.L50=reg.line+2*se)
> }
>
> #plots model results
> plot.fit=function(glmm.fit,limx,Spec.Text){
> varm=glmm.fit$varFix
> coeff=glmm.fit$coeff$fixed
> eta.matrix=get.sel.and.conf.band.catch(coeff,varm,min(limx),max(limx))
> xxx=c(eta.matrix[,1],rev(eta.matrix[,1]))
> yyy=c(ilogit(eta.matrix[,2]),rev(ilogit(eta.matrix[,4])))
> plot(eta.matrix[,1],ilogit(eta.matrix[,3])
> ,type="l"
> ,xlab="Length (cm)"
> ,ylab="Proportion: exp/(exp+control)"
> ,ylim=c(0,1), xlim=limx
> ,col="black",lwd=2,las=1
> ,axes=FALSE
> ,family="serif")
> axis(1,labels=T,cex.axis=0.9)
> #axis(3,labels=TRUE,cex.axis=0.9)
> axis(2,labels=TRUE,las=2,cex.axis=0.9)
> box()
> lines(eta.matrix [,1],ilogit(eta.matrix [,2]),col=4)
> lines(eta.matrix[,1],ilogit(eta.matrix[,4]),col=2)
> polygon(xxx,yyy,col="grey",density=NULL,border="black")
> lines(eta.matrix[,1],ilogit(eta.matrix[,3]),col=1,lwd=2)
> abline(h=0.5,lty=2)
> text(sum(limx)/2,0.9,Spec.Text,cex=1)
> }
>
> #model plot code
> plot.fit(m4,c(min(species$Length),max(species$Length)),paste(unique(species$Species)))
>
> I looked at other code online from Daniel Hocking and Ben Bolker at
> glmm.wikidot.com. There are differences between what I received and the
> examples from Daniel and Ben.
>
> Any help would be appreciated.
> Thanks Sally
>
> --
>
> Sally Roman
> Fisheries Biologist
> University of Massachusetts Dartmouth
> School for Marine Science and Technology
> 200 Mill Road Suite 325
> Fairhaven, MA 02917
> phone 508-910-6379
>
> 	[[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>