hwk4-1 - > library(MASS>> 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: > library(MASS) > > budworm <‐ read.table(file="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 > budworm$ldose <‐ log2(budworm$dose) > male <‐ budworm$sex==1 > budworm$male <‐ as.numeric(male) > budworm$female <‐ as.numeric(!male) > budworm$maleldose <‐ budworm$male*budworm$ldose > budworm$femaleldose <‐ budworm$female*budworm$ldose > > > m1 <‐ glm(p~male+maleldose+female+femaleldose‐1, + data=budworm,weights=rep(20,12), + family=binomial(link="logit")) > summary(m1) Call: glm(formula = p ~ male + maleldose + female + femaleldose ‐ 1, 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|) male ‐2.8186 0.5480 ‐5.143 2.70e‐07 *** maleldose 1.2589 0.2121 5.937 2.91e‐09 *** female ‐2.9935 0.5527 ‐5.416 6.09e‐08 *** femaleldose 0.9060 0.1671 5.422 5.89e‐08 *** ‐‐‐ Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 126.2269 on 12 degrees of freedom Residual deviance: 4.9937 on 8 degrees of freedom AIC: 43.104 Number of Fisher Scoring iterations: 4 > > 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) + } > logit <‐ function(x){ log(x/(1‐x)) } > varbetahat <‐ summary(m1)$cov.scaled > > #MALES > > #Fieller intervals > p <‐ .9 > z <‐ qnorm(1‐.05/2,0,1) > coefs <‐ coef(m1)[1:2] > v <‐ varbetahat[1:2,1:2] > > a <‐ coefs[2]^2‐(z^2)*v[2,2] > b <‐ 2*coefs[2]*coefs[1]‐2*(z^2)*v[1,2] ‐2*logit(p)*coefs[2] > c <‐ coefs[1]^2‐(z^2)*v[1,1]+logit(p)^2‐2*logit(p)*coefs[1] > lLD90interval.males.Fieller <‐ quad.form(a,b,c) > LD90interval.males.Fieller <‐ 2^lLD90interval.males.Fieller > lLD90interval.males.Fieller maleldose maleldose 3.432974 4.951387 > LD90interval.males.Fieller maleldose maleldose 10.80011 30.93970 > > #Wald intervals > lLD90 <‐ (logit(.9)‐coefs[1])/coefs[2] > se.lLD90 <‐ sqrt((coefs[2]^(‐2))* t(c(1,lLD90))%*%v%*%c(1,lLD90)) > lLD90 male 3.984099 > se.lLD90 [,1] [1,] 0.3516367 > Lower.LD90 <‐ 2^( lLD90‐1.96*se.lLD90 ) > Upper.LD90 <‐ 2^( lLD90+1.96*se.lLD90 ) > LD90interval.males.Wald <‐ cbind(Lower.LD90,2^(lLD90),Upper.LD90) > lLD90interval.males.Wald <‐ log2(LD90interval.males.Wald) > lLD90interval.males.Wald [,1] [,2] [,3] male 3.294891 3.984099 4.673307 > LD90interval.males.Wald [,1] [,2] [,3] male 9.81434 15.82462 25.51559 > dose.p(m1,cf=c(1,2),p=c(.9)) Dose SE p = 0.9: 3.984099 0.3516367 > > # Profile likelihood intervals > pl <‐ c(0)[‐1] > pl numeric(0) > > x0vec <‐ seq(from=3, to=5, by=.001)#from Wald & Fieller approach we know > #that the limits should be in this range > n <‐ length(x0vec) > for(i in 1:n){ + bud <‐ budworm + bud$off <‐ budworm$male*log(.9/.1) + bud$x1 <‐ (budworm$ldose‐x0vec[i])*budworm$male + modfit <‐ glm(p ~ ‐1+x1+female+femaleldose,data=bud,weights=rep(20,12), + family=binomial(link="logit"),offset=off) + pl <‐ c(pl,as.numeric(logLik(modfit))) + } > > plot(x0vec,pl,type="l", main="Profile loglikelihood for x0 (males)") > abline(h=max(pl)‐qchisq(.95,1)/2) > # MLE of x0: > x0vec[pl==max(pl)] [1] 3.984 > # MLE of z0: > 2^x0vec[pl==max(pl)] [1] 15.82353 > > max(pl)‐qchisq(.95,1)/2 [1] ‐19.47279 > mat <‐ cbind(2^x0vec,x0vec,pl) > mat[abs(pl‐(max(pl)‐qchisq(.95,1)/2))<.01,] x0vec pl [1,] 10.47588 3.389 ‐19.48122 [2,] 10.48315 3.390 ‐19.47377 [3,] 10.49042 3.391 ‐19.46633 [4,] 28.30529 4.823 ‐19.46607 [5,] 28.32492 4.824 ‐19.46976 [6,] 28.34456 4.825 ‐19.47346 [7,] 28.36421 4.826 ‐19.47717 [8,] 28.38388 4.827 ‐19.48087 > # limits for log2dose are 3.390 and 4.825 > # for dose they are 10.48 and 28.34 > > > > #FEMALES > > coefs <‐ coef(m1)[3:4] > v <‐ varbetahat[3:4,3:4] > > a <‐ coefs[2]^2‐(z^2)*v[2,2] > b <‐ 2*coefs[2]*coefs[1]‐2*(z^2)*v[1,2] ‐2*logit(p)*coefs[2] > c <‐ coefs[1]^2‐(z^2)*v[1,1]+logit(p)^2‐2*logit(p)*coefs[1] > lLD90interval.females.Fieller <‐ quad.form(a,b,c) > LD90interval.females.Fieller <‐ 2^lLD90interval.females.Fieller > lLD90interval.females.Fieller femaleldose femaleldose 4.886181 7.390556 > LD90interval.females.Fieller femaleldose femaleldose 29.57244 167.79505 > > #Wald intervals > lLD90 <‐ (logit(.9)‐coefs[1])/coefs[2] > se.lLD90 <‐ sqrt((coefs[2]^(‐2))* t(c(1,lLD90))%*%v%*%c(1,lLD90)) > lLD90 female 5.729092 > se.lLD90 [,1] [1,] 0.5629648 > Lower.LD90 <‐ 2^( lLD90‐1.96*se.lLD90 ) > Upper.LD90 <‐ 2^( lLD90+1.96*se.lLD90 ) > LD90interval.females.Wald <‐ cbind(Lower.LD90,2^(lLD90),Upper.LD90) > lLD90interval.females.Wald <‐ log2(LD90interval.females.Wald) > lLD90interval.females.Wald [,1] [,2] [,3] female 4.625681 5.729092 6.832503 > LD90interval.females.Wald [,1] [,2] [,3] female 24.68703 53.04307 113.9694 > dose.p(m1,cf=c(3,4),p=c(.9)) Dose SE p = 0.9: 5.729092 0.5629648 > > > # Profile likelihood intervals > pl <‐ c(0)[‐1] > pl numeric(0) > > x0vec <‐ seq(from=4.5, to=7.5, by=.001)#from Wald & Fielle r approach we know > #that the limits should be in this range e > n <‐ length(x0vec) > for(i in 1:n){ + bud <‐ budworm + bud$off <‐ budworm$female*log(.9/.1) + bud$x1 <‐ (budworm$ldose‐x0vec[i])*budworm$female + modfit <‐ glm(p ~ ‐1+x1+male+maleldose,data=bud,weigh ts=rep(20,12), 2 + family=binomial(link="logit"),offset=off) + pl <‐ c(pl,as.numeric(logLik(modfit))) + } > > plot(x0vec,pl,type="l", main="Profile loglikelihood for x0 (fe males)") > abline(h=max(pl)‐qchisq(.95,1)/2) > # MLE of x0: > x0vec[pl==max(pl)] [1] 5.729 > # MLE of z0: > 2^x0vec[pl==max(pl)] [1] 53.03967 > > max(pl)‐qchisq(.95,1)/2 [1] ‐19.47279 > mat <‐ cbind(2^x0vec,x0vec,pl) > mat[abs(pl‐(max(pl)‐qchisq(.95,1)/2))<.01,] x0vec pl [1,] 28.72032 4.844 ‐19.48143 [2,] 28.74024 4.845 ‐19.47605 [3,] 28.76016 4.846 ‐19.47069 [4,] 28.78011 4.847 ‐19.46533 [5,] 151.16706 7.240 ‐19.46423 [6,] 151.27188 7.241 ‐19.46605 [7,] 151.37677 7.242 ‐19.46788 [8,] 151.48173 7.243 ‐19.46970 [9,] 151.58677 7.244 ‐19.47153 [10,] 151.69187 7.245 ‐19.47336 [11,] 151.79706 7.246 ‐19.47518 [12,] 151.90231 7.247 ‐19.47701 [13,] 152.00764 7.248 ‐19.47884 [14,] 152.11304 7.249 ‐19.48066 [15,] 152.21851 7.250 ‐19.48249 > # limits for log2dose are 4.846 and 7.245 > # for dose they are 28.76 and 151.69 > ...
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