GNU R: Ordinale Regression mit R

Professor Hastig
R
Autor:in

produnis

Veröffentlichungsdatum

12. Juni 2015

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)

Coefficients: 
              Estimate Std. Error z value Pr(>|z|)    
(Intercept):1  0.62102    0.51277   1.211   0.2258    
(Intercept):2 -1.77021    0.51586  -3.432   0.0006 ***
(Intercept):3 -4.05769    0.54042  -7.508 5.99e-14 ***
Konflikt      -0.57611    0.07627  -7.553 4.24e-14 ***
Zufriedenh     1.36420    0.16034   8.508  < 2e-16 ***
---
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: 11 

No Hauck-Donner effect found in any of the estimates


Exponentiated coefficients:
  Konflikt Zufriedenh 
 0.5620803  3.9126029

Mit dem Modell können nun weitere Parameter errechnet werden:

# Koeffizienten
pom.coef <- (coef(summary(pom)));pom.coef
                Estimate Std. Error   z value     Pr(>|z|)
(Intercept):1  0.6210245 0.51276900  1.211119 2.258496e-01
(Intercept):2 -1.7702089 0.51585658 -3.431591 6.000515e-04
(Intercept):3 -4.0576853 0.54042146 -7.508372 5.986725e-14
Konflikt      -0.5761105 0.07627287 -7.553282 4.244248e-14
Zufriedenh     1.3642029 0.16033569  8.508417 1.763242e-17
# Odds Ratio
pom.odds <- exp(coef(pom));pom.odds
(Intercept):1 (Intercept):2 (Intercept):3      Konflikt    Zufriedenh 
   1.86083349    0.17029742    0.01728899    0.56208033    3.91260294 
# Devianz und AIC
pom.devi <- deviance(pom);pom.devi
[1] 920.5493
pom.aic  <- AIC(pom);pom.aic
[1] 930.5493
# logLikelihood
pom.ll <- logLik(pom);pom.ll
[1] -460.2746
# 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
# 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
# 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.

# 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:

# 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 "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.6612053     1.5371932     2.7332281     0.5705813    -1.9683001 
 Zufriedenh:2  Zufriedenh:3 
   -1.2804632    -0.8863414 

Degrees of Freedom: 1245 Total; 1238 Residual
Residual deviance: 914.4929 
Log-likelihood: -457.2464 
# 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.6612053     1.5371932     2.7332281     0.5705813    -1.9683001 
 Zufriedenh:2  Zufriedenh:3 
   -1.2804632    -0.8863414 

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

ppom.ce <- coef(summary(ppom));ppom.ce
                Estimate Std. Error    z value     Pr(>|z|)
(Intercept):1  0.6612053 0.78638478  0.8408165 4.004507e-01
(Intercept):2  1.5371932 0.61812089  2.4868812 1.288684e-02
(Intercept):3  2.7332281 0.89023625  3.0702278 2.138956e-03
Konflikt       0.5705813 0.07607737  7.5000132 6.381141e-14
Zufriedenh:1  -1.9683001 0.33252407 -5.9192710 3.233717e-09
Zufriedenh:2  -1.2804632 0.21301487 -6.0111449 1.842177e-09
Zufriedenh:3  -0.8863414 0.30450746 -2.9107380 3.605763e-03

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.5705813 0.4214696 0.7196929 0.07607737 7.500013 6.37268e-14 1.769295

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)

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.

# 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:

# 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
# 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:

# 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
# Devianz
crm.devi <- deviance(crm); crm.devi
[1] 920.4627
# 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

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