library(faraway) data(ctsib) ctsib$stable <- ifelse(ctsib$CTSIB==1,1,0) summary(ctsib) gf <- glm(stable ~ Sex+Age+Height+Weight+Surface+Vision,binomial,data=ctsib) summary(gf) gfs <- glm(stable ~ Sex+Age+Height+Weight+Surface+Vision+factor(Subject),binomial,data=ctsib) anova(gf,gfs,test="Chi") library(brlr) modbr <- brlr(stable ~ Sex+Age+Height+Weight+Surface+Vision+factor(Subject), data=ctsib) summary(modbr) library(MASS) gg <- glmmPQL(stable ~ Sex+Age+Height+Weight+Surface+Vision,random=~1|Subject, family=binomial,data=ctsib) summary(gg) sd(coef(modbr)[9:43]) data(ctsib) ctsib$stable <- ifelse(ctsib$CTSIB==1,1,0) library(gee) gg <- gee(stable ~ Sex+Age+Height+Weight+Surface+Vision,id=Subject,family=binomial,data=ctsib,corstr="exchangeable",scale.fix=TRUE) summary(gg) data(epilepsy) epilepsy[1:10,] with(epilepsy,by(seizures/timeadj,list(treat,expind),mean)) y <- matrix(epilepsy$seizures,nrow=5) matplot(1:4,sqrt(y[-1,]),type="l",lty=epilepsy$treat[5*(1:59)]+1,xlab="Period",ylab="Sqrt(Seizures)") my <- apply(y[-1,],2,mean)/2 plot(sqrt(epilepsy$seizures[epilepsy$expind == 0]/8),sqrt(my),pch=epilepsy$treat[5*(1:59)]+2,xlab="sqrt(Baseline seizures)",ylab="sqrt(Experiment seizures)") abline(0,1) g <- gee(seizures ~offset(log(timeadj))+expind+treat+I(expind*treat), id,family=poisson,corstr="AR-M",Mv=1,data=epilepsy,subset=(id!=49)) summary(g)