budwormOutput - > budworm<‐...

Info icon This 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

{[ snackBarMessage ]}

What students are saying

  • Left Quote Icon

    As a current student on this bumpy collegiate pathway, I stumbled upon Course Hero, where I can find study resources for nearly all my courses, get online help from tutors 24/7, and even share my old projects, papers, and lecture notes with other students.

    Student Picture

    Kiran Temple University Fox School of Business ‘17, Course Hero Intern

  • Left Quote Icon

    I cannot even describe how much Course Hero helped me this summer. It’s truly become something I can always rely on and help me. In the end, I was not only able to survive summer classes, but I was able to thrive thanks to Course Hero.

    Student Picture

    Dana University of Pennsylvania ‘17, Course Hero Intern

  • Left Quote Icon

    The ability to access any university’s resources through Course Hero proved invaluable in my case. I was behind on Tulane coursework and actually used UCLA’s materials to help me move forward and get everything together on time.

    Student Picture

    Jill Tulane University ‘16, Course Hero Intern