################################################### ### chunk number 1: ################################################### library(lattice) options(width=65,show.signif.stars=FALSE) options(SweaveHooks=list(fig=function() { ltheme <- canonical.theme(color = FALSE) ## in-built B&W theme ltheme$strip.background$col <- "transparent" ## change strip bg ltheme$plot.symbol$col <- "black" ltheme$axis.text$cex <- 1.2 ltheme$par.xlab.text$cex <- 1.2 ltheme$par.ylab.text$cex <- 1.2 lattice.options(default.theme = ltheme) ## set as default })) #options(width=65,show.signif.stars=FALSE) #library(xtable) #library(lattice) #ltheme <- canonical.theme(color = FALSE) ## in-built B&W theme #ltheme$strip.background$col <- "transparent" ## change strip bg #lattice.options(default.theme = ltheme) ## set as default ################################################### ### chunk number 2: ################################################### faults <- read.csv("Data/faults.csv") ################################################### ### chunk number 3: faults.glm ################################################### #data("faults",package="SMIR") faults <- transform(faults, ll = log(l)) faults.glm <- glm(n ~ ll, data = faults, family = poisson) round(coef(summary(faults.glm)),4) ################################################### ### chunk number 4: faults.glm1 ################################################### faults.glm1 <- glm(n ~ offset(ll), data = faults, family = poisson) anova(faults.glm,faults.glm1, test = "Chisq") ################################################### ### chunk number 5: faults.figure eval=FALSE ################################################### ## xval <- with(faults,seq(min(l),max(l),len=100)) ## fval <- exp(coef(faults.glm1)[1]+log(xval)) ## fval.sd <- sqrt(fval) ## tval <- qt(0.975, faults.glm1$df.residual) ## lb <- fval - tval * fval.sd ## ub <- fval + tval * fval.sd ## print(xyplot(fval+lb+ub~xval,col="black",type="l", ## xlab="length",ylab="number",lty=c(1,3,3),ylim=c(0,30), ## panel=function(...){ ## panel.xyplot(...) ## lpoints(faults$l,faults$n,col="black")})) ################################################### ### chunk number 6: ################################################### xval <- with(faults,seq(min(l),max(l),len=100)) fval <- exp(coef(faults.glm1)[1]+log(xval)) fval.sd <- sqrt(fval) tval <- qt(0.975, faults.glm1$df.residual) lb <- fval - tval * fval.sd ub <- fval + tval * fval.sd print(xyplot(fval+lb+ub~xval,col="black",type="l", xlab="length",ylab="number",lty=c(1,3,3),ylim=c(0,30), panel=function(...){ panel.xyplot(...) lpoints(faults$l,faults$n,col="black")})) ################################################### ### chunk number 7: claims ################################################### claims <- read.csv('Data/claims.csv') claims$age <- factor(claims$age, labels=c(" <25","25-29","30-35"," >35")) claims$dist <- factor(claims$dist, labels=c("rural","small towns","large towns","major cities")) claims$car <- factor(claims$car, labels=c("<1","1-1.5","1.5-2",">2")) ################################################### ### chunk number 8: eval=FALSE ################################################### ## data(claims,data="SMIR") ################################################### ### chunk number 9: claims.anova ################################################### claims <- transform(claims, ln=log(n)) claims.glm <- glm(c ~ offset(ln) + dist + car + age + car:age + dist:car + dist:age, data = claims, family = poisson) anova(claims.glm) ################################################### ### chunk number 10: ################################################### library(xtable) xtable(anova(claims.glm), caption="Analysis of deviance table for the claims data ", label="tab:claims.anova") ################################################### ### chunk number 11: claims.glm1 ################################################### claims.glm1 <- update(claims.glm, . ~ offset(ln) + dist + car + age) claims.res <- round(resid(claims.glm1, type = "pearson"), 3) claims.res[rev(order(abs(claims.res)))] ################################################### ### chunk number 12: claims.glm2 ################################################### claims.glm2 <- update(claims.glm1, . ~ ln + dist + car + age) print(anova(claims.glm1, claims.glm2),digits=4) round(coef(summary(claims.glm2))[2,,drop=FALSE],4) ################################################### ### chunk number 13: ################################################### round(coef(summary(claims.glm1)),4) ################################################### ### chunk number 14: claims.glm4 ################################################### claims <- transform(claims,d4=factor(dist == "major cities")) claims.glm4 <- update(claims.glm1, . ~ . - dist + d4) anova(claims.glm4, claims.glm1) round(coef(summary(claims.glm4)),4) ################################################### ### chunk number 15: claims.glm6 ################################################### claims <- transform(claims, lcar=unclass(car)) claims.glm5 <- update(claims.glm4, . ~ . - car + lcar) anova(claims.glm5,claims.glm4) round(coef(summary(claims.glm5)),4) claims <- transform(claims,lage=unclass(age)) claims.glm6 <- update(claims.glm5, . ~ . - age + lage) anova(claims.glm6, claims.glm5) round(coef(summary(claims.glm6)),4) ################################################### ### chunk number 16: claims.glm7 ################################################### claims <- transform(claims,score = lcar - lage) claims.glm7 <- update(claims.glm6, . ~ offset(ln)+d4 + score) round(coef(summary(claims.glm7)),4) ################################################### ### chunk number 17: claims.glm8 ################################################### claims <- transform(claims,s2 = as.numeric(d4)-1 + score) claims.glm8 <- update(claims.glm7, . ~ offset(ln)+s2) round(coef(summary(claims.glm8)),4) ################################################### ### chunk number 18: claims_table ################################################### source('claims.r') ################################################### ### chunk number 19: ################################################### miners <- read.csv("Data/miners.csv") miners <- transform(miners, t = n + m + s) miners <- transform(miners, np=n/t, mp=m/t, sp=s/t) ################################################### ### chunk number 20: eval=FALSE ################################################### ## data(miners) ## miners <- transform(miners, t = n + m + s) ## miners <- transform(miners, np=n/t, mp=m/t, sp=s/t) ################################################### ### chunk number 21: miners.plot eval=FALSE ################################################### ## miners.plot <- xyplot(np+mp+sp~years,data=miners,col=1, ## ylab="proportion",type="p",pch=c("n","m","s")) ## print(update(miners.plot,type="b",cex=1.5,lty=1:3)) ## #print(xyplot(np+mp+sp~years,data=miners,col=1,cex=2, ## #ylab="proportion",type="b",pch=c("n","m","s"),lty=c(1,2,3))) ################################################### ### chunk number 22: ################################################### miners.plot <- xyplot(np+mp+sp~years,data=miners,col=1, ylab="proportion",type="p",pch=c("n","m","s")) print(update(miners.plot,type="b",cex=1.5,lty=1:3)) #print(xyplot(np+mp+sp~years,data=miners,col=1,cex=2, #ylab="proportion",type="b",pch=c("n","m","s"),lty=c(1,2,3))) ################################################### ### chunk number 23: miners.glm ################################################### miners <- transform(miners,period=factor(years)) miners.glm <- glm(cbind(n, (t - n)) ~ period, family = binomial, data = miners) summary(miners.glm) ################################################### ### chunk number 24: miners.glm1 ################################################### miners.glm1 <- update(miners.glm,. ~ years) #summary(miners.glm1) print(anova(miners.glm1),digits=4) round(resid(miners.glm1,type='pearson'),3) ################################################### ### chunk number 25: miners.glm1a ################################################### miners.glm1a <- update(miners.glm1, family = binomial(link = "cloglog")) #summary(miners.glm1a) print(anova(miners.glm1a),digits=4) round(resid(miners.glm1a,type='pearson'),3) ################################################### ### chunk number 26: miners.glm2 ################################################### miners <- transform(miners, ly=log(years)) miners.glm2 <- update(miners.glm1a,.~ly) summary(miners.glm2) round(resid(miners.glm2,type="pearson"),3) ################################################### ### chunk number 27: ################################################### miners$cnp <- fitted(miners.glm2) ################################################### ### chunk number 28: miners.glm3 ################################################### miners.glm3 <- update(miners.glm2,.~ly,family=binomial) summary(miners.glm3) round(resid(miners.glm3,type='pearson'),3) ################################################### ### chunk number 29: miners.plot2 eval=FALSE ################################################### ## miners$lnp <- fitted(miners.glm3) ## print(xyplot(np+cnp+lnp~years,type=c("p","l","l"),data=miners,ylab="proportion", ## panel=panel.superpose.2,lty=c(1,2))) ################################################### ### chunk number 30: miners.plot3 eval=FALSE ################################################### ## print(xyplot(np~years,data=miners,ylab='proportion',ylim=c(0,1.2),col="black", ## panel=function(x,y){ ## panel.xyplot(x,y) ## panel.curve(1-exp(-exp(4.582-1.312*log(x)))) ## panel.curve((eta <- exp(9.609-2.576*log(x)))/(1+eta),lty=2) ## })) ################################################### ### chunk number 31: ################################################### print(xyplot(np~years,data=miners,ylab='proportion',ylim=c(0,1.2),col="black", panel=function(x,y){ panel.xyplot(x,y) panel.curve(1-exp(-exp(4.582-1.312*log(x)))) panel.curve((eta <- exp(9.609-2.576*log(x)))/(1+eta),lty=2) })) ################################################### ### chunk number 32: miners.empirical.logits ################################################### miners <- transform(miners,logit2 = log(m/n), logit3 = log(s/n)) print(miners[c("years","logit2","logit3")]) ################################################### ### chunk number 33: miners.logit.plot eval=FALSE ################################################### ## print(xyplot(logit2+logit3~years,type="b",pch=c("m","s"), ## data=miners,cex=1.5, ylab="logit",lty=c(1,2),lwd=1.8)) ################################################### ### chunk number 34: ################################################### miners$lnp <- fitted(miners.glm3) print(xyplot(np+cnp+lnp~years,type=c("p","l","l"),data=miners,ylab="proportion", panel=panel.superpose.2,lty=c(1,2))) ################################################### ### chunk number 35: ################################################### print(xyplot(logit2+logit3~years,type="b",pch=c("m","s"), data=miners,cex=1.5, ylab="logit",lty=c(1,2),lwd=1.8)) ################################################### ### chunk number 36: miners.log.plot eval=FALSE ################################################### ## print(xyplot(logit2+logit3~years,type="b",pch=c("m","s"),col="black",lty=c(2,3),lwd=1.5, ## data=miners,cex=1.5,ylab="logit",xlab="years (log scale)",scale=list(x=list(log=TRUE, ## at=c(10,20,30,40,50),labels=c("10","20","30","40","50"))))) ################################################### ### chunk number 37: ################################################### print(xyplot(logit2+logit3~years,type="b",pch=c("m","s"),col="black",lty=c(2,3),lwd=1.5, data=miners,cex=1.5,ylab="logit",xlab="years (log scale)",scale=list(x=list(log=TRUE, at=c(10,20,30,40,50),labels=c("10","20","30","40","50"))))) ################################################### ### chunk number 38: clinic ################################################### clinic <- data.frame(death=c(3,4,17,2),survive=c(176,293,197,23), clinic=gl(2,2),attend=gl(2,1)) clinic <- transform(clinic,total = death + survive) clinic.glm <- glm(cbind(death,survive)~1,data=clinic, family=binomial) clinic.glma <- update(clinic.glm,.~clinic) clinic.glmb <- update(clinic.glma,.~.+attend) clinic.glmc <- update(clinic.glm,~clinic:attend-1) # #anova(clinic.glm) #summary(update(clinic.glm,.~1))$coeff #summary(update(clinic.glm,.~clinic))$coeff ################################################### ### chunk number 39: clinic.glm4 ################################################### pclinic <- reshape(clinic,varying=list(c('survive','death')), v.names='count',timevar='resp',times=c('survive','death'), direction='long') pclinic <- transform(pclinic, resp = factor(resp)) clinic.glm1 <- glm(count~resp+clinic:attend,family=poisson,data=pclinic) round(coef(summary(clinic.glm1)),4) clinic.glm2 <- update(clinic.glm1,.~.+resp:clinic) #anova(clinic.glm1,clinic.glm2) round(coef(summary(clinic.glm2)),4) clinic.glm3 <- update(clinic.glm2,.~.+resp:attend) #anova(clinic.glm2,clinic.glm3) clinic.glm4 <- update(clinic.glm3,.~.+resp:attend:clinic) #anova(clinic.glm3,clinic.glm4) anova(clinic.glm3) ################################################### ### chunk number 40: ################################################### print(with(pclinic,tapply(fitted(clinic.glm1), list(clinic=clinic,attend=attend),sum))) print(with(pclinic,tapply(count, list(clinic=clinic,attend=attend),sum))) ################################################### ### chunk number 41: clinic.gnm ################################################### library(gnm) clinic.gnm <- gnm(count~resp,eliminate=clinic:attend,data=pclinic, family=poisson,verbose=FALSE) summary(clinic.gnm) clinic.gnm1 <- update(clinic.gnm, .~resp+resp:clinic) summary(clinic.gnm1) ################################################### ### chunk number 42: miners.multinomial ################################################### miners.long <- reshape(miners,drop=names(miners)[6:11], varying=list(names(miners)[2:4]), timevar="resp", times=c("n","m","s"), v.names="count", direction="long") miners.long <- transform(miners.long, resp=factor(resp)) miners.long <- transform(miners.long, resp = relevel(resp,ref="n"), lyr = log(years),period = factor(years)) miners.glm5 <- glm(count~period+resp,data=miners.long,family=poisson) #miners.gnm1 <- gnm(count~resp,eliminate=period,data=miners.long,family=poisson,ofInterest="resp") coef(summary(miners.glm5))[9:10,,drop=FALSE] #summary(miners.gnm1) miners.glm6 <- update(miners.glm5,.~.+resp:years) #miners.gnm2 <- update(miners.gnm1,.~.+resp:years) anova(miners.glm5,miners.glm6) #anova(miners.gnm1,miners.gnm2) round(summary(miners.glm6)$coef[9:12,],4) #round(coef(summary(miners.glm6))[9:12,,drop=FALSE],4) #summary(miners.gnm2) #getContrasts(miners.gnm2,ofInterest(miners.gnm2)[5:3]) ################################################### ### chunk number 43: miners.glm7 ################################################### miners.glm7 <- update(miners.glm6,.~.+years+resp:years) #miners.gnm3 <- update(miners.gnm2,.~.+years+resp:years) print(anova(miners.glm7,miners.glm6),digits=5) #anova(miners.6nm2,miners.gnm3) round(coef(summary(miners.glm7))[9:12,],4) #summary(miners.gnm3) miners.long$fval7 <- fitted(miners.glm7) miners.long$resid7 <- resid(miners.glm7) options(digits=3) cbind(miners.long[1:8,c('count','fval7','resid7')], miners.long[9:16,c('count','fval7','resid7')], miners.long[17:24,c('count','fval7','resid7')], row.names=NULL) ################################################### ### chunk number 44: miners.glm8 ################################################### miners.glm8 <- update(miners.glm6,.~period+resp+lyr+resp:lyr) #miners.gnm4 <- update(miners.gnm3,.~resp*lyr) print(anova(miners.glm6,miners.glm8),digits=4) #anova(miners.gnm2,miners.gnm4) print(coef(summary(miners.glm8))[9:12,,drop=FALSE],digits=5) #summary(miners.gnm4) miners.long$fval8 <- fitted(miners.glm8) miners.long$resid8 <- resid(miners.glm8) options(digits=3) cbind(miners.long[1:8,c('count','fval8','resid8')], miners.long[9:16,c('count','fval8','resid8')], miners.long[17:24,c('count','fval8','resid8')],row.names=NULL) ################################################### ### chunk number 45: miners.plot5 eval=FALSE ################################################### ## miners.long$fv8 <- fitted(miners.glm8) ## #miners.long$fv <- fitted(miners.gnm4) ## miners.long <- transform(miners.long, ## f2p8=fv8/t, ## op = count/t) ## print(miners.plot <- xyplot(op~years,data=miners.long,groups=resp, ## ylab='proportion',xlab='years',pch=c("n","m","s"),cex=1.5)) ## print(update(miners.plot, ## panel=function(x,y,...){ ## panel.xyplot(x,y,...) ## panel.xyplot(x,miners.long$f2p8,type="l",...)})) ################################################### ### chunk number 46: ################################################### miners.long$fv8 <- fitted(miners.glm8) #miners.long$fv <- fitted(miners.gnm4) miners.long <- transform(miners.long, f2p8=fv8/t, op = count/t) print(miners.plot <- xyplot(op~years,data=miners.long,groups=resp, ylab='proportion',xlab='years',pch=c("n","m","s"),cex=1.5)) print(update(miners.plot, panel=function(x,y,...){ panel.xyplot(x,y,...) panel.xyplot(x,miners.long$f2p8,type="l",...)})) ################################################### ### chunk number 47: ################################################### miners.long$respm <- relevel(miners.long$resp,ref="m") coef(summary(update(miners.glm8,.~period+respm*lyr)))[12,,drop=FALSE] #getContrasts(miners.gnm4,ofInterest(miners.gnm4)[3:4]) ################################################### ### chunk number 48: miners.glm9 ################################################### miners.long <- transform(miners.long,r23=ifelse(resp!="n",1,0)) miners.long <- transform(miners.long,r23lyr=r23*lyr) miners.glm9 <- update(miners.glm8,.~factor(years)+resp+r23lyr) #miners.gnm5 <- update(miners.gnm4,.~resp+r23:lyr) anova(miners.glm9,miners.glm8) #anova(miners.gnm5,miners.gnm4) print(coef(summary(miners.glm9))[9:11,,drop=FALSE],5) #coef(summary(miners.gnm5))[-(1:8),] ################################################### ### chunk number 49: miners.plot6 eval=FALSE ################################################### ## miners.long$fval9 <- fitted(miners.glm9) ## miners.long <- transform(miners.long, ## f2p9=fval9/t, ## op = count/t) ## print(update(miners.plot, ## panel=function(x,y,...){ ## panel.xyplot(x,y,...) ## panel.xyplot(x,miners.long$f2p9,type="l",...)})) ################################################### ### chunk number 50: ################################################### miners.long$fval9 <- fitted(miners.glm9) miners.long <- transform(miners.long, f2p9=fval9/t, op = count/t) print(update(miners.plot, panel=function(x,y,...){ panel.xyplot(x,y,...) panel.xyplot(x,miners.long$f2p9,type="l",...)})) ################################################### ### chunk number 51: miners.glm10 ################################################### #miners.gnm6 <- gnm(count~r23lyr+r23,eliminate=period,data=miners.long,family=poisson) #anova(miners.gnm5,miners.gnm6) #summary(miners.gnm6) miners.glm10 <- update(miners.glm9,.~.-resp+r23) anova(miners.glm10,miners.glm9) print(coef(summary(miners.glm10))[9:10,,drop=FALSE],digits=4) ################################################### ### chunk number 52: miners.plot7 eval=FALSE ################################################### ## miners.long <- transform(miners.long, ## fval10=fitted(miners.glm10)/t, ## resid10=resid(miners.glm10)) ## # ## print(update(miners.plot, ## panel=function(x,y,...){ ## panel.xyplot(x,y,...) ## panel.xyplot(x,miners.long$fval10,type="l",...)})) ################################################### ### chunk number 53: ################################################### miners.long <- transform(miners.long, fval10=fitted(miners.glm10)/t, resid10=resid(miners.glm10)) # print(update(miners.plot, panel=function(x,y,...){ panel.xyplot(x,y,...) panel.xyplot(x,miners.long$fval10,type="l",...)})) ################################################### ### chunk number 54: ################################################### miners.long$lres <- unclass(miners.long$resp)-1 ################################################### ### chunk number 55: mienrs.glm11 ################################################### miners.glm11 <- update(miners.glm10,.~factor(years)+resp+lyr+lres+lyr:lres) #miners.gnm7 <- gnm(count~resp+lyr:lres,data=miners.long, #eliminate=period,family=poisson) #anova(miners.gnm4,miners.gnm7) anova(miners.glm11,miners.glm8) print(coef(summary(miners.glm11))[9:11,,drop=FALSE],4) #summary(miners.gnm7) ################################################### ### chunk number 56: miners.glm12 ################################################### #miners.gnm8 <- update(miners.gnm7,.~.-resp+lres) miners.glm12 <- update(miners.glm11,.~.-resp+lres) anova(miners.glm12,miners.glm11) #anova(miners.gnm7,miners.gnm8) #summary(miners.gnm8) round(coef(summary(miners.glm12))[9:10,,drop=FALSE],4) ################################################### ### chunk number 57: ################################################### round(summary(update(miners.glm6,.~.-resp:years+lyr:lres), corr=TRUE)$corr[9:11,9:11,drop=FALSE],4) ################################################### ### chunk number 58: miners.cr ################################################### miners$notn <- miners$t-miners$n miners.cr <- reshape(miners,drop=names(miners)[c(6:8,10,11)], varying=list(names(miners)[2:3]),times=c("n","m"),timevar="resp", v.names="count",direction="long") miners.cr$other <- miners.cr$t-miners.cr$count miners.cr$other[miners.cr$resp=="m"] <- miners.cr$other[miners.cr$resp=="n"]- miners.cr$count[miners.cr$resp=="m"] miners.cr <- transform(miners.cr,lyears=log(years), resp=relevel(factor(resp),ref="n")) miners.cr[c("years","resp","count","other")] ################################################### ### chunk number 59: miners.cr.glm ################################################### miners.cr.glm <- glm(cbind(count,other)~resp, data=miners.cr,family=binomial) summary(miners.cr.glm) miners.cr.glm2 <- update(miners.cr.glm,.~.+lyears) summary(miners.cr.glm2) miners.cr.glm3 <- update(miners.cr.glm2,.~.+resp:lyears) summary(miners.cr.glm3) ################################################### ### chunk number 60: ################################################### miners.cr$fval <- fitted(miners.cr.glm3) round(cr.tab <- with(miners.cr,tapply(fval,list(Response=resp,Period=period),mean)),3) ################################################### ### chunk number 61: ################################################### fn <- cr.tab[1,] fm <- cr.tab[2,]*(1-fn) fs <- 1-fn-fm round(rbind(fn,fm,fs),3) ################################################### ### chunk number 62: miners.cr.plot eval=FALSE ################################################### ## #yval <- miners.cr$years[1:8] ## print(xyplot(np+mp+sp+fn+fm+fs~years,data=miners,type=c("p","p","p","l","l","l"),cex=1.5, ## panel=panel.superpose.2,ylab="proportion",xlab="year",pch=c("n","m","s"),lty=1:3,lwd=1.5)) ## #print(xyplot(fn+fm+fs~yval,col=1,cex=2, ## #ylab="proportion",xlab="year", ## #panel=function(x,y,...){ ## #panel.xyplot(x,y,type="l",lwd=1.5,...) ## #lpoints(miners$years,miners$np,type="p",pch="n",...) ## #lpoints(miners$years,miners$mp,type="p",pch="m",...) ## #lpoints(miners$years,miners$sp,type="p",pch="s",...) ## #})) ################################################### ### chunk number 63: ################################################### ################################################### ### chunk number 64: ################################################### #yval <- miners.cr$years[1:8] print(xyplot(np+mp+sp+fn+fm+fs~years,data=miners,type=c("p","p","p","l","l","l"),cex=1.5, panel=panel.superpose.2,ylab="proportion",xlab="year",pch=c("n","m","s"),lty=1:3,lwd=1.5)) #print(xyplot(fn+fm+fs~yval,col=1,cex=2, #ylab="proportion",xlab="year", #panel=function(x,y,...){ #panel.xyplot(x,y,type="l",lwd=1.5,...) #lpoints(miners$years,miners$np,type="p",pch="n",...) #lpoints(miners$years,miners$mp,type="p",pch="m",...) #lpoints(miners$years,miners$sp,type="p",pch="s",...) #})) ################################################### ### chunk number 65: ################################################### vietnam <- read.table('Data/vietnam.csv',header=TRUE) ################################################### ### chunk number 66: eval=FALSE ################################################### ## data(vietnam,package="SMIR") ################################################### ### chunk number 67: vietnam.observed.proportions eval=FALSE ################################################### ## oval <- aperm(array(vietnam$count,dim=c(4,5,2), ## dimnames=list(Policy=levels(vietnam$policy), ## Year=levels(vietnam$year), ## Sex=levels(vietnam$sex)[c(2,1)])),c(3,2,1)) ## vietnam.ob.tab <- prop.table(oval,c(1,2)) ## round(ftable(vietnam.ob.tab,row.vars=c(1,2)),3) ################################################### ### chunk number 68: ################################################### oval <- aperm(array(vietnam$count,dim=c(4,5,2), dimnames=list(Policy=levels(vietnam$policy), Year=levels(vietnam$year), Sex=levels(vietnam$sex)[c(2,1)])),c(3,2,1)) vietnam.ob.tab <- prop.table(oval,c(1,2)) round(ftable(vietnam.ob.tab,row.vars=c(1,2)),3) ################################################### ### chunk number 69: vietnam.gnm1 ################################################### #vietnam.glm1 <- glm(count~year:sex+policy,data=vietnam,family=poisson) vietnam.gnm1 <- gnm(count~policy,eliminate=year:sex,data=vietnam,family=poisson,verbose=FALSE) #anova(vietnam.glm1) #anova(vietnam.gnm1) ################################################### ### chunk number 70: vietnam.gnm2 ################################################### #vietnam.glm2 <- update(vietnam.glm1,.~.+policy:sex) vietnam.gnm2 <- update(vietnam.gnm1,.~.+policy:sex) #anova(vietnam.glm2,vietnam.glm1) #anova(vietnam.gnm1,vietnam.gnm2) #vietnam.glm3 <- update(vietnam.glm1,.~.+policy:year) vietnam.gnm3 <- update(vietnam.gnm1,.~.+policy:year) #anova(vietnam.glm3,vietnam.glm2) #anova(vietnam.gnm2,vietnam.gnm3) #vietnam.glm4 <- update(vietnam.glm3,.~.+policy:sex) vietnam.gnm4 <- update(vietnam.gnm3,.~.+policy:sex) #anova(vietnam.glm4,vietnam.glm3) print(anova(vietnam.gnm1,vietnam.gnm2,vietnam.gnm4),digits=5) print(anova(vietnam.gnm1,vietnam.gnm3,vietnam.gnm4),digits=5) ################################################### ### chunk number 71: vietnam.males ################################################### vietnam.males <- subset(vietnam,sex=="male") ################################################### ### chunk number 72: vietnam.males.gnm2 ################################################### #vietnam.males.glm1 <- glm(count~year+policy,data=vietnam.males, #family=poisson) vietnam.males.gnm1 <- gnm(count~policy,data=vietnam.males,eliminate=year, family=poisson,verbose=FALSE) #vietnam.males.glm2 <- update(vietnam.males.glm1,.~.+policy:year) vietnam.males.gnm2 <- update(vietnam.males.gnm1,.~.+policy:year) #anova(vietnam.males.glm1,vietnam.males.glm2) #anova(vietnam.males.gnm1,vietnam.males.gnm2) #round(coef(summary(vietnam.males.glm2))[6:20,,drop=FALSE],4) summary(vietnam.males.gnm2) ################################################### ### chunk number 73: vietnam.margins eval=FALSE ################################################### ## # par(mfrow=c(1,2)) ## # matplot(1:4, ## vietnam.mat1 <- data.frame(matrix(c(rep(0,4),coef(summary(vietnam.males.gnm2))[9:20]),nrow=4,byrow=T)) ## names(vietnam.mat1) <- c("Y1","Y2","Y3","Y4") ## vietnam.mat1$pol <- 1:4 ## #, ## # axes=FALSE,xlab='Policy',ylab='',type='l',lwd=1.5,col=1) ## # axis(2,las=1) ## # axis(1,at=1:4,LETTERS[1:4]) ## # legend(1.2,2,c(2,3,4,'Grad'),lty=1:4,col=1,lwd=1.5) ## # text(1,2.5,'(a)',pos=4) ## # box() ## # matplot(1:5, ## vietnam.mat2 <- data.frame(rbind(c(0,0,0),t(matrix(coef(summary(vietnam.males.gnm2))[9:20],ncol=4,byrow=T)))) ## #, ## # axes=FALSE,xlab='Year',ylab='',type='l',lwd=1.5,col=1,las=1) ## # axis(2,las=1) ## # axis(1,at=1:5,c(1,2,3,4,'Grad')) ## # legend(2,2,LETTERS[2:4],lty=1:3,col=1,lwd=1.5) ## # text(1,2.5,'(b)',pos=4) ## # box() ## vietnam.coef.plot1 <- xyplot(Y1+Y2+Y3+Y4~pol,data=vietnam.mat1,type="l", ## xlab="policy",ylab="",main="(a)",lty=1,col="black", ## scales=list(x=list(at=1:4,lab=LETTERS[1:4])), ## panel=function(...){ ## panel.xyplot(...) ## panel.text(3.8,c(0.4,0.8,1.2,2.4),c(2,3,4,"Grad"))}) ## vietnam.coef.plot2 <- xyplot(X1+X2+X3~1:5,data=vietnam.mat2,type="l", ## xlab="year",ylab="",main="(b)",lty=1,col="black", ## scales=list(x=list(at=1:5,labels=c(1,2,3,4,"Grad"))), ## panel=function(...){ ## panel.xyplot(...) ## panel.text(4.8,c(0.8,1.2,2.4),LETTERS[2:4])}) ## ## print(vietnam.coef.plot1,position=c(0,0,0.45,1),more=TRUE) ## print(vietnam.coef.plot2,position=c(0.42,0,1,1)) ################################################### ### chunk number 74: ################################################### # par(mfrow=c(1,2)) # matplot(1:4, vietnam.mat1 <- data.frame(matrix(c(rep(0,4),coef(summary(vietnam.males.gnm2))[9:20]),nrow=4,byrow=T)) names(vietnam.mat1) <- c("Y1","Y2","Y3","Y4") vietnam.mat1$pol <- 1:4 #, # axes=FALSE,xlab='Policy',ylab='',type='l',lwd=1.5,col=1) # axis(2,las=1) # axis(1,at=1:4,LETTERS[1:4]) # legend(1.2,2,c(2,3,4,'Grad'),lty=1:4,col=1,lwd=1.5) # text(1,2.5,'(a)',pos=4) # box() # matplot(1:5, vietnam.mat2 <- data.frame(rbind(c(0,0,0),t(matrix(coef(summary(vietnam.males.gnm2))[9:20],ncol=4,byrow=T)))) #, # axes=FALSE,xlab='Year',ylab='',type='l',lwd=1.5,col=1,las=1) # axis(2,las=1) # axis(1,at=1:5,c(1,2,3,4,'Grad')) # legend(2,2,LETTERS[2:4],lty=1:3,col=1,lwd=1.5) # text(1,2.5,'(b)',pos=4) # box() vietnam.coef.plot1 <- xyplot(Y1+Y2+Y3+Y4~pol,data=vietnam.mat1,type="l", xlab="policy",ylab="",main="(a)",lty=1,col="black", scales=list(x=list(at=1:4,lab=LETTERS[1:4])), panel=function(...){ panel.xyplot(...) panel.text(3.8,c(0.4,0.8,1.2,2.4),c(2,3,4,"Grad"))}) vietnam.coef.plot2 <- xyplot(X1+X2+X3~1:5,data=vietnam.mat2,type="l", xlab="year",ylab="",main="(b)",lty=1,col="black", scales=list(x=list(at=1:5,labels=c(1,2,3,4,"Grad"))), panel=function(...){ panel.xyplot(...) panel.text(4.8,c(0.8,1.2,2.4),LETTERS[2:4])}) print(vietnam.coef.plot1,position=c(0,0,0.45,1),more=TRUE) print(vietnam.coef.plot2,position=c(0.42,0,1,1)) ################################################### ### chunk number 75: vietnam.males.gnm3 ################################################### vietnam.males <- transform(vietnam.males,linp=unclass(policy)) #vietnam.males.glm3 <- update(vietnam.males.glm2,.~.-policy:year+linp+linp:year) vietnam.males.gnm3 <- update(vietnam.males.gnm2,.~.-policy:year+linp+linp:year) summary(vietnam.males.gnm3) #round(coef(summary(vietnam.males.glm3))[6:12,,drop=FALSE],4) ################################################### ### chunk number 76: vietnam.males.gnm4 ################################################### vietnam.males <- transform(vietnam.males,liny=ifelse(year=="G",7,unclass(year))) #vietnam.males.glm4 <- update(vietnam.males.glm3,.~.-year:linp+liny:linp) vietnam.males.gnm4 <- update(vietnam.males.gnm3,.~.-year:linp+liny:linp) #round(coef(summary(vietnam.males.glm4))[6:9,,drop=F],4) summary(vietnam.males.gnm4) ################################################### ### chunk number 77: vietnam.males.glm4.residuals ################################################### #round(residuals(vietnam.males.glm4,'pearson'),3) round(residuals(vietnam.males.gnm4,'pearson'),3) ################################################### ### chunk number 78: vietnam.males.glm5 ################################################### #vietnam.males.glm5 <- update(vietnam.males.glm4, # subset=(row.names(vietnam)!='14')) vietnam.males.gnm5 <- update(vietnam.males.gnm4, subset=(row.names(vietnam)!='14')) #round(coef(summary(vietnam.males.glm5))[6:9,,drop=F],4) summary(vietnam.males.gnm5) ################################################### ### chunk number 79: vietnam.females ################################################### vietnam.females <- subset(vietnam,sex=="female") #vietnam.females.glm1 <- glm(count~year+policy,data=vietnam.females, #family=poisson) vietnam.females.gnm1 <- gnm(count~policy,data=vietnam.females,eliminate=year, family=poisson) ################################################### ### chunk number 80: ################################################### summary(vietnam.females.gnm1) #round(coef(summary(vietnam.females.glm1))[6:8,,drop=FALSE],4) #round(residuals(vietnam.females.glm1,'pearson'),3) sort(round(residuals(vietnam.females.gnm1,'pearson'),3)) ################################################### ### chunk number 81: vietnam.gnm5 ################################################### vietnam <- transform(vietnam,liny=ifelse(year=="G",7,unclass(year)), linp=unclass(policy),female=sex=="female") vietnam <- transform(vietnam,v=(sex=="male")*liny*linp) #vietnam.glm5 <- glm(count~year*sex+policy+female:policy+v,data=vietnam, #family=poisson,subset=(row.names(vietnam)!='14')) vietnam.gnm5 <- gnm(count~policy+female:policy+v,data=vietnam,eliminate=year:sex, family=poisson,subset=(row.names(vietnam)!='14'),verbose=FALSE) summary(vietnam.gnm5) #anova(vietnam.glm5) #round(coef(summary(vietnam.glm5))[6:17,,drop=FALSE],4) ################################################### ### chunk number 82: ################################################### #fv <- predict(vietnam.glm5,new=vietnam[14,],type="response") #fv <- predict(vietnam.gnm5,new=vietnam[14,],type="response") #round((vietnam$count[14]-fv)/sqrt(fv),3) ################################################### ### chunk number 83: vietnam.fitted.plots eval=FALSE ################################################### ## vietnam$v <- as.numeric(vietnam$v) ## #vietnam$fv <- predict(vietnam.glm5,type="response",new=vietnam) ## vietnam.gnm5.14 <- update(vietnam.gnm5,data=vietnam,subset=row.names(vietnam)) ## #vietnam$fv <- predict(vietnam.glm5,type="response",new=vietnam) ## vietnam$fv <- fitted(vietnam.gnm5.14) ## #vietnam$t <- rep(t(tapply(vietnam$count,list(vietnam$sex,vietnam$year),sum)[c(2,1),]),rep(4,10)) ## vietnam$t <- rep(t(tapply(vietnam$count,list(vietnam$sex,vietnam$year),sum)[c(2,1),]),rep(4,10)) ## vietnam <- transform(vietnam,fp=fv/t,op=count/t) ## # ## print(xyplot(fp~liny|sex,data=vietnam, ## groups=policy, ## panel=function(x,y,type,subscripts,...) ## { ## lpoints(vietnam$liny[subscripts], ## vietnam$op[subscripts], ## pch=LETTERS[vietnam$policy[subscripts]],cex=1.5) ## # ## # pch=seq(1,4)[vietnam$policy[subscripts]]) ## panel.superpose(x,y,type="l",subscripts,...) ## #panel.text(x=5.5,c(0.1,0.2,0.3,0.45),labels=LETTERS[c(4,1,2,3)]) ## })) ## # ################################################### ### chunk number 84: ################################################### vietnam$v <- as.numeric(vietnam$v) #vietnam$fv <- predict(vietnam.glm5,type="response",new=vietnam) vietnam.gnm5.14 <- update(vietnam.gnm5,data=vietnam,subset=row.names(vietnam)) #vietnam$fv <- predict(vietnam.glm5,type="response",new=vietnam) vietnam$fv <- fitted(vietnam.gnm5.14) #vietnam$t <- rep(t(tapply(vietnam$count,list(vietnam$sex,vietnam$year),sum)[c(2,1),]),rep(4,10)) vietnam$t <- rep(t(tapply(vietnam$count,list(vietnam$sex,vietnam$year),sum)[c(2,1),]),rep(4,10)) vietnam <- transform(vietnam,fp=fv/t,op=count/t) # print(xyplot(fp~liny|sex,data=vietnam, groups=policy, panel=function(x,y,type,subscripts,...) { lpoints(vietnam$liny[subscripts], vietnam$op[subscripts], pch=LETTERS[vietnam$policy[subscripts]],cex=1.5) # # pch=seq(1,4)[vietnam$policy[subscripts]]) panel.superpose(x,y,type="l",subscripts,...) #panel.text(x=5.5,c(0.1,0.2,0.3,0.45),labels=LETTERS[c(4,1,2,3)]) })) # ################################################### ### chunk number 85: ################################################### policy.tab <- with(vietnam,tapply(count,list(year,sex,policy),sum)) Atab <- apply(policy.tab[,,1],c(1,2),sum) Btab <- apply(policy.tab[,,2],c(1,2),sum) Ctab <- apply(policy.tab[,,3],c(1,2),sum) Dtab <- apply(policy.tab[,,4],c(1,2),sum) BCDtab <- apply(policy.tab[,,2:4],c(1,2),sum) CDtab <- apply(policy.tab[,,3:4],c(1,2),sum) vietnam.cr <- expand.grid(year=factor(c(1:4,"G")),sex=c("female","male"),policy=c("A","B","C")) vietnam.cr$d1 <- c(Atab,Btab,Ctab) vietnam.cr$d2 <- c(BCDtab,CDtab,Dtab) vietnam.cr.glm <- glm(cbind(d1,d2)~policy+sex+year+sex:year,data=vietnam.cr,family=binomial) anova(vietnam.cr.glm) vietnam.cr.glm1 <- update(vietnam.cr.glm,.~.+policy:sex) summary(vietnam.cr.glm1) vietnam.cr.glm2 <- update(vietnam.cr.glm1,.~.+policy:year) summary(vietnam.cr.glm2) ################################################### ### chunk number 86: ################################################### #vietnam.cr$linyear <- as.numeric(as.character(vietnam.cr$year)) vietnam.cr$linyear <- with(vietnam.cr,ifelse(year=="G",7,unclass(year))) ################################################### ### chunk number 87: ################################################### vietnam.cr.glm3 <- update(vietnam.cr.glm2,.~.-policy:year+policy:linyear) anova(vietnam.cr.glm3,vietnam.cr.glm2) summary(vietnam.cr.glm3) ################################################### ### chunk number 88: vietnam.cr.glm3.residuals ################################################### round(resid(vietnam.cr.glm3),3) ################################################### ### chunk number 89: ################################################### vres <- resid(vietnam.cr.glm3) ################################################### ### chunk number 90: ################################################### vietnam.cr.glm4 <- update(vietnam.cr.glm3,.~.-policy:linyear) anova(vietnam.cr.glm4,vietnam.cr.glm3,test="Chisq") ################################################### ### chunk number 91: vietnam.cr.glm6 ################################################### vietnam.cr.glm5 <- update(vietnam.cr.glm4,.~.-sex:year+sex:linyear) anova(vietnam.cr.glm5,vietnam.cr.glm4) vietnam.cr.glm6 <- update(vietnam.cr.glm5,.~.-year+linyear) anova(vietnam.cr.glm6,vietnam.cr.glm5) coef(summary(vietnam.cr.glm6)) ################################################### ### chunk number 92: vietnam.cr.glm7 ################################################### vietnam.cr <- transform(vietnam.cr,male.linyear = ifelse(sex=="male",linyear,0)) vietnam.cr.glm7 <- update(vietnam.cr.glm6,.~.-linyear-sex:linyear+male.linyear) anova(vietnam.cr.glm7,vietnam.cr.glm6) ################################################### ### chunk number 93: vietnam.fitted ################################################### fval <- array(fitted(vietnam.cr.glm6),dim=c(5,2,3), dimnames= list(Year=levels(vietnam.cr$year), Sex=levels(vietnam.cr$sex)[c(2,1)], Policy=levels(vietnam.cr$policy))) round(ftable(fval,row.vars=c(2,1),col.vars=3),4) ################################################### ### chunk number 94: ################################################### fvalA <- fval[,,1] fvalB <- fval[,,2]*(1-fvalA) fvalC <- fval[,,3]*(1-fvalA-fvalB) fvalD <- (1-fvalA-fvalB-fvalC) fvalABCD <- array(c(fvalA,fvalB,fvalC,fvalD),dim=c(5,2,4)) dimnames(fvalABCD) <- list( Year=levels(vietnam$year), Sex=levels(vietnam$sex), Policy=levels(vietnam$policy)) round(ftable(fvalABCD,row.vars=c(2,1),col.vars=3),3) ################################################### ### chunk number 95: vietnam.females.glm ################################################### vietnam.females.glm <- glm(count~year+policy,data=vietnam, subset=sex=="female",family='poisson') vietnam.females.glm1 <- update(vietnam.females.glm,.~.+policy:liny) anova(vietnam.females.glm,vietnam.females.glm1) round(coef(summary(vietnam.females.glm1))[-c(1:5),],4) vietnam.females.glm2 <- update(vietnam.females.glm1,.~.-policy:liny+linp:liny) anova(vietnam.females.glm2,vietnam.females.glm1) round(coef(summary(vietnam.females.glm2))[-c(1:5),],4) vietnam.females.glm3 <- update(vietnam.females.glm2,.~.-linp:liny) anova(vietnam.females.glm3,vietnam.females.glm2) ################################################### ### chunk number 96: vietnam.cr.plot eval=FALSE ################################################### ## #vietnam$cr.fval <- as.vector(aperm(fvalABCD,c(2,3,1))) ## vietnam.cr.df <- as.data.frame.table(aperm(fvalABCD,c(2,1,3))) ## vietnam.cr.df$pobs <- as.vector(vietnam.ob.tab[c(2,1),,]) ## vietnam.cr.df$year <- with(vietnam.cr.df,ifelse(Year=="G",7,unclass(vietnam.cr.df$Year))) ## # ## print(xyplot(Freq~year|Sex,data=vietnam.cr.df, ## groups=Policy, ## panel=function(x,y,type,subscripts,...){ ## lpoints(vietnam.cr.df$year[subscripts], ## vietnam.cr.df$pobs[subscripts], ## pch=LETTERS[vietnam.cr.df$Policy[subscripts]],cex=1.5) ## panel.superpose(x,y,type="l",subscripts,...)})) ################################################### ### chunk number 97: ################################################### #vietnam$cr.fval <- as.vector(aperm(fvalABCD,c(2,3,1))) vietnam.cr.df <- as.data.frame.table(aperm(fvalABCD,c(2,1,3))) vietnam.cr.df$pobs <- as.vector(vietnam.ob.tab[c(2,1),,]) vietnam.cr.df$year <- with(vietnam.cr.df,ifelse(Year=="G",7,unclass(vietnam.cr.df$Year))) # print(xyplot(Freq~year|Sex,data=vietnam.cr.df, groups=Policy, panel=function(x,y,type,subscripts,...){ lpoints(vietnam.cr.df$year[subscripts], vietnam.cr.df$pobs[subscripts], pch=LETTERS[vietnam.cr.df$Policy[subscripts]],cex=1.5) panel.superpose(x,y,type="l",subscripts,...)})) ################################################### ### chunk number 98: ################################################### toxaemia <- read.csv("Data/toxaemia.csv") ################################################### ### chunk number 99: eval=FALSE ################################################### ## data(toxaemia,package="SMIR") ################################################### ### chunk number 100: toxaemia.prop.table ################################################### tox.prop.table1 <- with(toxaemia,prop.table(tapply(count, list(class=class,response=response,smoke=smoke),sum),c(1,3))[,c(2,1,4,3),1:2]) tox.prop.table2 <- with(toxaemia,prop.table(tapply(count, list(class=class,response=response,smoke=smoke),sum),c(1,3))[,c(2,1,4,3),3,drop=FALSE]) ################################################### ### chunk number 101: ################################################### round(ftable(tox.prop.table1,row.vars=c(1),col.vars=c(3,2)),3) round(ftable(tox.prop.table2,row.vars=c(1),col.vars=c(3,2)),3) ################################################### ### chunk number 102: hypertension ################################################### toxaemia$H <- factor(toxaemia$response%in%c("HN","HU")) round(ftable(with(toxaemia,prop.table(tapply(count,list(class,H,smoke),sum), margin=c(1,3)))[,2,]),4) ################################################### ### chunk number 103: hypertension.glm ################################################### toxaemia1 <- data.frame(xtabs(count~class+H+smoke,data=toxaemia)) hypertension <- cbind(toxaemia1[toxaemia1$H=="TRUE",-2],t=toxaemia1[toxaemia1$H=="FALSE",4]) hypertension.glm <- glm(cbind(Freq,t)~1,data=hypertension,family=binomial) #summary(hypertension.glm) hypertension.glm1 <- update(hypertension.glm,.~.+class) #anova(hypertension.glm,hypertension.glm1) hypertension.glm2 <- update(hypertension.glm1,.~.+smoke) #anova(hypertension.glm1,hypertension.glm2) hypertension.glm3 <- update(hypertension.glm2,.~.+class*smoke) #anova(hypertension.glm2,hypertension.glm3) anova(hypertension.glm,hypertension.glm1,hypertension.glm2,hypertension.glm3) round(coef(summary(hypertension.glm3)),4) ################################################### ### chunk number 104: hypertension.glm4 ################################################### hypertension <- transform(hypertension, sn1=ifelse(smoke=="0",0,1), cn1=ifelse(class=="I",0,1), s3=ifelse(smoke=="20+",1,0)) hypertension.glm4 <- update(hypertension.glm,.~ class + smoke + sn1 + cn1 + class:sn1+cn1:s3) summary(hypertension.glm4) ################################################### ### chunk number 105: hypertension.glm5 ################################################### hypertension <- transform(hypertension, c2 = ifelse(class=="II",1,0)) hypertension.glm5 <- update(hypertension.glm4, .~.-class:sn1 + cn1:sn1 + c2:sn1) coef(summary(hypertension.glm5)) ################################################### ### chunk number 106: hypertension.glm6 ################################################### hypertension.glm6 <- update(hypertension.glm5, .~ . -smoke) summary(hypertension.glm6) hypertension.glm7 <- update(hypertension.glm6, .~ . - cn1:s3) summary(hypertension.glm7) hypertension.glm8 <- update(hypertension.glm7, .~ . -class) summary(hypertension.glm8) hypertension.glm9 <- update(hypertension.glm8, .~ . -sn1) summary(hypertension.glm9) hypertension.glm10 <- update(hypertension.glm9, .~ . -cn1) summary(hypertension.glm10) hypertension$fval <- fitted(hypertension.glm10) print(with(hypertension,tapply(fval,list(Class=class,Smoke=smoke),mean)),digits=3) ################################################### ### chunk number 107: proteinurea.response ################################################### toxaemia$U <- factor(toxaemia$response%in%c("NU","HU")) round(ftable(with(toxaemia,prop.table(tapply(count, list(class=class,hypertension=U,smoke=smoke),sum),c(1,3))), row.vars=1,col.vars=c(3,2)),4)[,c(2,4,6)] ################################################### ### chunk number 108: proteinurea.glm ################################################### proteinurea <- data.frame(xtabs(count~class+U+smoke,data=toxaemia)) proteinurea2 <- cbind(proteinurea[proteinurea$U=="TRUE",-2],t=proteinurea[proteinurea$U=="FALSE",4]) proteinurea.glm <- glm(cbind(Freq,t)~class+smoke,data=proteinurea2,family=binomial) anova(proteinurea.glm) round(coef(summary(proteinurea.glm)),4) ################################################### ### chunk number 109: proteinurea.glm1 ################################################### proteinurea.glm1 <- update(proteinurea.glm, .~.-smoke) summary(proteinurea.glm1) ################################################### ### chunk number 110: proteinurea.glm2 ################################################### proteinurea2 <- transform(proteinurea2, c2 = ifelse(class=="II",1,0), c3 = ifelse(class=="III",1,0)) proteinurea.glm2 <- update(proteinurea.glm1, .~. -class + c2 + c3) summary(proteinurea.glm2) proteinurea2$fval <- fitted(proteinurea.glm2) round(with(proteinurea2,tapply(fval,class,mean)),3) ################################################### ### chunk number 111: correlated.outcomes ################################################### table1 <- with(toxaemia,tapply(count,list(class,smoke,response),sum)) lodds <- log(table1[,,1]*table1[,,4]/(table1[,,2]*table1[,,3])) round(-lodds,2) ################################################### ### chunk number 112: toxaemia.glm11 ################################################### toxaemia.glm11 <- glm(count~class:smoke+H*U,data=toxaemia,family=poisson) print(anova(toxaemia.glm11),digits=7) coef(summary(toxaemia.glm11))[c(2,3,18),] ################################################### ### chunk number 113: toxaemia.glm14 ################################################### toxaemia.glm12 <- update(toxaemia.glm11,.~H*U+class:smoke+ H:class+H:smoke+H:class:smoke+ U:class+U:smoke+U:class:smoke+H:U:class+H:U:smoke) toxaemia.glm13 <- update(toxaemia.glm12,.~.-H:U:class) toxaemia.glm14 <- update(toxaemia.glm13,.~.-H:U:smoke) anova(toxaemia.glm12,toxaemia.glm13,toxaemia.glm14) ################################################### ### chunk number 114: toxaemia.glm13 ################################################### round(coef(summary(toxaemia.glm13))[47:48,],4) ################################################### ### chunk number 115: toxaemia.glm15 ################################################### toxaemia$lsmoke <- unclass(toxaemia$smoke) toxaemia.glm15 <- update(toxaemia.glm13,.~.-H:U:smoke+H:U:lsmoke) anova(toxaemia.glm13,toxaemia.glm15) coef(summary(toxaemia.glm15))[c(4,47),] ################################################### ### chunk number 116: toxaemia.glm17 ################################################### toxaemia.glm16 <- update(toxaemia.glm15,.~.-U:class:smoke) toxaemia.glm17 <- update(toxaemia.glm16,.~.-U:smoke) anova(toxaemia.glm15,toxaemia.glm16,toxaemia.glm17) ################################################### ### chunk number 117: toxaemia.glm18 ################################################### toxaemia.glm18 <- update(toxaemia.glm17,.~.+lsmoke+H:lsmoke+U:lsmoke) round(coef(summary(toxaemia.glm18))[23:24,],4) ################################################### ### chunk number 118: toxaemia.glm19 ################################################### toxaemia.glm19 <- update(toxaemia.glm18, .~. -H:smoke) anova(toxaemia.glm18,toxaemia.glm19) round(coef(summary(toxaemia.glm19))[29:37,],4) ################################################### ### chunk number 119: toxaemia.glm20 ################################################### toxaemia <- transform(toxaemia, cn1=ifelse(class!="I",1,0), sn1=ifelse(smoke!="0",1,0), s3=ifelse(smoke=="20+",1,0)) toxaemia.glm20 <- update(toxaemia.glm19,.~.-H:class:smoke+cn1+sn1+s3+class:sn1+ cn1:s3+H:(cn1+sn1+class:sn1+cn1:s3)) anova(toxaemia.glm19,toxaemia.glm20) round(coef(summary(toxaemia.glm20))[35,,drop=F],4) toxaemia.glm21 <- update(toxaemia.glm20,.~.-H:cn1:s3) anova(toxaemia.glm20,toxaemia.glm21) round(coef(summary(toxaemia.glm21))[31:34,],4) toxaemia <- transform(toxaemia,c2 =ifelse(class=="II",1,0)) toxaemia.glm22 <- update(toxaemia.glm21,.~.-H:class:sn1+H:cn1+H:c2+cn1:sn1+ c2:sn1+H:cn1:sn1+H:c2:sn1) anova(toxaemia.glm21,toxaemia.glm22) round(coef(summary(toxaemia.glm22))[23:26,],4) toxaemia <- transform(toxaemia,c3=ifelse(class=="III",1,0)) toxaemia.glm23 <- update(toxaemia.glm22,.~.-U:class++c3+U:(c2+c3)) anova(toxaemia.glm22,toxaemia.glm23) round(coef(summary(toxaemia.glm23))[19:22,],4) toxaemia.glm24 <- update(toxaemia.glm23,.~.-H:class+H:cn1) anova(toxaemia.glm23,toxaemia.glm24) round(coef(summary(toxaemia.glm24))[c(2,3,10,19:28),],4) ################################################### ### chunk number 120: H.ptab ################################################### H.tab <- with(toxaemia,tapply(fitted(toxaemia.glm24),list(Class=class,Smoke=smoke,H),mean)) H.ptab <- H.tab[,,2]/(H.tab[,,1]+H.tab[,,2]) round(H.ptab,3) ################################################### ### chunk number 121: U.ptab ################################################### U.tab <- with(toxaemia,tapply(fitted(toxaemia.glm24),list(Class=class,U),mean)) U.ptab <- U.tab[,2]/(U.tab[,1]+U.tab[,2]) round(U.ptab,3)