Logistic Regression Fit
require(mosaic)
require(mosaicData)
require(lmtest)
cols = trellis.par.get()$superpose.symbol$colTwo-Way Tables and Chi-Squared Tests
Recall the logistic model for \(Acceptance\) to medical school that we built previously, based on \(Sex\). This is a simple logistic regression model in which both the response and explanatory variable are binary.
MedGPA = read.csv("http://www.math.smith.edu/~bbaumer/mth247/MedGPA.csv")
logm = glm(Acceptance ~ Sex, data=MedGPA, family=binomial)
summary(logm)##
## Call:
## glm(formula = Acceptance ~ Sex, family = binomial, data = MedGPA)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.435 -1.084 0.940 0.940 1.274
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.5878 0.3944 1.490 0.136
## SexM -0.8109 0.5528 -1.467 0.142
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 75.791 on 54 degrees of freedom
## Residual deviance: 73.594 on 53 degrees of freedom
## AIC: 77.594
##
## Number of Fisher Scoring iterations: 4As we can see, the \(p\)-value for the coefficient of \(Sex\) is not significant at the 10% level.
Since both variables in this case are binary, an essentially equivalent model is to simply compute a 2-way table.
two.way = tally(~Acceptance | Sex, data=MedGPA, format="count")
two.way## Sex
## Acceptance F M
## 0 10 15
## 1 18 12Note that the logistic model only has two possible imputs (and therefore two possible outputs), and that they match the entries in the 2-way table.
fit.accept = makeFun(logm)
fit.accept(Sex="M")## 1
## 0.4444444fit.accept(Sex="F")## 1
## 0.6428571Furthermore, the odds ratio from the 2-way table is the same as the odds ratio from the logistic model.
oddsRatio(two.way)## [1] 2.25# Since the coefficients is negative, we add a negative here to match the 2-way table
exp(-coef(logm))## (Intercept) SexM
## 0.5555556 2.2500000When we built the logistic regression model, we obtained a \(p\)-value for the coefficient of \(Sex\). This was helpful, because it gave us some information about the extent of the statistical significance of the relationship between \(Sex\) and \(Acceptance\). Can we obtain something similar from the 2-way table? Indeed, the \(\chi^2\)-test gives essentially the same information.
chisq.test(two.way[1:2,1:2], correct=FALSE)##
## Pearson's Chi-squared test
##
## data: two.way[1:2, 1:2]
## X-squared = 2.1825, df = 1, p-value = 0.1396Classification
We already learned how to compare models using the deviance, but how do we know how well our model works? One techinque for assessing the goodness-of-fit in a logistic regression model is to examine the percentage of the time that our model was “right.”
Recall the logistic model that we built for \(isAlive\) as a function of \(age\) and \(smoker\).
Whickham = Whickham %>%
mutate(isAlive = 2 - as.numeric(outcome))
logm = glm(isAlive ~ age + smoker, data=Whickham, family=binomial(logit))
summary(logm)##
## Call:
## glm(formula = isAlive ~ age + smoker, family = binomial(logit),
## data = Whickham)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.2795 -0.4381 0.2228 0.5458 1.9581
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 7.599221 0.441231 17.223 <2e-16 ***
## age -0.123683 0.007177 -17.233 <2e-16 ***
## smokerYes -0.204699 0.168422 -1.215 0.224
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1560.32 on 1313 degrees of freedom
## Residual deviance: 945.02 on 1311 degrees of freedom
## AIC: 951.02
##
## Number of Fisher Scoring iterations: 6When we visualize the model, it’s hard to assess how often our model is correct.
plotPoints(jitter(isAlive) ~ age, groups=smoker, data=Whickham, alpha=0.3, pch=19, cex=2,
ylab="Probability of Being Alive (units)",
xlab="Age (years)", main="Whickham Study Outcomes",
sub=paste("Number of Cases = ", nrow(Whickham)),
auto.key=TRUE)
fit.outcome = makeFun(logm)
plotFun(fit.outcome(age=x, smoker="Yes") ~ x, add=TRUE)
plotFun(fit.outcome(age=x, smoker="No") ~ x, col=cols[2], add=TRUE)
Moreover, what does it even mean to be correct? The response variable is binary, but the predictions generated by the model are quantities on \([0,1]\). A simple way to classify the fitted values of our model is to simply round them off. Once we do this, we can tabulate how often the rounded-off probability from the model agrees with the actual response variable.
Whickham = Whickham %>%
mutate(fit.alive = ifelse(logm$fitted.values >= 0.5, 1, 0))
tally(~isAlive | fit.alive, data=Whickham)## fit.alive
## isAlive 0 1
## 0 251 118
## 1 94 851tbl = tally(~isAlive | fit.alive, data=Whickham, format="count")
sum(diag(tbl)) / nrow(Whickham)## [1] 0.8386606Thus, our model was correct nearly 84% of the time. Is that good? It depends on the application…
We might want to compare this with what would happen if we used the null model, or if we had just flipped a coin instead of using a model at all.
mean(~isAlive, data=Whickham)## [1] 0.7191781Whickham = Whickham %>%
mutate(fit.alive = sample(c(0,1), size=1314, replace=TRUE))
tally(~isAlive | fit.alive, data=Whickham)## fit.alive
## isAlive 0 1
## 0 166 203
## 1 482 463Concordance
A more sophisticated technique for assessing the quality of fit is to examine the \(C\)-statistic, where \(C\) stands for concordance. The idea here is to pair up each actual success with each actual failure, and then compute whether the fitted probability for the success was greater than the fitted probability for the failure. The percentage of pairs statisfying this condition makes up the \(C\)-statistic.
Note that even if our model had no predictive power, it would still be right about half of the time. Thus, the \(C\)-statistic for two randomly generated vectors is about 0.5.
X = data.frame(a = runif(10000), b = runif(10000))
require(Hmisc)## Loading required package: Hmisc## Loading required package: survival## Loading required package: Formula##
## Attaching package: 'Hmisc'## The following objects are masked from 'package:dplyr':
##
## src, summarize## The following objects are masked from 'package:base':
##
## format.pval, unitsrcorrcens(a~b, data=X)##
## Somers' Rank Correlation for Censored Data Response variable:a
##
## C Dxy aDxy SD Z P n
## b 0.502 0.004 0.004 0.007 0.54 0.587 10000In this case our model is much more effective.
rcorrcens(isAlive ~ fitted.values(logm), data=Whickham)##
## Somers' Rank Correlation for Censored Data Response variable:isAlive
##
## C Dxy aDxy SD Z P n
## fitted.values(logm) 0.893 0.786 0.786 0.02 39.6 0 1314