UCLA (University of California, Los Angeles)
A researcher is interested in how variables, such as GRE (Graduate Record Exam scores), GPA (grade point average) and prestige of the undergraduate institution, effect admission into graduate school. The response variable, admit/don’t admit, is a binary variable.
Variable | Bedeutung |
---|---|
admit | Zulassung zu graduate school (binär) |
gre | graduate record exmainations (Zulassungsprüfungspunkte) |
gpa | grade point average (in etwa Notendurchschnitt: hier nach descriptives mit Maximum 4 (Bestnote) (A=4, B=3, C=2, D=1, F=0)) |
rank | Bewertung der Qualtiät Institution an der man vorher gelernt hat (vier Stufen, 1..4)) |
Fragen:
[res_begin:res_1]
# original data location
# df <- read.csv("http://www.ats.ucla.edu/stat/dataa/binary.csv")
# alternatively
df <- read.csv("http://md.psych.bio.uni-goettingen.de/mv/data/div/example_ucla.csv")
# rank might not be considered intervall-scaled
df$rank <- factor(df$rank)
# we fit a model
m.1 <- glm(admit ~ gre + gpa + rank, data = df, family = "binomial")
# and take a look at the result summary
summary(m.1)
##
## Call:
## glm(formula = admit ~ gre + gpa + rank, family = "binomial",
## data = df)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.6268 -0.8662 -0.6388 1.1490 2.0790
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.989979 1.139951 -3.500 0.000465 ***
## gre 0.002264 0.001094 2.070 0.038465 *
## gpa 0.804038 0.331819 2.423 0.015388 *
## rank2 -0.675443 0.316490 -2.134 0.032829 *
## rank3 -1.340204 0.345306 -3.881 0.000104 ***
## rank4 -1.551464 0.417832 -3.713 0.000205 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 499.98 on 399 degrees of freedom
## Residual deviance: 458.52 on 394 degrees of freedom
## AIC: 470.52
##
## Number of Fisher Scoring iterations: 4
Im Summary:
df <- read.csv("http://md.psych.bio.uni-goettingen.de/mv/data/div/example_ucla.csv")
df$rank <- factor(df$rank)
m.1 <- glm(admit ~ gre + gpa + rank, data = df, family = "binomial")
# overall test:
m.0 <- glm(admit ~ 1, data = df, family = "binomial")
anova(m.0, m.1, test="Chisq")
## Analysis of Deviance Table
##
## Model 1: admit ~ 1
## Model 2: admit ~ gre + gpa + rank
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 399 499.98
## 2 394 458.52 5 41.459 7.578e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# test predictive variables (null model is automatically generated, deviance of null model can be found in model$null.deviance)
m.1$null.deviance
## [1] 499.9765
m.1$deviance
## [1] 458.5175
anova(m.1, test="Chisq")
## Analysis of Deviance Table
##
## Model: binomial, link: logit
##
## Response: admit
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 399 499.98
## gre 1 13.9204 398 486.06 0.0001907 ***
## gpa 1 5.7122 397 480.34 0.0168478 *
## rank 3 21.8265 394 458.52 7.088e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
In den Modell-Tests - das Gesamtmodell unterscheidet sich signifikant vom Null-Modell - jede der drei Vorhersagevariablen hat eine signifikante Vorhersagekraft für die Reaktionsvariable
Um die Konfidenzintervalle der Koeffizienten zu erhalten, können wir confint()
benutzen
df <- read.csv("http://md.psych.bio.uni-goettingen.de/mv/data/div/example_ucla.csv")
df$rank <- factor(df$rank)
m.1 <- glm(admit ~ gre + gpa + rank, data = df, family = "binomial")
## CIs using profiled log-likelihood
confint(m.1)
## Waiting for profiling to be done...
## 2.5 % 97.5 %
## (Intercept) -6.2716202334 -1.792547080
## gre 0.0001375921 0.004435874
## gpa 0.1602959439 1.464142727
## rank2 -1.3008888002 -0.056745722
## rank3 -2.0276713127 -0.670372346
## rank4 -2.4000265384 -0.753542605
## CIs using standard errors
confint.default(m.1)
## 2.5 % 97.5 %
## (Intercept) -6.2242418514 -1.755716295
## gre 0.0001202298 0.004408622
## gpa 0.1536836760 1.454391423
## rank2 -1.2957512650 -0.055134591
## rank3 -2.0169920597 -0.663415773
## rank4 -2.3703986294 -0.732528724
Wir können den Gesamteffekt von “rank” mit der Funktion aod::wald.test()
prüfen. Die Reihenfolge, mit der die Koeffizienten in der Tabelle der Koeffizienten ausgegeben wird, ist dieselbe, wie die Reihenfolge der Terme im Modell. Das ist wichtig, weil die Funktion wald.test()
Bezug nimmt auf die Reihenfolge der Koeffizienten im Modell. Wir benutzen wald.test()
b
zeigt die Koeffizienten, unter Sigma
finden wir die Varianz-Kovarianz-Matrix der Error-Terme und schließlich sagt uns R, welche Terme im Modell getestet werden. In unserem Fall sind das die Terme 4, 5 und 6 für die Levels der Variable “rank” vom Typ factor
.
Wir können auch Einzelvergleiche von bestimmten Ausprägungen der Variable “rank” anfordern. Im Beispiel unten testen wir, ob die Koeffizienten für “rank”=2 und “rank”=3 verschieden sind. Hierzu erzeugen wir einen Vektor “l”, der einen Kontrast definiert, den wir prüfen wollen. In vorliegenden Fall prüfen wird die Differenz zwischen den Termen für “rank”=2 und “rank”=3, also zwischen dem 4. und 5. Term im Modell. Die interessierenden Terme bekommen 1 bzw. -1, der Rest jeweils eine 0 als Kontrastkoeffizient. Danach wird “l” im wald.test()
als zusätzlicher Parameter übergeben und wir sagen R, dass “l” an Stelle von “Terms” geprüft werden soll, wie vorher.
df <- read.csv("http://md.psych.bio.uni-goettingen.de/mv/data/div/example_ucla.csv")
df$rank <- factor(df$rank)
m.1 <- glm(admit ~ gre + gpa + rank, data = df, family = "binomial")
# aod::wald.test() gives us an overall effect of "rank"
require(aod)
## Loading required package: aod
wald.test(b = coef(m.1), Sigma = vcov(m.1), Terms = 4:6)
## Wald test:
## ----------
##
## Chi-squared test:
## X2 = 20.9, df = 3, P(> X2) = 0.00011
# we can also test specific differences of levels of "rank"
l <- cbind(0,0,0,1,-1,0)
wald.test(b = coef(m.1), Sigma = vcov(m.1), L = l)
## Wald test:
## ----------
##
## Chi-squared test:
## X2 = 5.5, df = 1, P(> X2) = 0.019
Für die \(\chi^2\)-Statistik des overall-Tests von 20.9 mit 3 Freiheitsgraden erhalten wir einen p-Wert von 0.00011, was uns sagt, dass der Gesamteffekt von “rank” statistisch signifikant ist.
Einem \(\chi^2\)-Test-Statistik von 5.5 mit einem Freiheitsgrad ist ein p-Wert von 0.019 zugeordnet. Wir interpretieren also, dass sich Level 2 und 3 der Variable “rank” statistisch signifikant voneinander unterscheiden.
Wenn wir die Koeffizienten exponieren, kommen wir direkt an die Odds der Vorhersagevariable. Hierzu dient exp()
. Die Koeffizienten können wir aus dem Ergebnisobjekt abrufen oder coef()
benutzen und das Ergebnisobjekt des Modells übergeben. Das können wir auch für die Konfidenzintervalle tun. Darüber hinaus können wir auch eine kleine Ergebnistabelle generieren.
df <- read.csv("http://md.psych.bio.uni-goettingen.de/mv/data/div/example_ucla.csv")
df$rank <- factor(df$rank)
m.1 <- glm(admit ~ gre + gpa + rank, data = df, family = "binomial")
## odds ratios only
exp(coef(m.1))
## (Intercept) gre gpa rank2 rank3 rank4
## 0.0185001 1.0022670 2.2345448 0.5089310 0.2617923 0.2119375
## odds ratios and 95% CI
exp(cbind(OR = coef(m.1), confint(m.1)))
## Waiting for profiling to be done...
## OR 2.5 % 97.5 %
## (Intercept) 0.0185001 0.001889165 0.1665354
## gre 1.0022670 1.000137602 1.0044457
## gpa 2.2345448 1.173858216 4.3238349
## rank2 0.5089310 0.272289674 0.9448343
## rank3 0.2617923 0.131641717 0.5115181
## rank4 0.2119375 0.090715546 0.4706961
Nun können wir sagen, dass für eine Anstieg von “gpa” die odds der Zulassung um einen Faktor von 2.23 steigen. Wegen der Nicht-Linearität, logarithmierte Odds werden vorhergesagt, ist bei der Interpretation allerdings Vorsicht geboten. Auch wenn er mit ausgegeben wird, interpretieren wir den Koeffizienten für den Intercept normalerweise nicht.
Wir können auch die Probabilities direkt vorhersagen, um das Modell besser zu verstehen. Vorhergesagte Wahrscheinlichkeiten können für stetige, aber auch für kategoriale Vorhersagewerte berechnet werden. Hierzu müssen wir zunächst einen Data-Frame erstellen, in dem die Werte enthalten sind, die wir für die Vorhersage brauchen.
Wir fangen damit an, die vorhergesagten Wahrscheinlichkeiten der Zulassung für jeden Level der Variable “rank” an ihrem Mittelwert zu berechnen. Wir kreieren ein neues Dataframe.
These objects must have the same names as the variables in your logistic regression above (e.g. in this example the mean for gre must be named gre). Now that we have the data frame we want to use to calculate the predicted probabilities, we can tell R to create the predicted probabilities. The first line of code below is quite compact, we will break it apart to discuss what various components do. The newdata1$rankP tells R that we want to create a new variable in the dataset (data frame) newdata1 called rankP, the rest of the command tells R that the values of rankP should be predictions made using the predict( ) function. The options within the parentheses tell R that the predictions should be based on the analysis mylogit with values of the predictor variables coming from newdata1 and that the type of prediction is a predicted probability (type=“response”). The second line of the code lists the values in the data frame newdata1. Although not particularly pretty, this is a table of predicted probabilities.
In the above output we see that the predicted probability of being accepted into a graduate program is 0.52 for students from the highest prestige undergraduate institutions (rank=1), and 0.18 for students from the lowest ranked institutions (rank=4), holding gre and gpa at their means. We can do something very similar to create a table of predicted probabilities varying the value of gre and rank. We are going to plot these, so we will create 100 values of gre between 200 and 800, at each value of rank (i.e., 1, 2, 3, and 4).
The code to generate the predicted probabilities (the first line below) is the same as before, except we are also going to ask for standard errors so we can plot a confidence interval. We get the estimates on the link scale and back transform both the predicted values and confidence limits into probabilities.
It can also be helpful to use graphs of predicted probabilities to understand and/or present the model. We will use the ggplot2 package for graphing. Below we make a plot with the predicted probabilities, and 95% confidence intervals.
df <- read.csv("http://md.psych.bio.uni-goettingen.de/mv/data/div/example_ucla.csv")
df$rank <- factor(df$rank)
m.1 <- glm(admit ~ gre + gpa + rank, data = df, family = "binomial")
df.n <- with(df, data.frame(gre = mean(gre), gpa = mean(gpa), rank = factor(1:4)))
## view data frame
df.n
## gre gpa rank
## 1 587.7 3.3899 1
## 2 587.7 3.3899 2
## 3 587.7 3.3899 3
## 4 587.7 3.3899 4
# df.n$rankP <- predict(mylogit, newdata = newdata1, type = "response")
df.n$rankP <- predict(m.1, newdata = df.n, type = "response")
df.n
## gre gpa rank rankP
## 1 587.7 3.3899 1 0.5166016
## 2 587.7 3.3899 2 0.3522846
## 3 587.7 3.3899 3 0.2186120
## 4 587.7 3.3899 4 0.1846684
# see comments below
with(m.1, null.deviance - deviance)
## [1] 41.45903
with(m.1, df.null - df.residual)
## [1] 5
with(m.1, pchisq(null.deviance - deviance, df.null - df.residual, lower.tail = FALSE))
## [1] 7.578194e-08
logLik(m.1)
## 'log Lik.' -229.2587 (df=6)
We may also wish to see measures of how well our model fits. This can be particularly useful when comparing competing models. The output produced by summary(mylogit) included indices of fit (shown below the coefficients), including the null and deviance residuals and the AIC. One measure of model fit is the significance of the overall model. This test asks whether the model with predictors fits significantly better than a model with just an intercept (i.e., a null model). The test statistic is the difference between the residual deviance for the model with predictors and the null model. The test statistic is distributed chi-squared with degrees of freedom equal to the differences in degrees of freedom between the current and the null model (i.e., the number of predictor variables in the model). To find the difference in deviance for the two models (i.e., the test statistic) we can use the command:
41.4590251
The degrees of freedom for the difference between the two models is equal to the number of predictor variables in the mode, and can be obtained using:
5
Finally, the p-value can be obtained using:
7.5781942^{-8}
The chi-square of 41.46 with 5 degrees of freedom and an associated p-value of less than 0.001 tells us that our model as a whole fits significantly better than an empty model. This is sometimes called a likelihood ratio test (the deviance residual is -2*log likelihood). To see the model’s log likelihood, we type:
-229.2587462
# and a plot
newdata2 <- with(df,
data.frame(gre = rep(seq(from = 200, to = 800, length.out = 100), 4),
gpa = mean(gpa), rank = factor(rep(1:4, each = 100))))
# todo
source zu Logistischer Regression auf der UCLA Website.
[res_end]