[R] logistic regression

[Frank E Harrell Jr]

Mary,

The 10-group approach results in a low-resolution and fairly arbitrary 
calibration curve.  Also, it is the basis of the original 
Hosmer-Lemeshow goodness of fit statistic which has been superceded by 
the Hosmer et al single degree of freedom GOF test that does not require 
any binning.  The Design package handles both.  Do ?calibrate.lrm, 
?residuals.lrm, ?lrm for details.

Frank Harrell


Sullivan, Mary M wrote:
> Greetings,
>  
> 
> I am working on a logistic regression model in R and I am struggling with the code, as it is a relatively new program for me.  In searching Google for 'logistic regression diagnostics' I came Elizabeth Brown's Lecture 14 from her Winter 2004 Biostatistics 515 course  (http://courses.washington.edu/b515/l14.pdf) .  I found most of the code to be very helpful, but I am struggling with the lines on to calculate the observed and expected values in the 10 groups created by the cut function.  I get error messages in trying to create the E and O matrices:  R won't accept assignment of "fi1c==j" and it won't calculate the sum.  
> 
>  
> 
> I am wondering whether someone might be able to offer me some assistance...my search of the archives was not fruitful.
> 
>  
> 
> Here is the code that I adapted from the lecture notes:
> 
>  
> 
> fit <- fitted(glm.lyme)
> 
> fitc <- cut(fit, br = c(0, quantile(fit, p = seq(.1, .9, .1)),1))  
> 
> t<-table(fitc)
> 
> fitc <- cut(fit, br = c(0, quantile(fit, p = seq(.1, .9, .1)), 1), labels = F)
> 
> t<-table(fitc)
> 
>  
> 
> #Calculate observed and expected values in ea group
> 
> E <- matrix(0, nrow=10, ncol = 2)
> 
> O <- matrix(0, nrow=10, ncol=2)
> 
> for (j in 1:10) {
> 
>       E[j, 2] = sum(fit[fitc==j])
> 
>       E[j, 1] = sum((1- fit)[fitc==j])
> 
>       O[j, 2] = sum(pcdata$lymdis[fitc==j])
> 
>       O[j, 1] = sum((1-pcdata$lymdis)[fitc==j])
> 
>             
> 
> }
> 
>  
> 
> Here is the error message:  Error in Summary.factor(..., na.rm = na.rm) : 
>         sum not meaningful for factors
> 
>  
> 
>  
> 
> I understand what it means; I just can't figure out how to get around it or how to get the output printed in table form.  Thank you in advance for any assistance.
> 
>  
> 
> Mary Sullivan
> 
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
> 


-- 
Frank E Harrell Jr   Professor and Chair           School of Medicine
                      Department of Biostatistics   Vanderbilt University