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 die R-Pakete „rms“ und „VGAM„, welche zunächst in R installiert werden müssen
install.packages(c("car", "rms","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) |
Die Variable „Mood“ dient hier als ordinales Ziel mit den Ausprägungen „schlecht“ < „maessig“ < „gut“ < „sehr gut“.
Die Variable „JobFamily“ dient als Prädiktor. Sie gibt an, inwieweit Familie und Job konfligieren, wobei hohe Werte für viele Konflikte stehen.
Die Variable „JobSatisfaction“ beschreibt, wie zufrieden Probanden mit ihrem Job sind.
Die Variable „Mood“ soll nun durch „JobFamily“ und „JobSatisfaction“ beschrieben werden.
Proportional Odds Modell
Ein Proportional Odds Modell errechnet sich leicht mit dem VGAM-Paket:
require(VGAM) pom < - vglm(Mood ~ JobFamily+JobSatisfaction, data=mydata, family=propodds) summary(pom)
Call:
vglm(formula = Mood ~ JobFamily + JobSatisfaction, family = propodds,
data = mydata)
Pearson Residuals:
Min 1Q Median 3Q Max
logit(P[Y>=2]) -5.2907 0.11698 0.19564 0.41275 1.4296
logit(P[Y>=3]) -3.1412 -0.76542 0.23575 0.76141 2.5996
logit(P[Y>=4]) -1.3935 -0.46309 -0.21527 -0.11447 7.2513
Coefficients:
Estimate Std. Error z value
(Intercept):1 0.62100 0.648218 0.9580
(Intercept):2 -1.77027 0.654044 -2.7067
(Intercept):3 -4.05777 0.675867 -6.0038
JobFamily -0.57616 0.097004 -5.9395
JobSatisfaction 1.36429 0.202090 6.7509
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.5493 on 1240 degrees of freedom
Log-likelihood: -460.2747 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
# Odds Ratio pom.odds < - exp(coef(pom));pom.odds
# Devianz und AIC pom.devi < - deviance(pom);pom.devi pom.aic <- AIC(pom);pom.aic
# logLikelihood pom.ll < - logLik(pom);pom.ll
# 0-Modell (für pseudo R^2) p0 < - vglm(WAIo ~ 1, data=mydata, family=propodds);p0 p0.ll <- logLik(p0);p0.ll # 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(Mood ~ JobFamily+JobSatisfaction, data=mydata, family=propodds) ist eine Abkürzung für vglm(Mood ~ JobFamily+JobSatisfaction, 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(Mood ~ JobFamily+JobSatisfaction, data=mydata, family=cumulative(parallel=F)) # log-likelihood-Test auf Equal Slopes Assumption lrtest(pom,npom)# Test
Likelihood ratio test Model 1: Mood ~ JobFamily + JobSatisfaction Model 2: Mood ~ JobFamily + JobSatisfaction #Df LogLik Df Chisq Pr(>Chisq) 1 1236 -456.46 2 1240 -460.27 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:
# 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:
# Equal Slopes NUR für "JobFamily" ppom < - vglm(Mood ~ JobFamily+JobSatisfaction, data=mydata, family=cumulative(parallel=T~JobFamily-1));ppom # Nochmals EqualSlopes NUR für "JobFamily" ppom <- vglm(Mood ~ JobFamily+JobSatisfaction, data=mydata, family=cumulative(parallel=F~JobSatisfaction));ppom
parallel=T~JobFamily-1 bedeutet, dass equal slopes nur für JobFamily angenommen wird.
parallel=F~JobSatisfaction bedeutet, dass equal slopes nur für JobSatisfaction nicht angenommen wird.
Beide Befehle bedeuten also das selbe. Daher lautet die Ausgabe bei beiden Befehlen:
Call:
vglm(formula = Mood ~ JobFamily + JobSatisfaction, family = cumulative(parallel = F ~
JobSatisfaction), data = mydata)
Coefficients:
(Intercept):1 (Intercept):2 (Intercept):3 JobFamily JobSatisfaction:1 JobSatisfaction:2
0.6614816 1.5373699 2.7337389 0.5705897 -1.9684047 -1.2805318
JobSatisfaction:3
-0.8865225
Degrees of Freedom: 1245 Total; 1238 Residual
Residual deviance: 914.4929
Log-likelihood: -457.2464
Eine Koeffizientenübersicht erhält man per
ppom.ce < - coef(summary(ppom));ppom.ce
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.5705897 0.3801667 0.7610126 0.09715456 5.873009 4.279541e-09 1.76931
