Wie mache ich eine Regression mit Effektcodierung anstelle von Dummy-Codierung in R?

14

Ich arbeite derzeit an einem Regressionsmodell, bei dem ich nur kategoriale / Faktor-Variablen als unabhängige Variablen habe. Meine abhängige Variable ist ein logit transformiertes Verhältnis.

Es ist ziemlich einfach, eine normale Regression in R auszuführen, da R automatisch weiß, wie Dummies codiert werden, sobald sie vom Typ "Faktor" sind. Diese Art der Codierung impliziert jedoch auch, dass eine Kategorie aus jeder Variablen als Basis verwendet wird, was die Interpretation erschwert.

Mein Professor hat mir geraten, stattdessen nur die Effektcodierung (-1 oder 1) zu verwenden, da dies die Verwendung des großen Mittels für den Achsenabschnitt impliziert.

Weiß jemand, wie man damit umgeht?

Bis jetzt habe ich versucht:

gm <- mean(tapply(ds$ln.crea, ds$month,  mean))
model <- lm(ln.crea ~ month + month*month + year + year*year, data = ds, contrasts = list(gm = contr.sum))

Call:
lm(formula = ln.crea ~ month + month * month + year + year * 
    year, data = ds, contrasts = list(gm = contr.sum))

Residuals:
     Min       1Q   Median       3Q      Max 
-0.89483 -0.19239 -0.03651  0.14955  0.89671 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -3.244493   0.204502 -15.865   <2e-16 ***
monthFeb    -0.124035   0.144604  -0.858   0.3928    
monthMar    -0.365223   0.144604  -2.526   0.0129 *  
monthApr    -0.240314   0.144604  -1.662   0.0993 .  
monthMay    -0.109138   0.144604  -0.755   0.4520    
monthJun    -0.350185   0.144604  -2.422   0.0170 *  
monthJul     0.050518   0.144604   0.349   0.7275    
monthAug    -0.206436   0.144604  -1.428   0.1562    
monthSep    -0.134197   0.142327  -0.943   0.3478    
monthOct    -0.178182   0.142327  -1.252   0.2132    
monthNov    -0.119126   0.142327  -0.837   0.4044    
monthDec    -0.147681   0.142327  -1.038   0.3017    
year1999     0.482988   0.200196   2.413   0.0174 *  
year2000    -0.018540   0.200196  -0.093   0.9264    
year2001    -0.166511   0.200196  -0.832   0.4073    
year2002    -0.056698   0.200196  -0.283   0.7775    
year2003    -0.173219   0.200196  -0.865   0.3887    
year2004     0.013831   0.200196   0.069   0.9450    
year2005     0.007362   0.200196   0.037   0.9707    
year2006    -0.281472   0.200196  -1.406   0.1625    
year2007    -0.266659   0.200196  -1.332   0.1855    
year2008    -0.248883   0.200196  -1.243   0.2164    
year2009    -0.153083   0.200196  -0.765   0.4461    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Residual standard error: 0.3391 on 113 degrees of freedom
Multiple R-squared: 0.3626, Adjusted R-squared: 0.2385 
F-statistic: 2.922 on 22 and 113 DF,  p-value: 0.0001131 
Kasper Christensen
quelle
1
Schauen Sie sich? Kontraste Ich denke, es ist kontr.sum, um gegen den großen Mittelwert zu testen - überprüfen Sie die R-Hilfedateien
user20650
2
Dies könnte hilfreich sein: unc.edu/courses/2006spring/ecol/145/001/docs/lectures/…
mark999

Antworten:

13

Grundsätzlich gibt es zwei Arten der Kontrastcodierung, mit denen der Achsenabschnitt den Grand Mean abschätzt. Dies sind Summenkontraste und wiederholte Kontraste (gleitende Unterschiede).

Hier ist ein Beispieldatensatz:

set.seed(42)
x <- data.frame(a = c(rnorm(100,2), rnorm(100,1),rnorm(100,0)),
                b = rep(c("A", "B", "C"), each = 100))

Die Bedingungen bedeuten:

tapply(x$a, x$b, mean)
         A           B           C 
2.03251482  0.91251629 -0.01036817 

Der große Mittelwert:

mean(tapply(x$a, x$b, mean))
[1] 0.978221

Sie können die Art der Kontrastcodierung mit dem contrastsParameter in angeben lm.

Summenkontraste

lm(a ~ b, x, contrasts = list(b = contr.sum))

Coefficients:
(Intercept)           b1           b2  
     0.9782       1.0543      -0.0657 

Der Achsenabschnitt ist der Grand Mean. Die erste Steigung ist die Differenz zwischen der ersten Faktorstufe und dem Grand Mean. Die zweite Steigung ist die Differenz zwischen der zweiten Faktorstufe und dem Grand Mean.

Wiederholte Kontraste

Die Funktion zum Erzeugen wiederholter Kontraste ist Bestandteil des MASSPakets.

lm(a ~ b, x, contrasts = list(b = MASS::contr.sdif))

Coefficients:
(Intercept)         b2-1         b3-2  
     0.9782      -1.1200      -0.9229 

Der Achsenabschnitt ist der Grand Mean. Die Steigungen zeigen die Unterschiede zwischen aufeinanderfolgenden Faktorstufen (2 gegen 1, 3 gegen 2).

Sven Hohenstein
quelle
Hmm, habe gerade ausprobiert, was Sie vorgeschlagen haben, aber ich bin nicht sicher, ob einer der Codes das tut, was ich will. Die Sache ist, dass ich Jahre {1998, ..., 2007} in einer IV und Monate {Jan, ..., Dec} in einer anderen IV habe. Da es sich nun um die lm-Funktion handelt, wird der April automatisch zum Abfang sowie 1998. Stattdessen möchte ich nur, dass der Abfang ein allgemeiner Mittelwert ist ... Ich weiß nicht wirklich, ob er sinnvoll ist, wenn ich darüber nachdenke ...
Kasper Christensen
@KasperChristensen Wenn Sie die Kontraste wie in den Beispielen angeben, ist der Achsenabschnitt der große Mittelwert. Bitte geben Sie ein reproduzierbares Beispiel dafür, was Sie versucht haben.
Sven Hohenstein
@SvenHohenstein warum gibt es keinen b3-Koeffizienten für C-Kategoriewert in Summenkontrasten? Es sollte -0,9885891 sein.
Vivaldi
@Vivaldi Der Wert von b3 wird durch den Achsenabschnitt und b1, b2 bestimmt. Es gibt keinen Freiheitsgrad mehr für einen anderen Kontrast.
Sven Hohenstein
@SvenHohenstein Ist dies nicht eher ein Kollinearitätsproblem, da b3 direkt als lineare Kombination anderer Variablen ausgedrückt werden kann: (3 * Mittelwert - b1 - b2)?
Vivaldi
6

Nitpicking: Wenn Sie von Ihrem Professor angewiesen wurden, Ihre Variablen mit (-1, 1)zu codieren , sollten Sie die Effektcodierung und nicht die Effektgrößen verwenden . @ User20650 ist jedenfalls richtig. Wie üblich bietet die UCLA-Statistikhilfe-Website eine nützliche Seite, auf der erklärt wird, wie dies mit R geschehen kann.

gung - Wiedereinsetzung von Monica
quelle