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
install.packages("VGAM", dependencies = T) |
Als Beispiel dient der folgende Datensatz:
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:
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)
##
## Pearson Residuals:
## Min 1Q Median 3Q Max
## logit(P[Y>=2]) -5.3 0.12 0.20 0.41 1.4
## logit(P[Y>=3]) -3.1 -0.77 0.24 0.76 2.6
## logit(P[Y>=4]) -1.4 -0.46 -0.22 -0.11 7.3
##
## Coefficients:
## Estimate Std. Error z value
## (Intercept):1 0.62 0.648 0.96
## (Intercept):2 -1.77 0.654 -2.71
## (Intercept):3 -4.06 0.676 -6.00
## Konflikt -0.58 0.097 -5.94
## Zufriedenh 1.36 0.202 6.75
##
## Number of linear predictors: 3
##
## Names of linear predictors:
## logit(P[Y>=2]), logit(P[Y>=3]), logit(P[Y>=4])
##
## Dispersion Parameter for cumulative family: 1
##
## Residual deviance: 920.5 on 1240 degrees of freedom
##
## Log-likelihood: -460.3 on 1240 degrees of freedom
##
## Number of iterations: 5 |
Mit dem Modell können nun weitere Parameter errechnet werden:
# Koeffizienten
pom.coef =(coef(summary(pom)))
pom.coef |
## Estimate Std. Error z value
## (Intercept):1 0.6210 0.6482 0.958
## (Intercept):2 -1.7703 0.6540 -2.707
## (Intercept):3 -4.0578 0.6759 -6.004
## Konflikt -0.5762 0.0970 -5.940
## Zufriedenh 1.3643 0.2021 6.751 |
# Odds Ratio
pom.odds =exp(coef(pom))
pom.odds |
## (Intercept):1 (Intercept):2 (Intercept):3 Konflikt Zufriedenh
## 1.86078 0.17029 0.01729 0.56205 3.91293 |
# Devianz
pom.devi =deviance(pom)
pom.devi |
# AIC
pom.aic =AIC(pom)
pom.aic |
# logLikelihood
pom.ll =logLik(pom)
pom.ll |
# 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 |
# 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 |
# R^2 Nagelkerke
pom.nagel =as.vector((1 - exp((2/N) * (p0.ll - pom.ll)))/(1 - exp(p0.ll)^(2/N)))
pom.nagel |
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.
# 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
## 2 1236 -456 -4 7.62 0.11 |
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:
# Devianz-Test auf Equal Slopes Assumption
pom.pdevi = (1 - pchisq(deviance(pom) - deviance(npom), df = df.residual(pom) -
df.residual(npom)))
pom.pdevi |
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:
# 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.6615 1.5374 2.7337 0.5706 -1.9684
## Zufriedenh:2 Zufriedenh:3
## -1.2805 -0.8865
##
## Degrees of Freedom: 1245 Total; 1238 Residual
## Residual deviance: 914.5
## Log-likelihood: -457.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.6615 1.5374 2.7337 0.5706 -1.9684
## Zufriedenh:2 Zufriedenh:3
## -1.2805 -0.8865
##
## Degrees of Freedom: 1245 Total; 1238 Residual
## Residual deviance: 914.5
## Log-likelihood: -457.2 |
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
ppom.ce =coef(summary(ppom))
ppom.ce |
## Estimate Std. Error z value
## (Intercept):1 0.6615 0.85866 0.7704
## (Intercept):2 1.5374 0.71998 2.1353
## (Intercept):3 2.7337 0.98084 2.7871
## Konflikt 0.5706 0.09715 5.8730
## Zufriedenh:1 -1.9684 0.34535 -5.6997
## Zufriedenh:2 -1.2805 0.23585 -5.4294
## Zufriedenh:3 -0.8865 0.32543 -2.7242 |
Mit einem kleinen Script können Konfidenzintervalle und andere Werte ausgegeben werden:
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:
get.ci(as.numeric(ppom.ce[4, ])) |
## Estimate lCI uCI SD z p-value OR
## [1,] 0.5706 0.3802 0.761 0.09715 5.873 4.28e-09 1.769 |
Continuation Ratio Model
Um ein Continuation Ratio Modell zu rechnen, muss der Parameter family auf cratio gesetzt werden.
# 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)
##
## Pearson Residuals:
## Min 1Q Median 3Q Max
## logit(P[Y>1|Y>=1]) -4.9 0.18 2.6e-01 3.7e-01 1.2
## logit(P[Y>2|Y>=2]) -3.0 -0.86 3.2e-01 7.2e-01 2.4
## logit(P[Y>3|Y>=3]) -1.5 -0.63 1.1e-05 6.2e-05 11.8
##
## Coefficients:
## Estimate Std. Error z value
## (Intercept):1 0.93 0.574 1.6
## (Intercept):2 -1.14 0.579 -2.0
## (Intercept):3 -3.02 0.607 -5.0
## Konflikt -0.54 0.086 -6.2
## Zufriedenh 1.16 0.179 6.5
##
## Number of linear predictors: 3
##
## Names of linear predictors:
## logit(P[Y>1|Y>=1]), logit(P[Y>2|Y>=2]), logit(P[Y>3|Y>=3])
##
## Dispersion Parameter for cratio family: 1
##
## Residual deviance: 924.4 on 1240 degrees of freedom
##
## Log-likelihood: -462.2 on 1240 degrees of freedom
##
## Number of iterations: 5 |
Dies berechnet standardmäßig die Vorwärts-Methode. Möchte man die Rückwärts-Methode rechnen, setzt man den Parameter reverse auf TRUE.
# 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)
##
## Pearson Residuals:
## Min 1Q Median 3Q Max
## logit(P[Y<2|Y< =2]) -1.6 -0.57 1.5e-05 4.1e-05 7.0
## logit(P[Y<3|Y<=3]) -2.5 -0.88 7.3e-05 7.2e-01 3.1
## logit(P[Y<4|Y<=4]) -6.3 0.18 3.0e-01 4.4e-01 1.1
##
## Coefficients:
## Estimate Std. Error z value
## (Intercept):1 0.13 0.573 0.22
## (Intercept):2 2.09 0.581 3.59
## (Intercept):3 4.04 0.606 6.67
## Konflikt 0.45 0.085 5.33
## Zufriedenh -1.26 0.181 -6.94
##
## Number of linear predictors: 3
##
## Names of linear predictors:
## logit(P[Y<2|Y<=2]), logit(P[Y<3|Y<=3]), logit(P[Y<4|Y<=4])
##
## Dispersion Parameter for cratio family: 1
##
## Residual deviance: 920.5 on 1240 degrees of freedom
##
## Log-likelihood: -460.2 on 1240 degrees of freedom
##
## Number of iterations: 5 |
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:
# 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)
##
## Pearson Residuals:
## Min 1Q Median 3Q Max
## logit(P[Y<2|Y< =2]) -1.8 -0.57 -1.6e-05 0.00015 8.5
## logit(P[Y<3|Y<=3]) -2.5 -0.87 5.3e-07 0.72073 3.1
## logit(P[Y<4|Y<=4]) -5.4 0.18 2.9e-01 0.43071 1.0
##
## Coefficients:
## Estimate Std. Error z value
## (Intercept):1 2.36 1.17 2.01
## (Intercept):2 2.00 0.81 2.46
## (Intercept):3 2.47 1.12 2.20
## Konflikt:1 0.15 0.18 0.81
## Konflikt:2 0.49 0.12 4.07
## Konflikt:3 0.63 0.17 3.74
## Zufriedenh:1 -1.80 0.39 -4.60
## Zufriedenh:2 -1.26 0.26 -4.89
## Zufriedenh:3 -0.84 0.34 -2.47
##
## Number of linear predictors: 3
##
## Names of linear predictors:
## logit(P[Y<2|Y<=2]), logit(P[Y<3|Y<=3]), logit(P[Y<4|Y<=4])
##
## Dispersion Parameter for cratio family: 1
##
## Residual deviance: 914.8 on 1236 degrees of freedom
##
## Log-likelihood: -457.4 on 1236 degrees of freedom
##
## Number of iterations: 5 |
# 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
## 2 1240 -460 4 5.67 0.23 |
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:
# 0-Modell (fuer pseudo R^2) Vorbereitungen
c0 =vglm(Stimmung ~ 1, data = mydata, family = cratio(parallel = T))
c0.ll = logLik(c0)
crm.ll = logLik(crm)
N = length(mydata[, 1]) # Anzahl der Fälle |
# R^2 McFadden
crm.mcfad =as.vector(1 - (crm.ll/c0.ll))
crm.mcfad |
# R^2 Cox&Snell
crm.cox =as.vector(1 - exp((2/N) * (c0.ll - crm.ll)))
crm.cox |
# R^2 Nagelkerke
crm.nagel =as.vector((1 - exp((2/N) * (c0.ll - crm.ll)))/(1 - exp(c0.ll)^(2/N)))
crm.nagel |
# Devianz
crm.devi =deviance(crm)
crm.devi |
# AIC
crm.aic =AIC(crm)
crm.aic |
Mit unserer Funktion von oben können wir uns die Modellparameter der Koeffizienten ausgeben lassen
# Konfidenzinztervalle
crm.ce =coef(summary(crm))
crm.ce |
## Estimate Std. Error z value
## (Intercept):1 0.1282 0.57259 0.2239
## (Intercept):2 2.0875 0.58094 3.5933
## (Intercept):3 4.0425 0.60636 6.6668
## Konflikt 0.4523 0.08488 5.3289
## Zufriedenh -1.2563 0.18108 -6.9378 |
get.ci(as.numeric(crm.ce[4, ])) |
## Estimate lCI uCI SD z p-value OR
## [1,] 0.4523 0.286 0.6187 0.08488 5.329 9.881e-08 1.572 |
get.ci(as.numeric(crm.ce[5, ])) |
## Estimate lCI uCI SD z p-value OR
## [1,] -1.256 -1.611 -0.9014 0.1811 -6.938 3.982e-12 0.2847 |
Links
Literatur
- große Schlarmann, J.; Galatsch, M. (2014): “Regressionsmodelle für ordinale Zielvariablen”. GMS Medizinische Informatik, Biometrie und Epidemiologie, 10(1):2-10
- 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.
