Ordinale Regression mit R

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:

library(VGAM)
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
## [1] 920.5

# AIC
pom.aic <- AIC(pom)
pom.aic
## [1] 930.5
# logLikelihood
pom.ll <- logLik(pom)
pom.ll
## [1] -460.3
# 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.1241
# 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.2696
# 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.2929

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
## [1] 0.1065

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

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
## [1] 0.1241

# R^2 Cox&Snell
crm.cox <- as.vector(1 - exp((2/N) * (c0.ll - crm.ll)))
crm.cox
## [1] 0.2698

# 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.293

# Devianz
crm.devi <- deviance(crm)
crm.devi
## [1] 920.5

# AIC
crm.aic <- AIC(crm)
crm.aic
## [1] 930.5

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 und Literatur