Zur Bestimmung des Einflusses von Prädiktoren auf ein ordinales Ziel stehen verschiedene Modelle der ordinalen Regression zur Verfügung. Namentlich werden hier vor allem Proportional Odds Modelle und Continuation Ratio Modelle verwendet.
Dieser Blogbeitrag beschreibt, wie sich diese Modelle mit Hilfe der freien Statistikumgebung R errechnen lassen.
Vorbereitung
Wir verwenden hier das R-Paket “VGAM”, welches zunächst in R installiert werden muss
1
|
install.packages("VGAM", dependencies=T)
|
Als Beispiel dient der folgende Datensatz:
1
2
3
|
load(url("http://www.produnis.de/R/OrdinalSample.RData"))
mydata <- ordinalSample
head(mydata)
|
## Konflikt Zufriedenh Geschlecht Stimmung
## 1 4.2 2.50 1 maessig
## 6 1.8 3.00 1 sehr gut
## 15 3.4 1.75 1 maessig
## 16 4.0 2.50 2 maessig
## 22 2.2 2.50 1 gut
## 23 3.4 2.25 1 maessig
Die Variable “Stimmung” dient hier als ordinales Ziel mit den Ausprägungen “schlecht” < “maessig” < “gut” < “sehr gut”.
Die Variable “Konflikt” dient als Prädiktor. Sie gibt an, wieviele Konflikte derzeit im Arbeitsleben vorliegen.
Die Variable “Zufriedenh” beschreibt, wie zufrieden Probanden mit ihrem Job sind.
Die Variable “Stimmung” soll nun durch “Konflikt” und “Zufriedenh” beschrieben werden.
Proportional Odds Modell
Ein Proportional Odds Modell errechnet sich leicht mit dem VGAM-Paket:
1
2
3
4
5
|
pom <- vglm(Stimmung ~ Konflikt+Zufriedenh, data=mydata, family=cumulative(parallel=T))
# oder abgekürzt
pom <- vglm(Stimmung ~ Konflikt+Zufriedenh, data=mydata, family=propodds)
summary(pom)
|
##
## Call:
## vglm(formula = Stimmung ~ Konflikt + Zufriedenh, family = propodds,
## data = mydata)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept):1 0.6210 0.6482 0.958 0.3381
## (Intercept):2 -1.7703 0.6540 -2.707 0.0068 **
## (Intercept):3 -4.0578 0.6759 -6.004 1.93e-09 ***
## Konflikt -0.5762 0.0970 -5.940 2.86e-09 ***
## Zufriedenh 1.3643 0.2021 6.751 1.47e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Names of linear predictors: logitlink(P[Y>=2]), logitlink(P[Y>=3]),
## logitlink(P[Y>=4])
##
## Residual deviance: 920.5493 on 1240 degrees of freedom
##
## Log-likelihood: -460.2746 on 1240 degrees of freedom
##
## Number of Fisher scoring iterations: 5
##
## No Hauck-Donner effect found in any of the estimates
##
##
## Exponentiated coefficients:
## Konflikt Zufriedenh
## 0.5620549 3.9129275
Mit dem Modell können nun weitere Parameter errechnet werden:
1
2
|
# Koeffizienten
pom.coef <- (coef(summary(pom)));pom.coef
|
## Estimate Std. Error z value Pr(>|z|)
## (Intercept):1 0.6209962 0.64821816 0.9580049 3.380603e-01
## (Intercept):2 -1.7702692 0.65404383 -2.7066523 6.796539e-03
## (Intercept):3 -4.0577657 0.67586741 -6.0037895 1.927645e-09
## Konflikt -0.5761558 0.09700358 -5.9395317 2.858375e-09
## Zufriedenh 1.3642858 0.20208974 6.7508909 1.469400e-11
1
2
|
# Odds Ratio
pom.odds <- exp(coef(pom));pom.odds
|
## (Intercept):1 (Intercept):2 (Intercept):3 Konflikt Zufriedenh
## 1.8607808 0.1702871 0.0172876 0.5620549 3.9129275
1
2
|
# Devianz und AIC
pom.devi <- deviance(pom);pom.devi
|
## [1] 920.5493
1
|
pom.aic <- AIC(pom);pom.aic
|
## [1] 930.5493
1
2
|
# logLikelihood
pom.ll <- logLik(pom);pom.ll
|
## [1] -460.2746
1
2
3
4
5
|
# 0-Modell (fuer pseudo R^2)
p0 <- vglm(Stimmung ~ 1, data=mydata, family=propodds)
p0.ll <- logLik(p0)
# R^2 McFadden
pom.mcfad <- as.vector(1- (pom.ll/p0.ll));pom.mcfad
|
## [1] 0.1240613
1
2
3
|
# R^2 Cox&Snell
N <- length(mydata[,1]) # Anzahl der Fälle
pom.cox <- as.vector(1 - exp((2/N) * (p0.ll - pom.ll)));pom.cox
|
## [1] 0.2696034
1
2
|
# R^2 Nagelkerke
pom.nagel <- as.vector((1 - exp((2/N) * (p0.ll - pom.ll))) / (1 - exp(p0.ll)^(2/N))) ;pom.nagel
|
## [1] 0.2928789
Das Proportional Odds Modell geht von der “equal slopes assumption” (auch: “proportional odds assumption”) aus. Diese Annahme muss geprüft werden, bevor das Modell als gültig angesehen werden kann. Dies geht mit dem VGAM-Paket recht einfach. Der Befehl vglm(Stimmung ~ Konflikt+Zufriedenh, data=mydata, family=propodds) ist eine Abkürzung für vglm(Stimmung ~ Konflikt+Zufriedenh, data=mydata, family=cumulative(parallel=T)). Der Parameter parallel=TRUE/FALSE stellt ein, ob das Modell mit equal slopes (=TRUE), oder ohne equal slopes assumption (=FALSE) erstellt werden soll. Zur Überprüfung der “Equal Slopes Assumption” erstellt man jeweils ein TRUE und ein FALSE-Modell, und führt dann einen Likelihood-Test durch. Die Nullhypothese lautet, dass equal slopes vorliegen.
1
2
3
4
5
|
# Modell OHNE equal slopes
npom <- vglm(Stimmung ~ Konflikt+Zufriedenh, data=mydata, family=cumulative(parallel=F))
# log-likelihood-Test auf Equal Slopes Assumption
lrtest(pom,npom)# Test
|
## Likelihood ratio test
##
## Model 1: Stimmung ~ Konflikt + Zufriedenh
## Model 2: Stimmung ~ Konflikt + Zufriedenh
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 1240 -460.27
## 2 1236 -456.46 -4 7.6202 0.1065
Der Test ist mit p=0,1065 nicht signifikant. Das heisst, es liegen equal slopes vor. Das Modell darf also verwendet werden.
Alternativ kann dieser Test auch auf die Devianz der Modelle bestimmt werden:
1
2
|
# Devianz-Test auf Equal Slopes Assumption
pom.pdevi = (1 - pchisq(deviance(pom) - deviance(npom), df=df.residual(pom)-df.residual(npom)));pom.pdevi
|
## [1] 0.106526
Ist der Test signifikant, liegen keine “equal slopes” vor. Hier kann ein Partial Proportional Odds Modell erstellt werden, welches die “equal slopes” nur für bestimmte Prädiktoren annimmt. Mit dem VGAM-Paket kann recht einfach bestimmt werden, wie ein solches Modell erstellt werden soll. Hierfür erweitertet man den “parallel”-Switch wie folgt:
1
2
|
# Equal Slopes NUR für "Konflikt"
ppom <- vglm(Stimmung ~ Konflikt+Zufriedenh, data=mydata, family=cumulative(parallel=T~Konflikt-1));ppom
|
##
## Call:
## vglm(formula = Stimmung ~ Konflikt + Zufriedenh, family = cumulative(parallel = T ~
## Konflikt - 1), data = mydata)
##
##
## Coefficients:
## (Intercept):1 (Intercept):2 (Intercept):3 Konflikt Zufriedenh:1
## 0.6614816 1.5373699 2.7337389 0.5705897 -1.9684047
## Zufriedenh:2 Zufriedenh:3
## -1.2805318 -0.8865225
##
## Degrees of Freedom: 1245 Total; 1238 Residual
## Residual deviance: 914.4929
## Log-likelihood: -457.2464
1
2
|
# Nochmals EqualSlopes NUR für "Konflikt"
ppom <- vglm(Stimmung ~ Konflikt+Zufriedenh, data=mydata, family=cumulative(parallel=F~Zufriedenh));ppom
|
##
## Call:
## vglm(formula = Stimmung ~ Konflikt + Zufriedenh, family = cumulative(parallel = F ~
## Zufriedenh), data = mydata)
##
##
## Coefficients:
## (Intercept):1 (Intercept):2 (Intercept):3 Konflikt Zufriedenh:1
## 0.6614816 1.5373699 2.7337389 0.5705897 -1.9684047
## Zufriedenh:2 Zufriedenh:3
## -1.2805318 -0.8865225
##
## Degrees of Freedom: 1245 Total; 1238 Residual
## Residual deviance: 914.4929
## Log-likelihood: -457.2464
parallel=T~Konflikt-1 bedeutet, dass equal slopes nur für Konflikt angenommen wird.parallel=F~Zufriedenh bedeutet, dass equal slopes nur für Zufriedenh nicht angenommen wird.
Beide Befehle bedeuten also das selbe. Daher ist die R-Ausgabe bei beiden Befehlen gleich.
Eine Koeffizientenübersicht erhält man per
1
|
ppom.ce <- coef(summary(ppom));ppom.ce
|
## Estimate Std. Error z value Pr(>|z|)
## (Intercept):1 0.6614816 0.85866412 0.7703613 4.410856e-01
## (Intercept):2 1.5373699 0.71997824 2.1353005 3.273647e-02
## (Intercept):3 2.7337389 0.98084390 2.7871294 5.317723e-03
## Konflikt 0.5705897 0.09715456 5.8730094 4.279541e-09
## Zufriedenh:1 -1.9684047 0.34535219 -5.6997024 1.200167e-08
## Zufriedenh:2 -1.2805318 0.23585229 -5.4293805 5.654999e-08
## Zufriedenh:3 -0.8865225 0.32542569 -2.7241933 6.445877e-03
Mit einem kleinen Script können Konfidenzintervalle und andere Werte ausgegeben werden:
1
2
3
4
5
6
7
8
9
10
|
get.ci <- function(x){
back <-cbind( x[1], # estimate
(x[1] - (1.96 * x[2])), #lCI
(x[1] + (1.96 * x[2])), #uCI
x[2], x[3], # SD und z
(2*(1 -pnorm(abs(x[3]))) ),# p-wert
exp(x[1]))
colnames(back) <- c("Estimate", "lCI", "uCI","SD","z","p-value", "OR")
return(back)
}
|
Der Aufruf erfolgt z.B. so:
1
|
get.ci(as.numeric(ppom.ce[4,]))
|
## Estimate lCI uCI SD z p-value OR
## [1,] 0.5705897 0.3801667 0.7610126 0.09715456 5.873009 4.279541e-09 1.76931
Continuation Ratio Model
Um ein Continuation Ratio Modell zu rechnen, muss der Parameter family auf cratio gesetzt werden.
1
2
3
|
# CR-Modell MIT PO-Assumption (Vorwärts)
crm<- vglm(Stimmung ~ Konflikt+Zufriedenh, data=mydata, family=cratio(parallel=T,reverse=F))
summary(crm)
|
##
## Call:
## vglm(formula = Stimmung ~ Konflikt + Zufriedenh, family = cratio(parallel = T,
## reverse = F), data = mydata)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept):1 0.92610 0.57405 1.613 0.1067
## (Intercept):2 -1.14132 0.57924 -1.970 0.0488 *
## (Intercept):3 -3.01660 0.60706 -4.969 6.72e-07 ***
## Konflikt -0.53645 0.08639 -6.209 5.32e-10 ***
## Zufriedenh 1.16010 0.17944 6.465 1.01e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Names of linear predictors: logitlink(P[Y>1|Y>=1]), logitlink(P[Y>2|Y>=2]),
## logitlink(P[Y>3|Y>=3])
##
## Residual deviance: 924.3824 on 1240 degrees of freedom
##
## Log-likelihood: -462.1912 on 1240 degrees of freedom
##
## Number of Fisher scoring iterations: 6
##
## No Hauck-Donner effect found in any of the estimates
Dies berechnet standardmäßig die Vorwärts-Methode. Möchte man die Rückwärts-Methode rechnen, setzt man den Parameter reverse auf TRUE.
1
2
3
|
# CR-Modell MIT PO-Assumption (Rückwärts)
crm<- vglm(Stimmung ~ Konflikt+Zufriedenh, data=mydata, family=cratio(parallel=T,reverse=TRUE))
summary(crm)
|
##
## Call:
## vglm(formula = Stimmung ~ Konflikt + Zufriedenh, family = cratio(parallel = T,
## reverse = TRUE), data = mydata)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept):1 0.12822 0.57259 0.224 0.822809
## (Intercept):2 2.08751 0.58094 3.593 0.000326 ***
## (Intercept):3 4.04245 0.60636 6.667 2.62e-11 ***
## Konflikt 0.45234 0.08488 5.329 9.88e-08 ***
## Zufriedenh -1.25631 0.18108 -6.938 3.98e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Names of linear predictors: logitlink(P[Y<2|Y<=2]), logitlink(P[Y<3|Y<=3]),
## logitlink(P[Y<4|Y<=4])
##
## Residual deviance: 920.4627 on 1240 degrees of freedom
##
## Log-likelihood: -460.2314 on 1240 degrees of freedom
##
## Number of Fisher scoring iterations: 5
##
## No Hauck-Donner effect found in any of the estimates
Hat man die passende Methode gewählt, folgen die selben Befehle wie beim Proportional Odds Modell.
Zunächst muss überprüft werden, ob “equal slopes” vorliegen:
1
2
|
# CR-Modell OHNE PO-Assumption (Rückwärts)
ncrm<- vglm(Stimmung ~ Konflikt+Zufriedenh, data=mydata, family=cratio(parallel=F,reverse=T));summary(ncrm)
|
##
## Call:
## vglm(formula = Stimmung ~ Konflikt + Zufriedenh, family = cratio(parallel = F,
## reverse = T), data = mydata)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept):1 2.3649 1.1740 2.014 0.043969 *
## (Intercept):2 1.9998 0.8117 2.464 0.013753 *
## (Intercept):3 2.4670 1.1190 2.205 0.027474 *
## Konflikt:1 0.1465 0.1811 0.809 0.418633
## Konflikt:2 0.4864 0.1194 4.073 4.65e-05 ***
## Konflikt:3 0.6308 0.1689 3.736 0.000187 ***
## Zufriedenh:1 -1.8008 0.3919 -4.596 4.32e-06 ***
## Zufriedenh:2 -1.2602 0.2578 -4.889 1.01e-06 ***
## Zufriedenh:3 -0.8351 0.3382 -2.469 0.013543 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Names of linear predictors: logitlink(P[Y<2|Y<=2]), logitlink(P[Y<3|Y<=3]),
## logitlink(P[Y<4|Y<=4])
##
## Residual deviance: 914.7924 on 1236 degrees of freedom
##
## Log-likelihood: -457.3962 on 1236 degrees of freedom
##
## Number of Fisher scoring iterations: 6
##
## No Hauck-Donner effect found in any of the estimates
1
2
|
# loglikelihood test auf equal slopes assumption
lrtest(ncrm,crm)# Test
|
## Likelihood ratio test
##
## Model 1: Stimmung ~ Konflikt + Zufriedenh
## Model 2: Stimmung ~ Konflikt + Zufriedenh
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 1236 -457.40
## 2 1240 -460.23 4 5.6703 0.2252
Der Test ist nicht signifikant (p=0,23). Das bedeutet, dass die Annahme der equal slopes für die hier getesteten Rückwärts-Modelle beibehalten werden kann.
Nun können weitere Modellparameter bestimmt werden:
1
2
3
4
5
6
7
|
# 0-Modell (fuer pseudo R^2)
c0 <- vglm(Stimmung ~ 1, data=mydata, family=cratio(parallel=T))
c0.ll <- logLik(c0)
crm.ll <- logLik(crm)
# R^2 Nagelkerke
crm.nagel <- as.vector((1 - exp((2/N) * (c0.ll - crm.ll))) / (1 - exp(c0.ll)^(2/N))) ;crm.nagel
|
## [1] 0.2930443
1
2
|
# Devianz
crm.devi <- deviance(crm); crm.devi
|
## [1] 920.4627
1
2
|
# AIC
crm.aic <- AIC(crm); crm.aic
|
## [1] 930.4627
Mit unserer Funktion von oben können wir uns die Modellparameter der Koeffizienten ausgeben lassen
1
2
|
# Konfidenzinztervalle
crm.ce <- coef(summary(crm));crm.ce
|
## Estimate Std. Error z value Pr(>|z|)
## (Intercept):1 0.1282227 0.57259420 0.2239329 8.228095e-01
## (Intercept):2 2.0875126 0.58094099 3.5933299 3.264788e-04
## (Intercept):3 4.0424508 0.60636000 6.6667504 2.615293e-11
## Konflikt 0.4523353 0.08488357 5.3288915 9.881398e-08
## Zufriedenh -1.2563060 0.18108066 -6.9378255 3.981811e-12
1
|
get.ci(as.numeric(crm.ce[4,]))
|
## Estimate lCI uCI SD z p-value OR
## [1,] 0.4523353 0.2859635 0.6187071 0.08488357 5.328892 9.881398e-08 1.571979
1
|
get.ci(as.numeric(crm.ce[5,]))
|
## Estimate lCI uCI SD z p-value OR
## [1,] -1.256306 -1.611224 -0.9013879 0.1810807 -6.937825 3.981704e-12 0.2847038
Links und Literatur
- R Data Analysis Examples: Ordinal Logistic Regression
- Thomas Yee (2010): The VGAM Package for Categorical Data Analysis
- Harrell, Frank: Regression Modeling Strategies: With Applications to Linear Models, Logistic Regression, and Survival Analysis. New York: Springer 2001.
- Agresti, A: Analysis of Ordinal Categorical Data. New York: Wiley & Sons, 2nd edition 2010.
- große Schlarmann, J; Galatsch, M (2014): Regressionsmodelle für ordinale Zielvariablen, GMS Medizinische Informatik, Biometrie und Epidemiologie, 10(1):2-10, doi: 10.3205/mibe000154, https://www.egms.de/static/en/journals/mibe/2014-10/mibe000154.shtml
