budwormOutput - > budworm <‐

Info iconThis preview shows page 1. Sign up to view the full content.

View Full Document Right Arrow Icon
This is the end of the preview. Sign up to access the rest of the document.

Unformatted text preview: > budworm <‐ read.table(file="N:\\courses\\stat8620\\Fall 08\\budworm.dat",header=T) > #budworm <‐ read.table(file="C:\\Documents and Settings\\dhall\\My Documents\\Dan's Work Stuff\\courses\\STAT8620\\Fall 08\\budworm.dat",header=T) > budworm[1:3,] sex dose y m 1 1 1 1 20 2 1 2 4 20 3 1 4 9 20 > budworm$p <‐ budworm$y/budworm$m > > male <‐ budworm$sex==1 > budworm$male <‐ as.numeric(male) > plot(budworm$dose[male],budworm$p[male], + main="Sample prop dead vs. dose",pch=1, sub="Plot (a)") > points(budworm$dose[!male],budworm$p[!male],pch=2) > legend(locator(1),pch=c(1,2),legend=c("male","female")) > > budworm$samplogit <‐ log((budworm$y+0.5)/(budworm$m‐budworm$y+0.5)) > budworm$ldose <‐ log2(budworm$dose) > > plot(budworm$dose[male],budworm$samplogit[male], + main="Sample logit vs. dose",pch=1, sub="Plot (b)") > points(budworm$dose[!male],budworm$samplogit[!male],pch=2) > legend(locator(1),pch=c(1,2),legend=c("male","female")) > > plot(budworm$ldose[male],budworm$samplogit[male], + main="Sample logit vs. log2(dose)",pch=1, sub="Plot (c)") > points(budworm$ldose[!male],budworm$samplogit[!male],pch=2) > legend(locator(1),pch=c(1,2),legend=c("male","female")) > > > m1 <‐ glm(p~ldose*male,data=budworm,weights=rep(20,12), + family=binomial(link="logit")) > dead.alive <‐ cbind(budworm$y,rep(20,12)‐budworm$y) > m1alt <‐ glm(dead.alive~ldose*male,data=budworm, + family=binomial(link="logit")) > > bud2 <‐ rbind(budworm,budworm) > bud2$dead <‐ c(rep(1,12),rep(0,12)) > bud2$freq <‐ c(budworm$y,20‐budworm$y) > dead.alt <‐ rep(bud2$dead,times=bud2$freq) > ldose.alt <‐ rep(bud2$ldose,times=bud2$freq) > male.alt <‐ rep(bud2$male,times=bud2$freq) > budworm2 <‐ data.frame(cbind(dead.alt,ldose.alt,male.alt)) > budworm2[1:6,] dead.alt ldose.alt male.alt 1 1 0 1 2 1 1 1 3 1 1 1 4 1 1 1 5 1 1 1 6 1 2 1 > budworm2[121:126,] dead.alt ldose.alt male.alt 121 0 0 1 122 0 0 1 123 0 0 1 124 0 0 1 125 0 0 1 126 0 0 1 > m1alt2 <‐ glm(dead.alt~ldose.alt*male.alt,data=budworm2, + family=binomial(link="logit")) > > m1 Call: glm(formula = p ~ ldose * male, family = binomial(link = "logit"), data = budworm, weights = rep(20, 12)) Coefficients: (Intercept) ldose male ldose:male ‐2.9935 0.9060 0.1750 0.3529 Degrees of Freedom: 11 Total (i.e. Null); 8 Residual Null Deviance: 124.9 Residual Deviance: 4.994 AIC: 43.1 > m1alt Call: glm(formula = dead.alive ~ ldose * male, family = binomial(link = "logit"), data = budworm) Coefficients: (Intercept) ldose male ldose:male ‐2.9935 0.9060 0.1750 0.3529 Degrees of Freedom: 11 Total (i.e. Null); 8 Residual Null Deviance: 124.9 Residual Deviance: 4.994 AIC: 43.1 > m1alt2 Call: glm(formula = dead.alt ~ ldose.alt * male.alt, family = binomial(link = "logit"), data = budworm2) Coefficients: (Intercept) ldose.alt male.alt ldose.alt:male.alt ‐2.9935 0.9060 0.1750 0.3529 Degrees of Freedom: 239 Total (i.e. Null); 236 Residual Null Deviance: 331.4 Residual Deviance: 211.5 AIC: 219.5 > logLik(m1) 'log Lik.' ‐17.55206 (df=4) > logLik(m1alt) 'log Lik.' ‐17.55206 (df=4) > logLik(m1alt2) 'log Lik.' ‐105.7388 (df=4) > logLik(m1)‐sum(log( choose(budworm$m,budworm$y) ) ) 'log Lik.' ‐105.7388 (df=4) > > summary(m1,corr=F) Call: glm(formula = p ~ ldose * male, family = binomial(link = "logit"), data = budworm, weights = rep(20, 12)) Deviance Residuals: Min 1Q Median 3Q Max ‐1.39849 ‐0.32094 ‐0.07592 0.38220 1.10375 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) ‐2.9935 0.5527 ‐5.416 6.09e‐08 *** ldose 0.9060 0.1671 5.422 5.89e‐08 *** male 0.1750 0.7783 0.225 0.822 ldose:male 0.3529 0.2700 1.307 0.191 ‐‐‐ Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 124.8756 on 11 degrees of freedom Residual deviance: 4.9937 on 8 degrees of freedom AIC: 43.104 Number of Fisher Scoring iterations: 4 > > #saturated model: > m2 <‐ glm(p~factor(ldose)*factor(male),data=budworm,weights=rep(20,12), + family=binomial(link="logit")) > > summary(m2,corr=F) Call: glm(formula = p ~ factor(ldose) * factor(male), family = binomial(link = "logit"), data = budworm, weights = rep(20, 12)) Deviance Residuals: [1] 0 0 0 0 0 0 0 0 0 0 0 0 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) ‐25.752 52998.328 ‐4.86e‐04 1 factor(ldose)1 23.555 52998.328 4.44e‐04 1 factor(ldose)2 24.904 52998.328 4.70e‐04 1 factor(ldose)3 25.752 52998.328 4.86e‐04 1 factor(ldose)4 26.157 52998.328 4.94e‐04 1 factor(ldose)5 27.138 52998.328 0.001 1 factor(male)1 22.807 52998.328 4.30e‐04 1 factor(ldose)1:factor(male)1 ‐21.996 52998.328 ‐4.15e‐04 1 factor(ldose)2:factor(male)1 ‐22.161 52998.328 ‐4.18e‐04 1 factor(ldose)3:factor(male)1 ‐22.188 52998.328 ‐4.19e‐04 1 factor(ldose)4:factor(male)1 ‐21.016 52998.328 ‐3.97e‐04 1 factor(ldose)5:factor(male)1 1.558 74951.009 2.08e‐05 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 1.2488e+02 on 11 degrees of freedom Residual deviance: 5.2388e‐10 on 0 degrees of freedom AIC: 54.11 Number of Fisher Scoring iterations: 22 > > #deviance of model m1 is GOF statistic: > summary(m1) Call: glm(formula = p ~ ldose * male, family = binomial(link = "logit"), data = budworm, weights = rep(20, 12)) Deviance Residuals: Min 1Q Median 3Q Max ‐1.39849 ‐0.32094 ‐0.07592 0.38220 1.10375 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) ‐2.9935 0.5527 ‐5.416 6.09e‐08 *** ldose 0.9060 0.1671 5.422 5.89e‐08 *** male 0.1750 0.7783 0.225 0.822 ldose:male 0.3529 0.2700 1.307 0.191 ‐‐‐ Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 124.8756 on 11 degrees of freedom Residual deviance: 4.9937 on 8 degrees of freedom AIC: 43.104 Number of Fisher Scoring iterations: 4 > deviance(m1) [1] 4.993727 > 2*(logLik(m2)‐logLik(m1)) [1] 4.993727 attr(,"df") [1] 12 attr(,"class") [1] "logLik" > #Pearson X^2 statistic: > sum(resid(m1,type="pearson")^2) [1] 3.504694 > > # The value commented out below reproduces the loglikelihood result given > # for model m2 in budworm.sas > #logLik(m2)‐sum(log( choose(budworm$m,budworm$y) ) ) > > ldose0 <‐ seq(from=0,to=5,by=.1) > expit <‐ function(x) {1/(1+exp(‐x)) } > plot(budworm$ldose[male],budworm$p[male],xlab="Log2(dose)", + ylab="Mortality Probability", sub="Plot (d)", + main="Fitted probability from model m1",pch=1) > points(budworm$ldose[!male],budworm$p[!male],pch=2) > lines(ldose0, + expit( (coef(m1)[1]+coef(m1)[3])+(coef(m1)[2]+coef(m1)[4])*ldose0),lty=2) > lines(ldose0, + expit( (coef(m1)[1])+(coef(m1)[2])*ldose0),lty=1) > legend(locator(1),legend=c("males","females"),lty=c(2,1),pch=c(1,2)) > > > plot(budworm$ldose[male],budworm$samplogit[male],xlab="Log2(dose)", + ylab="Log Odds", + main="Fitted log odds from model m1",pch=1, sub="Plot (e)") > points(budworm$ldose[!male],budworm$samplogit[!male],pch=2) > lines(ldose0, + ( (coef(m1)[1]+coef(m1)[3])+(coef(m1)[2]+coef(m1)[4])*ldose0),lty=2) > lines(ldose0, + ( (coef(m1)[1])+(coef(m1)[2])*ldose0),lty=1) > legend(locator(1),legend=c("males","females"),lty=c(2,1),pch=c(1,2)) > > > budworm$ldosec <‐ budworm$ldose‐2.5 > m3 <‐ glm(p~ldosec*male,data=budworm,weights=rep(20,12), + family=binomial(link="logit")) > > summary(m3) Call: glm(formula = p ~ ldosec * male, family = binomial(link = "logit"), data = budworm, weights = rep(20, 12)) Deviance Residuals: Min 1Q Median 3Q Max ‐1.39849 ‐0.32094 ‐0.07592 0.38220 1.10375 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) ‐0.7285 0.2455 ‐2.967 0.00301 ** ldosec 0.9060 0.1671 5.422 5.89e‐08 *** male 1.0573 0.3581 2.952 0.00316 ** ldosec:male 0.3529 0.2700 1.307 0.19117 ‐‐‐ Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 124.8756 on 11 degrees of freedom Residual deviance: 4.9937 on 8 degrees of freedom AIC: 43.104 Number of Fisher Scoring iterations: 4 > > m4 <‐ glm(p~male+ldose,data=budworm,weights=rep(20,12), + family=binomial(link="logit")) > anova(m4,m1,test="Chisq") Analysis of Deviance Table Model 1: p ~ male + ldose Model 2: p ~ ldose * male Resid. Df Resid. Dev Df Deviance P(>|Chi|) 1 9 6.7571 2 8 4.9937 1 1.7633 0.1842 > anova(m4,m3,test="Chisq") Analysis of Deviance Table Model 1: p ~ male + ldose Model 2: p ~ ldosec * male Resid. Df Resid. Dev Df Deviance P(>|Chi|) 1 9 6.7571 2 8 4.9937 1 1.7633 0.1842 > > summary(m4) Call: glm(formula = p ~ male + ldose, family = binomial(link = "logit"), data = budworm, weights = rep(20, 12)) Deviance Residuals: Min 1Q Median 3Q Max ‐1.10540 ‐0.65343 ‐0.02225 0.48471 1.42944 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) ‐3.4732 0.4685 ‐7.413 1.23e‐13 *** male 1.1007 0.3558 3.093 0.00198 ** ldose 1.0642 0.1311 8.119 4.70e‐16 *** ‐‐‐ Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 124.876 on 11 degrees of freedom Residual deviance: 6.757 on 9 degrees of freedom AIC: 42.867 Number of Fisher Scoring iterations: 4 > > m5 <‐ glm(p~factor(male)+ldose‐1,data=budworm,weights=rep(20,12), + family=binomial(link="logit")) > summary(m5) Call: glm(formula = p ~ factor(male) + ldose ‐ 1, family = binomial(link = "logit"), data = budworm, weights = rep(20, 12)) Deviance Residuals: Min 1Q Median 3Q Max ‐1.10540 ‐0.65343 ‐0.02225 0.48471 1.42944 Coefficients: Estimate Std. Error z value Pr(>|z|) factor(male)0 ‐3.4732 0.4685 ‐7.413 1.23e‐13 *** factor(male)1 ‐2.3724 0.3855 ‐6.154 7.56e‐10 *** ldose 1.0642 0.1311 8.119 4.70e‐16 *** ‐‐‐ Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 126.227 on 12 degrees of freedom Residual deviance: 6.757 on 9 degrees of freedom AIC: 42.867 Number of Fisher Scoring iterations: 4 > > coef(m4) (Intercept) male ldose ‐3.473155 1.100743 1.064214 > coef(m5) factor(male)0 factor(male)1 ldose ‐3.473155 ‐2.372412 1.064214 > > > #LD50 and Wald confidence interval: females > # if natural log of dose hade been used instead of log2(dose), then > # use exp() instead of 2^() below to invert the logarithm > varbetahat <‐ summary(m5)$cov.scaled > > lLD50 <‐ (0‐coef(m5)[1])/coef(m5)[3] > se.lLD50 <‐ sqrt((coef(m5)[3]^(‐2))* t(c(1,0,lLD50))%*%varbetahat%*%c(1,0,lLD50)) > lLD50 factor(male)0 3.263587 > se.lLD50 [,1] [1,] 0.2297539 > > #The same result can be obtained using the dose.p function in the MASS library > library(MASS) > dose.p(m5,cf=c(1,3),p=c(.5,.9)) Dose SE p = 0.5: 3.263587 0.2297539 p = 0.9: 5.328233 0.3611835 > > Lower.LD50 <‐ 2^( lLD50‐1.96*se.lLD50 ) > Upper.LD50 <‐ 2^( lLD50+1.96*se.lLD50 ) > LD50interval.females.Wald <‐ cbind(Lower.LD50,2^(lLD50),Upper.LD50) > LD50interval.females.Wald [,1] [,2] [,3] factor(male)0 7.028758 9.60368 13.12190 > > # if you wanted the LD90, you'd do the following > # lLD90 <‐ (logit(.9)‐coef(m5)[1])/coef(m5)[3] > # se.lLD90 <‐ sqrt((coef(m5)[3]^(‐2))* t(c(1,0,lLD90))%*%varbetahat%*%c(1,0,lLD90)) > > #LD50 and Wald confidence interval: males > dose.p(m5,cf=c(2,3),p=c(.5,.9)) Dose SE p = 0.5: 2.229262 0.2259649 p = 0.9: 4.293908 0.3336275 > > lLD50 <‐ (0‐coef(m5)[2])/coef(m5)[3] > se.lLD50 <‐ sqrt((coef(m5)[3]^(‐2))* t(c(0,1,lLD50))%*%varbetahat%*%c(0,1,lLD50)) > > Lower.LD50 <‐ 2^( lLD50‐1.96*se.lLD50 ) > Upper.LD50 <‐ 2^( lLD50+1.96*se.lLD50 ) > LD50interval.males.Wald <‐ cbind(Lower.LD50,2^(lLD50),Upper.LD50) > LD50interval.males.Wald [,1] [,2] [,3] factor(male)1 3.449461 4.688941 6.373798 > > > quad.form <‐ function(a,b,c){ + l <‐ (‐b ‐ sqrt(b^2‐4*a*c))/(2*a) + u <‐ (‐b + sqrt(b^2‐4*a*c))/(2*a) + c(l,u) + } > z <‐ qnorm(1‐.05/2,0,1) > a <‐ coef(m5)[3]^2‐(z^2)*varbetahat[3,3] > b <‐ 2*coef(m5)[3]*coef(m5)[2]‐2*(z^2)*varbetahat[2,3] > c <‐ coef(m5)[2]^2‐(z^2)*varbetahat[2,2] > lLD50interval.males.Fieller <‐ quad.form(a,b,c) > LD50interval.males.Fieller <‐ 2^lLD50interval.males.Fieller > lLD50interval.males.Fieller ldose ldose 1.768493 2.681297 > LD50interval.males.Fieller ldose ldose 3.406979 6.414322 > > ...
View Full Document

This note was uploaded on 11/13/2011 for the course STAT 8620 taught by Professor Hall during the Fall '11 term at University of Florida.

Ask a homework question - tutors are online