ESR.R.out - > ESR.R>> library(arm>...

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: > # ESR.R > > library(arm) > library(Design) > library(car) > > options(digits=5) > > ESR <‐ read.table(file="ESR.dat",header=T) > ESR[1:3,] individual fib glob y 1 1 2.52 38 0 2 2 2.56 31 0 3 3 2.19 33 0 > > m0 <‐ glm(y~1,data=ESR,family=binomial(link="logit")) > m1 <‐ glm(y~fib,data=ESR,family=binomial(link="logit")) > m2 <‐ glm(y~glob,data=ESR,family=binomial(link="logit")) > m3 <‐ glm(y~ fib + glob,data=ESR,family=binomial(link="logit")) > > anova(m0,m1) Analysis of Deviance Table Model 1: y ~ 1 Model 2: y ~ fib Resid. Df Resid. Dev Df Deviance 1 31 30.9 2 30 24.8 1 6.04 > anova(m0,m2) Analysis of Deviance Table Model 1: y ~ 1 Model 2: y ~ glob Resid. Df Resid. Dev Df Deviance 1 31 30.9 2 30 28.9 1 1.94 > anova(m0,m3) Analysis of Deviance Table Model 1: y ~ 1 Model 2: y ~ fib + glob Resid. Df Resid. Dev Df Deviance 1 31 30.9 2 29 23.0 2 7.91 > anova(m3) Analysis of Deviance Table Model: binomial, link: logit Response: y Terms added sequentially (first to last) Df Deviance Resid. Df Resid. Dev NULL 31 30.9 fib 1 6.04 30 24.8 glob 1 1.87 29 23.0 > > coef(m3) (Intercept) fib glob ‐12.79208 1.91037 0.15578 > coef(glm(y~ I(10*fib) + glob,data=ESR,family=binomial(link="logit"))) (Intercept) I(10 * fib) glob ‐12.79208 0.19104 0.15578 > summary(glm(y~ I(10*fib) + glob,data=ESR,family=binomial(link="logit"))) Call: glm(formula = y ~ I(10 * fib) + glob, family = binomial(link = "logit"), data = ESR) Deviance Residuals: Min 1Q Median 3Q Max ‐0.968 ‐0.612 ‐0.346 ‐0.212 2.264 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) ‐12.7921 5.7963 ‐2.21 0.027 * I(10 * fib) 0.1910 0.0971 1.97 0.049 * glob 0.1558 0.1195 1.30 0.193 ‐‐‐ Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 30.885 on 31 degrees of freedom Residual deviance: 22.971 on 29 degrees of freedom AIC: 28.97 Number of Fisher Scoring iterations: 5 > > f0 <‐ seq(from=min(ESR$fib),to=max(ESR$fib),length=100) > expit <‐ function(x) {1/(1+exp(‐x)) } > > quantile(ESR$fib) 0% 25% 50% 75% 100% 2.0900 2.2900 2.6000 3.1675 5.0600 > quantile(ESR$glob) 0% 25% 50% 75% 100% 28.00 31.75 36.00 38.00 46.00 > > plot(ESR$fib,ESR$y,type="p",xlab="Level of fibrinogen",ylab="y = I( ESR >= 20 )", + main="Fitted probability from model m3") > lines(f0,expit( coef(m3)[1]+coef(m3)[2]*f0+coef(m3)[3]*32),lty=2) > lines(f0,expit( coef(m3)[1]+coef(m3)[2]*f0+coef(m3)[3]*36),lty=1) > lines(f0,expit( coef(m3)[1]+coef(m3)[2]*f0+coef(m3)[3]*37),lty=3) > legend(locator(1),legend=c("glob=32 (1st quartile)","glob=36 (median)", + "glob=38 (3rd quartile)"),lty=c(2,1,3)) > > > # effect of .1 increase in fib at the median of fib, median of glob > expit( coef(m3)[1]+coef(m3)[2]*(2.6+.1)+coef(m3)[3]*36 )‐ + expit( coef(m3)[1]+coef(m3)[2]*(2.6)+coef(m3)[3]*36 ) (Intercept) 0.018268 > > # effect of .1 increase in fib at the 1st quartile of fib, median of glob > expit( coef(m3)[1]+coef(m3)[2]*(2.3+.1)+coef(m3)[3]*36 )‐ + expit( coef(m3)[1]+coef(m3)[2]*(2.3)+coef(m3)[3]*36 ) (Intercept) 0.011336 > > # effect of .1 increase in fib at the 3rd quartile of fib, median of glob > expit( coef(m3)[1]+coef(m3)[2]*(3.2+.1)+coef(m3)[3]*36 )‐ + expit( coef(m3)[1]+coef(m3)[2]*(3.2)+coef(m3)[3]*36 ) (Intercept) 0.037974 > > > #Now let's focus on model m1: > display(m1) glm(formula = y ~ fib, family = binomial(link = "logit"), data = ESR) coef.est coef.se (Intercept) ‐6.85 2.77 fib 1.83 0.90 ‐‐‐ n = 32, k = 2 residual deviance = 24.8, null deviance = 30.9 (difference = 6.0) > > # Nagelkerke's R^2 statistic given by the lrm() function, which is part of > # the Design package. > m1a <‐ lrm(y~fib,data=ESR) > m1a Logistic Regression Model lrm(formula = y ~ fib, data = ESR) Frequencies of Responses 0 1 26 6 Obs Max Deriv Model L.R. d.f. P C Dxy 32 1e‐07 6.04 1 0.0139 0.715 0.429 Gamma Tau‐a R2 Brier 0.432 0.135 0.278 0.109 Coef S.E. Wald Z P Intercept ‐6.845 2.7703 ‐2.47 0.0135 fib 1.827 0.9009 2.03 0.0425 > > > pred2.m1 <‐ predict(m1,type="response") > ESR$predict <‐ pred2.m1 > > ESR$prob.group <‐ cut(pred2.m1,breaks=quantile(pred2.m1, + seq(0,1,.1)), include.lowest=T ) > Hosmer.GOF <‐ sum(unlist(by(ESR, ESR$prob.group, function(x){ + p<‐sum(x$predict) + ((sum(x$y)‐p)^2)/(p*(1‐p/nrow(x))) + }))) > Hosmer.GOF [1] 8.9002 > pchisq(Hosmer.GOF,df=8,lower.tail=F) [1] 0.35079 > > # The scheme that SAS uses to form the groups for the H‐L GOF test is > # somewaht obscure, but I used the following set of cut‐points to create > # the same groups as SAS. Uncommenting this > # code will give the SAS value of the Hosmer‐Lemshow GOF stat > #ESR$prob.group<‐cut(pred2.m1, > # breaks=c(0,.0545,.06,.066,.09,.105,.12,.18,.26,.33,1),include.lowest=T) > #cbind(pred2.m1[order(pred2.m1)],ESR$y[order(pred2.m1)], > # ESR$prob.group[order(pred2.m1)]) > #Hosmer.GOF <‐ sum(unlist(by(ESR, ESR$prob.group, function(x){ > #p<‐sum(x$predict) > #((sum(x$y)‐p)^2)/(p*(1‐p/nrow(x))) > #}))) > #Hosmer.GOF > #pchisq(Hosmer.GOF,df=8,lower.tail=F) > > > > plot(ESR$fib,ESR$y,type="p",xlab="Level of fibrinogen",ylab="y = I( ESR >= 20 )", + main="Fitted probability from model m1") > lines(f0,expit( coef(m1)[1]+coef(m1)[2]*f0 )) > > #rate of change in pi at x such that pi(x)=.5 > # when pi(x)=1/2, x= ‐alpha/beta, and rate of change is beta/4 > slope <‐coef(m1)[2]/4 > abline(a=0.5‐slope*(‐coef(m1)[1]/coef(m1)[2]), b=slope,lty=2) > > #rate of change at x=2.3, the first quartile of fibrinogen > # in general, rate of change at x is beta*pi(x)*(1‐pi(x)) > slope <‐coef(m1)[2]*expit( coef(m1)[1]+coef(m1)[2]*2.3)*(1‐expit( coef(m1)[1]+coef(m1)[2]*2.3)) > abline(a=expit( coef(m1)[1]+coef(m1)[2]*2.3)‐slope*(2.3), b=slope,lty=3) > > ESR$fibc <‐ ESR$fib‐mean(ESR$fib) > m1b <‐ glm(y~fibc,data=ESR,family=binomial(link="logit")) > coef(m1) (Intercept) fib ‐6.8451 1.8271 > coef(m1b) (Intercept) fibc ‐1.7498 1.8271 > > fbar <‐ mean(ESR$fib) > fbar [1] 2.7887 > expit(coef(m1)[1]+coef(m1)[2]*fbar) (Intercept) 0.14807 > expit(coef(m1b)[1]) (Intercept) 0.14807 > > summary(m1b) Call: glm(formula = y ~ fibc, family = binomial(link = "logit"), data = ESR) Deviance Residuals: Min 1Q Median 3Q Max ‐0.930 ‐0.540 ‐0.438 ‐0.336 2.479 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) ‐1.750 0.554 ‐3.16 0.0016 ** fibc 1.827 0.901 2.03 0.0425 * ‐‐‐ Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 30.885 on 31 degrees of freedom Residual deviance: 24.840 on 30 degrees of freedom AIC: 28.84 Number of Fisher Scoring iterations: 5 > linear.hypothesis(m1b, + hypothesis.matrix=matrix( c(0,1),nrow=1,ncol=2,byrow=T), + rhs=c(0),verbose=T) Hypothesis matrix: [,1] [,2] [1,] 0 1 Right‐hand‐side vector: [1] 0 Estimated linear function (hypothesis.matrix %*% coef ‐ rhs) [1] 1.8271 Linear hypothesis test Hypothesis: fibc = 0 Model 1: restricted model Model 2: y ~ fibc Res.Df Df Chisq Pr(>Chisq) 1 31 2 30 1 4.11 0.043 * ‐‐‐ Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Warning message: 'linear.hypothesis' is deprecated. Use 'linearHypothesis' instead. See help("Deprecated") and help("car‐deprecated"). > > > > > pred.m1 <‐ predict(m1,data.frame(fib=f0),se.fit=T,type="link") > L <‐ expit(pred.m1$fit‐1.96*pred.m1$se.fit) > U <‐ expit(pred.m1$fit+1.96*pred.m1$se.fit) > plot(ESR$fib,ESR$y,type="p",xlab="Level of fibrinogen",ylab="y = I( ESR >= 20 )", + main="Fitted probability from model m1") > lines(f0,expit( coef(m1)[1]+coef(m1)[2]*f0 )) > lines(f0,L,lty=4) > lines(f0,U,lty=4) > legend(locator(1),lty=c(1,4),legend=c("Fitted probability","Approx 95% conf. limits")) > > > #Alternatively, the predictions and their standard errors ca n be > # made directly on the probability scale. In this case, standa rd errors > # are estimated with the delta‐method. This approach (com mented out below) > # is not recommended unless the sample size is large > #pred2.m1 <‐ predict(m1,data.frame(fib=f0,y=rep(NA,lengt h=length(f0))),se.fit=T,type="response") > #L2 <‐ pred2.m1$fit‐1.96*pred2.m1$se.fit > #U2 <‐ pred2.m1$fit+1.96*pred2.m1$se.fit > #lines(f0,expit( coef(m1)[1]+coef(m1)[2]*f0 )) > #lines(f0,L2,lty=5) > #lines(f0,U2,lty=5) > > confint.default(m1) 2.5 % 97.5 % (Intercept) ‐12.274734 ‐1.4154 fib 0.061437 3.5927 > confint(m1) Waiting for profiling to be done... 2.5 % 97.5 % (Intercept) ‐13.65654 ‐2.3273 fib 0.33876 3.9985 ...
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