Produnis uninteressant – aber wahr

24. November 2013

[GNU R]: ordinale Regression

Filed under: Professor Hastig — Schlagwörter: — produnis @ 11:35

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

20. Juli 2013

[GNU R]: Logistische Regression mit R

Filed under: Professor Hastig — Schlagwörter: — produnis @ 12:54

Als Beispiel dienen Daten der Challenger-Katastrophe [Am 28. Januar 1986 explodierte die
Raumfähre Challenger 73 Sekunden nach dem Start. Der Ausfall eines oder mehrerer Dichtungsringe in einer der seitlichen Feststoffrakten wurde als Grund der Katastrophe ermittelt.]

Probleme mit den Dichtungsringen sind auch vorher am Boden aufgetreten. So liegen Erhebungsdaten vor, die den Defekt eines Rings in Zusammenhang mit der Außentemperatur darstellen können: Je kälter es draußen ist, desto höher die Ausfallwahrscheinlichkeit eines Rings. Bedauerlicher Weise wurden diese Daten erst rückwirkend ausgewertet.

Diese Daten übertragen wir wie folgt in R:

Temp = c(66,67,68,70,72,75,76,79,53,58,70,75,67,67,69,70,73,76,78,81,57,63,70)
Ausfall = factor(c(0,0,0,0,0,0,0,0,1,1,1,1,0,0,0,0,0,0,0,0,1,1,1))
shuttle = data.frame(Temp, Ausfall)
colnames(shuttle) = c("Temp", "Ausfall")

Die Kodierung des factor Ausfall lautet: 0=kein Ausfall | 1=Ausfall
Die Einheit für Temp lautet Fahrenheit.

Der Aufruf zur logistischen Regression [Ausfall erklärt durch Temperatur] lautet:

fit = glm(Ausfall~Temp,data=shuttle, family=binomial); fit

Dies liefert folgenden Output:

Call:  glm(formula = Ausfall ~ Temp, family = binomial, data = shuttle)

Coefficients:
(Intercept)         Temp  
    15.0429      -0.2322  

Degrees of Freedom: 22 Total (i.e. Null);  21 Residual
Null Deviance:	    28.27 
Residual Deviance: 20.32 	AIC: 24.32 

Die Odds Ratio (OR) errechnet sich wie folgt:

exp(-0.2322)
[1] 0.7927875

Dies liefert das Ergebnis OR = 0.7927875. Die Interpretation lautet: Mit jedem höhren Grad an Temparatur ändert sich die Chance auf einem Ausfall um den Faktor 0.7927. Das bedeutet, dass die Chancen eines Ausfalls niedriger werden, je höher die Außentemperatur steigt.
In der Nacht vor dem Challengerstart konnte eine Temperatur von 37°Fahrenheit gemessen werden. Auf Grundlage der logistischen Regression lässt sich nun die Ausfallwahrscheinlichkeit bei 37° Außentemperatur berechnen

# zuerst eine Hilfsfunktion zum umrechnen
y = function(x){
     blup = exp(15.0429 + (x*(-0.2322)))
     bla = blup/(1+blup)
     return(bla)}
# jetzt setzen wir 37 ein:
y(37)
[1] 0.9984243

Das Modell sagt bei einer Außentemperatur von 37° eine Ausfallwahrscheinlichkeit von 99,84% voraus.

Die Funktion hätten wir uns auch sparen können, da der Befehl predict diese Arbeit übernimmt.

predict(fit, data.frame(Temp=37),type="response")

Jetzt kann man dies für jede mögliche Temperatur durchmachen, und schon erhält man dieses schöne Plot.

y = predict(fit, data.frame(Temp=30:90),type="response")
plot(30:90, y, type="l", ylab="Ausfallswahrscheinlichkeit", xlab="Temperatur in F")
abline(v=seq(20,100,1),h=seq(-0.2,1.2,0.05),col="gray") #Mathegitter
abline(v=37,col="red") # 37-Grad
abline(h=0.5,col="black") # Mittellinie
points(30:90, y, type="l",lwd=2,col="blue") # Schätzwerte in "blau"

logReg-Shuttle1

24. Mai 2013

Häufigkeiten von factor-levels für Gruppen plotten

Filed under: dev/null — Schlagwörter: , , — produnis @ 09:29

Gegeben ist folgende Tabelle:

JAHR LAND
1999 UK
1999 UK
1999 UK
1999 UK
2000 UK
2000 UK
2001 UK
2001 UK
2001 UK
2001 UK
2003 UK
2003 UK
2003 UK
2003 UK
2004 UK
2005 UK
2005 UK
2005 UK
2005 UK
2006 UK
2006 UK
2008 UK
2008 UK
2008 UK
2008 UK
2008 UK
2009 UK
2001 GER
2002 GER
2003 GER
2003 GER
2004 GER
2004 GER
2005 GER
2005 GER
2006 GER
2006 GER
2006 GER
2007 GER
2008 GER
2009 GER
2010 GER
2011 GER
2011 GER
2011 GER
2006 UK
2011 GER
2012 GER
2012 GER
2012 GER
2005 USA
2007 USA
2007 USA
2009 USA
2010 USA
2010 USA
2011 USA
2012 USA
2012 USA

Jetzt möchte ich mit ggplot eine Übersichtsgraphik erstellen. Hierbei soll die X-Achse die Jahreszahlen wiederspiegeln. Auf der Y-Achse sollen die Häufigkeiten der einzelnen Länder pro Jahr angezeigt werden.

Da die Länder-Variable ein factor ist, benötige ich einen ggplot-Code, welcher die Häufigkeiten der einzelnen factor-levels (Länder) gruppiert nach Jahren ausrechnet…

Das funktioniert hier „manuell“ per

# Datensatz einlesen
dat < - read.table("http://www.produnis.de/R/Faktorfrequenz.txt", sep="\t", colClasses = "character", header=T)

# Berechnung
with(dat,table(JAHR,LAND))

Aber wie "übersetze" ich das in einen ggplot-Code? Da ggplot selbst statistische Berechnungen durchführen kann, läge der Vorteil eines "reinen" ggplot-Codes darin, dass neue factor-levels (in diesem Fall also ein neues LAND) direkt in der Grafik eingebunden wäre, ohne dass man den Code ändern müsste.

ggplot-Code

Und so funktioniert es:

# Datensatz einlesen
dat < - read.table("http://www.produnis.de/R/Faktorfrequenz.txt", sep="\t", colClasses = "character", header=T)

# plotten
#------------
# Variante 1
ggplot(df, aes(x=factor(JAHR), fill=LAND)) + 
  geom_bar(position="dodge")

FreqPlot1

# Variante 2
ggplot(df, aes(JAHR)) + 
  geom_freqpoly(aes(y= ..count.., fill=LAND, group=LAND, colour=LAND),stat="bin",binwidth=1)

FreqPlot2

# Variante 3
ggplot(df, aes(JAHR)) + 
  geom_area(aes(y= ..count.., fill=LAND, group=LAND),stat="bin",binwidth=1)

FreqPlot3

Weblinks

Siehe auch

Dieser Artikel ist Teil des Handbuchs

22. Mai 2013

Ich mache XY mit Ubuntu

Filed under: Jean Pütz,Ubuntu,ubuntuusers — Schlagwörter: , , , — produnis @ 20:36

Ubuntu

  1. System
  2. Netzwerk
  3. Multimedia
  4. Software

LaTeX

GNU R

  1. Grafiken

Powered by WordPress