[GNU R]: ordinale Regression

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

Links und Literatur

This entry was posted in Professor Hastig and tagged . Bookmark the permalink. Follow any comments here with the RSS feed for this post. Post a comment or leave a trackback.

Schreibe einen Kommentar

Deine E-Mail-Adresse wird nicht veröffentlicht. Erforderliche Felder sind mit * markiert.

Your email address will never be published.

Diese Website verwendet Akismet, um Spam zu reduzieren. Erfahre mehr darüber, wie deine Kommentardaten verarbeitet werden.