# # Epage 266 # # Written by: # -- # John L. Weatherwax 2009-04-21 # # email: wax@alum.mit.edu # # Please send comments and especially bug reports to the # above email address. # #----- if(is.null(version$language) == FALSE){ require(alr3) }else{ library(alr3) } if(is.null(version$language) == FALSE) data(downer) # This should be applied to ONLY the features we will be using in the regressions below. # We don't want to delete a sample if some unrelated/unused feature is NA # downer = na.omit( downer ) attach(downer) postscript("../../WriteUp/Graphics/Chapter12/prob_1_scatterplot.eps", onefile=FALSE, horizontal=FALSE) plot( jitter(Myopathy,amount=0.05), jitter(Outcome,amount=0.05), xlab="Myopathy", ylab="Outcome" ) dev.off() # calculate the percentage of surviving cows given each setting of Myopathy: # inds = Myopathy==0 fSurvivingM0 = sum(Outcome[inds])/length(inds) # one if survived inds = Myopathy==1 fSurvivingM1 = sum(Outcome[inds])/length(inds) # fit a logistic regression model to this data: # m0 = glm( Outcome ~ Myopathy, family=binomial(link="logit"), data=downer ) s0 = summary(m0) s0 # compare this to the values obtained by the bootstrap in Chapter 4: # alpha = 0.95 nConf <- qnorm( 1 - (1-alpha)/2 ) estimate = s0$coefficients[2,1] se = s0$coefficients[2,2] print( estimate - nConf * se ) print( estimate + nConf * se ) # obtain the estimated value of $\theta(x)$ probability of survival for the two cases of Myopathy # fSurvLRM0 = predict( m0, data.frame(Myopathy=0), type="response" ) fSurvLRM1 = predict( m0, data.frame(Myopathy=1), type="response" ) print( fSurvLRM0 ) print( fSurvLRM1 ) plot( CK, jitter(Outcome,amount=0.05), xlab="CK", ylab="Outcome" ) plot( log(CK), jitter(Outcome,amount=0.05), xlab="log(CK)", ylab="Outcome" ) # fit a logistic regression model to this data: # downer$logCK = log(CK) m1 = glm( Outcome ~ logCK, family=binomial(link="logit"), data=downer ) s1 = summary(m1) s1 # WWX finish ... # m2 = glm( Outcome ~ logCK + Myopathy + Myopathy*logCK, family=binomial(link="logit"), data=downer ) s2 = summary(m2) s2