Greeting R Community, I'm trying to learn Logistic Regression on my own and am using An Introduction to Logistic Regression Analysis and Reporting (Peng, C., Lee, K., & Ingersoll, G. ,2002). This article uses a Score Test Stat as a measure of overall fit for a logistic regression model. The author calculates this statistic using SAS. I am looking for an [R] function that can compute this stat and a p=value (please don't critique the stat itself). As far as I can tell glm.scoretest() from library(statmod) is the only function for this but it does not give a pvalue and appears to not correspond with the author's example. Some chat on the discussion board a while back indicated that this was a well known test Click here for that thread but difficult to calculate in [R]. I'm wondering if there is a way to do this fairly easily in [R] at this time? I am working with the data set from the study and am looking to get close to the 9.5177 score test stat and .0086 p-value with 2 degrees freedom. Below is the code for the data set and some descriptives so you can see the data set is the same as from the study. I highlighted my attempt to get a score test statistic and am curious if this is it (minus the p-value). #BEGINNING OF CODE id<-factor(1:189) gender<-factor(c("Boy","Boy","Girl","Girl","Girl","Boy","Girl","Boy","Girl","Girl","Boy","Boy","Boy","Boy","Girl","Girl","Boy","Girl","Boy","Girl","Boy","Girl","Girl", "Boy","Boy","Girl","Girl","Girl","Boy","Boy","Boy","Girl","Boy","Girl","Boy","Girl","Girl","Girl","Girl","Girl","Boy","Girl","Boy","Girl","Girl","Girl", "Girl","Boy","Girl","Boy","Girl","Boy","Girl","Girl","Boy","Boy","Boy","Boy","Boy","Boy","Boy","Boy","Boy","Girl","Boy","Boy","Boy","Boy","Girl","Boy", "Girl","Boy","Boy","Boy","Girl","Boy","Girl","Girl","Boy","Girl","Girl","Girl","Boy","Boy","Boy","Boy","Boy","Girl","Girl","Girl","Girl","Boy","Girl", "Girl","Girl","Girl","Girl","Girl","Girl","Girl","Girl","Girl","Boy","Girl","Boy","Boy","Girl","Girl","Girl","Boy","Girl","Boy","Girl","Girl","Girl","Boy", "Girl","Boy","Girl","Boy","Girl","Boy","Girl","Girl","Girl","Girl","Girl","Girl","Girl","Girl","Boy","Girl","Boy","Boy","Boy","Boy","Boy","Boy","Boy","Girl", "Girl","Girl","Boy","Boy","Girl","Girl","Boy","Girl","Boy","Boy","Boy","Girl","Girl","Girl","Girl","Boy","Boy","Girl","Boy","Boy","Girl","Boy","Boy","Boy", "Boy","Girl","Boy","Boy","Girl","Girl","Boy","Boy","Boy","Boy","Boy","Girl","Girl","Girl","Girl","Boy","Boy","Boy","Girl","Boy","Girl","Boy","Boy","Boy","Girl")) reading.score<-c(91.0,77.5,52.5,54.0,53.5,62.0,59.0,51.5,61.5,56.5,47.5,75.0,47.5,53.5,50.0,50.0,49.0,59.0,60.0,60.0, 60.5,50.0,101.0,60.0,60.0,83.5,61.0,75.0,84.0,56.5,56.5,45.0,60.5,77.5,62.5,70.0,69.0,62.0,107.5,54.5,92.5,94.5,65.0, 80.0,45.0,45.0,66.0,66.0,57.5,42.5,60.0,64.0,65.0,47.5,57.5,55.0,55.0,76.5,51.5,59.5,59.5,59.5,55.0,70.0,66.5,84.5, 57.5,125.0,70.5,79.0,56.0,75.0,57.5,56.0,67.5,114.5,70.0,67.0,60.5,95.0,65.5,85.0,55.0,63.5,61.5,60.0,52.5,65.0,87.5, 62.5,66.5,67.0,117.5,47.5,67.5,67.5,77.0,73.5,73.5,68.5,55.0,92.0,55.0,55.0,60.0,120.5,56.0,84.5,60.0,85.0,93.0,60.0, 65.0,58.5,85.0,67.0,67.5,65.0,60.0,47.5,79.0,80.0,57.5,64.5,65.0,60.0,85.0,60.0,58.0,61.5,60.0,65.0,93.5,52.5,42.5, 75.0,48.5,64.0,66.0,82.5,52.5,45.5,57.5,65.0,46.0,75.0,100.0,77.5,51.5,62.5,44.5,51.0,56.0,58.5,69.0,65.0,60.0,65.0, 65.0,40.0,55.0,52.5,54.5,74.0,55.0,60.5,50.0,48.0,51.0,55.0,93.5,61.0,52.5,57.5,60.0,71.0,65.0,60.0,55.0,60.0,77.0, 52.5,95.0,50.0,47.5,50.0,47.0,71.0,65.0) reading.recommendation<-as.factor(c(rep("no",130),rep("yes",59))) DF<-data.frame(id,gender,reading.score,reading.recommendation) head(DF) #================================================================================= # DESCRIPTIVES #================================================================================= table(DF[2:4]) #BREAKDOWN OF SCORES BY GENDER AND REMEDIAL READING RECOMENDATIONS table(DF[c(2)]) #TABLE OF GENDER print(table(DF[c(2)])/sum(table(DF[c(4)]))*100,digits=4)#PERCENT GENDER BREAKDOWN table(DF[c(4)]) #TABLE RECOMENDDED FOR REMEDIAL READING print(table(DF[c(4)])/sum(table(DF[c(4)]))*100,digits=4)#Probability of Reccomended table(DF[c(2,4)]) #TABLE OF GIRLS AND BOYS RECOMENDDED FOR REMEDIAL READING print(prop.table(table(DF[c(2,4)]),1)*100,digits=4)#Probability of Reccomended #================================================================================= # ANALYSIS #================================================================================= (mod1<-glm(reading.recommendation~reading.score+gender,family=binomial,data=DF)) library(statmod) with(DF,glm.scoretest(mod1, c(0,2,3), dispersion=NULL))^2 #If I move the decimal over 2 to the right I get close to the 9.5177 from the study (with(DF,glm.scoretest(mod1, c(0,2,3), dispersion=NULL))^2)*100 #is this it? #END OF CODE I am running R 2.13.0 in a windows 7 machine. Peng, C., Lee, K., & Ingersoll, G. (2002). An Introduction to Logistic Regression Analysis and Reporting. The Journal of Educational Research, 96(1), 3-14. Heldref Publications. Retrieved from http://www.informaworld.com/openurl?genre=article&doi=10.1080/00220670209598786&magic=crossref [[alternative HTML version deleted]]
______________________________________________ R-help@r-project.org 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.