Kapitel 52 Faktoren in logistischen Regressionen
Problem: Wir haben eine Variable, die 0 und 1 codiert ist. Diese Variable brauchen wir in einer logistischen Regression als abhängige oder als unabhängige Variable. Wenn wir nun diese Variable zu einer Faktor-Variable transformieren, sind die Values 1 und 2. Aber keine Sorge: In der logistischen Regression bleibt diese Variable wie eine 0 und 1 codierte Variable.
Wir erstellen eine numerische Variable mit den Werten 0 und 1.
numeric<-sample(c(0,1), 100, replace=TRUE, prob=(c(0.3,0.7)))
Hmisc::describe(numeric)## numeric
## n missing distinct Info Sum Mean Gmd
## 100 0 2 0.642 69 0.69 0.4321Aus dieser Variable erstellen wir eine Faktor Variable. Die ursprünglichen Levels (Werte) sind 0 und 1.
yes_no_factor<-factor(numeric, levels=c(0,1), labels=c("nein", "Ja"))
levels(yes_no_factor)## [1] "nein" "Ja"R schützt uns vor Dummheiten. So können wir keinen Mittelwert mit der Faktorvariable rechnen.
mean(yes_no_factor)## Warning in mean.default(yes_no_factor): argument is not numeric or logical:
## returning NA## [1] NAManchman möchten wir jedoch mit Variablen, die 0 und 1 kodiert sind einen Mittelwert berechnen - weil dieser Mittelwert der Proportion (der Anteil) der 1-Antworten gibt. Wir könnten dies tun, indem wir die Faktorvariable zu einer numerischen Variable transformieren. Diese Umwandlung, wie im nächsten Code-Abschnitt gezeigt wird, wandelt jedoch die Variable in 1 und 2 um.
mean(as.numeric(yes_no_factor)) # simple functions do not treat the factor as intended.## [1] 1.69Wir müssten also die Variable wieder in 0 und 1 umwandeln.
yes1_no0_numeric<-as.numeric(yes_no_factor)-1
mean(yes1_no0_numeric)## [1] 0.69Take home message: Achtung wenn wir numerische Faktoren in Faktoren und zurück umwandeln.
Hier noch ein interessantes Beispiel:
Wenn wir eine 0 / 1 Variable zu einem Faktor umwandeln, werden die labels der Werte 0 und 1 bleiben, aber wenn wir diese Variable wieder zu einer numerischen Variable umwandeln, werden die Werte 1 und 2 sein. Dies geschieht, weil R Faktoren die Levels im Hintergrund von 0 / 1 auf 1 / 2 ändern. Das tiefste Level ist immer 1 und nicht 0.
numeric_0_1<-sample(c(0,1), 100, replace=TRUE, prob=c(0.5,0.5))
factor_0_1<-factor(numeric_0_1)
levels(factor_0_1)## [1] "0" "1"summarytools::freq(factor_0_1)## Frequencies
## factor_0_1
## Type: Factor
##
## Freq % Valid % Valid Cum. % Total % Total Cum.
## ----------- ------ --------- -------------- --------- --------------
## 0 53 53.00 53.00 53.00 53.00
## 1 47 47.00 100.00 47.00 100.00
## <NA> 0 0.00 100.00
## Total 100 100.00 100.00 100.00 100.00numeric_from_factor<-as.numeric(factor_0_1)
summarytools::freq(numeric_from_factor)## Frequencies
## numeric_from_factor
## Type: Numeric
##
## Freq % Valid % Valid Cum. % Total % Total Cum.
## ----------- ------ --------- -------------- --------- --------------
## 1 53 53.00 53.00 53.00 53.00
## 2 47 47.00 100.00 47.00 100.00
## <NA> 0 0.00 100.00
## Total 100 100.00 100.00 100.00 100.00psych::describe(numeric_from_factor)## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 100 1.47 0.5 1 1.46 0 1 2 1 0.12 -2.01 0.05Dafür gibt es aber eine einfache Lösunge: Wir müssen den Faktor zuerst zu einer charakter Variable umfunktionieren, und erst danach zu einer numrischen. Wir zeigen das hier in zwei Varianten:
character_from_factor_0_1<-as.character(factor_0_1)
summarytools::freq(character_from_factor_0_1)## Frequencies
## character_from_factor_0_1
## Type: Character
##
## Freq % Valid % Valid Cum. % Total % Total Cum.
## ----------- ------ --------- -------------- --------- --------------
## 0 53 53.00 53.00 53.00 53.00
## 1 47 47.00 100.00 47.00 100.00
## <NA> 0 0.00 100.00
## Total 100 100.00 100.00 100.00 100.00numeric_from_character_from_factor_0_1<-as.numeric(character_from_factor_0_1)
mean(numeric_from_character_from_factor_0_1)## [1] 0.47psych::describe(numeric_from_character_from_factor_0_1)## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 100 0.47 0.5 0 0.46 0 0 1 1 0.12 -2.01 0.05Jetzt das Ganze in einem Rutsch:
numeric_from_character_from_factor_0_1<-as.numeric(as.character(factor_0_1))
mean(numeric_from_character_from_factor_0_1)## [1] 0.47psych::describe(numeric_from_character_from_factor_0_1)## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 100 0.47 0.5 0 0.46 0 0 1 1 0.12 -2.01 0.0552.1 Faktoren in einer logistischen Regression.
Wir erstellen eine abhängige Variable Krank (0 / 1) und eine unabhängige Variable Variable Raucher (0 / 1).
Für die Simulation Dank an (Link)
id<-1:1000
Raucher<-sample(c(0,1), 100, replace = TRUE, prob=c(0.2, 0.8))
data<-data.frame(id, Raucher)
set.seed(666)
Raucher<-sample(c(0,1), 1000, replace = TRUE, prob=c(0.2, 0.8))
Alter = rnorm(1000, 60, 12)
z = 0.01*Alter + 0.2*Raucher# linear combination with a bias
pr = 1/(1+exp(-z)) # pass through an inv-logit function
Krank = rbinom(1000,1,pr) # bernoulli response variable
psych::describe(z)## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 1000 0.76 0.14 0.77 0.76 0.14 0.34 1.14 0.8 -0.21 -0.17 0psych::describe(((Alter-mean(Alter))/sd(Alter)))## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 1000 0 1 0.04 0.01 1.01 -3.34 2.94 6.29 -0.08 -0.01 0.03data<-data.frame(id,Krank, Raucher, Alter)
result_lr<-glm(Krank~Raucher+Alter, data=data, family="binomial")
jtools::summ(result_lr, confint=TRUE, exp=TRUE)| Observations | 1000 |
| Dependent variable | Krank |
| Type | Generalized linear model |
| Family | binomial |
| Link | logit |
| χ²(2) | 2.58 |
| Pseudo-R² (Cragg-Uhler) | 0.00 |
| Pseudo-R² (McFadden) | 0.00 |
| AIC | 1284.16 |
| BIC | 1298.89 |
| exp(Est.) | 2.5% | 97.5% | z val. | p | |
|---|---|---|---|---|---|
| (Intercept) | 1.23 | 0.58 | 2.59 | 0.54 | 0.59 |
| Raucher | 0.92 | 0.66 | 1.28 | -0.48 | 0.63 |
| Alter | 1.01 | 1.00 | 1.02 | 1.52 | 0.13 |
| Standard errors: MLE |
So weit so gut, jetzt transformieren wir die numerische (0 / 1 ) Raucher Variable zu einer Faktor Variable.
data<-data %>%
mutate(Raucher.factor=factor(Raucher))
result_lr2<-glm(Krank~Raucher.factor+Alter, data=data, family="binomial")
jtools::summ(result_lr2, confint=TRUE, exp=TRUE)| Observations | 1000 |
| Dependent variable | Krank |
| Type | Generalized linear model |
| Family | binomial |
| Link | logit |
| χ²(2) | 2.58 |
| Pseudo-R² (Cragg-Uhler) | 0.00 |
| Pseudo-R² (McFadden) | 0.00 |
| AIC | 1284.16 |
| BIC | 1298.89 |
| exp(Est.) | 2.5% | 97.5% | z val. | p | |
|---|---|---|---|---|---|
| (Intercept) | 1.23 | 0.58 | 2.59 | 0.54 | 0.59 |
| Raucher.factor1 | 0.92 | 0.66 | 1.28 | -0.48 | 0.63 |
| Alter | 1.01 | 1.00 | 1.02 | 1.52 | 0.13 |
| Standard errors: MLE |
Immer noch alles gut. Aber Achtung: Jetzt transformieren wir die Faktor Variable zurück zu einer numerischen Variable:
data<-data %>%
mutate(Raucher_numeric_from_factor=as.numeric(Raucher.factor))
result_lr3<-glm(Krank~Raucher_numeric_from_factor+Alter, data=data, family="binomial")
jtools::summ(result_lr3, confint=TRUE, exp=TRUE)| Observations | 1000 |
| Dependent variable | Krank |
| Type | Generalized linear model |
| Family | binomial |
| Link | logit |
| χ²(2) | 2.58 |
| Pseudo-R² (Cragg-Uhler) | 0.00 |
| Pseudo-R² (McFadden) | 0.00 |
| AIC | 1284.16 |
| BIC | 1298.89 |
| exp(Est.) | 2.5% | 97.5% | z val. | p | |
|---|---|---|---|---|---|
| (Intercept) | 1.33 | 0.53 | 3.36 | 0.61 | 0.54 |
| Raucher_numeric_from_factor | 0.92 | 0.66 | 1.28 | -0.48 | 0.63 |
| Alter | 1.01 | 1.00 | 1.02 | 1.52 | 0.13 |
| Standard errors: MLE |
Ups! Jetzt hat sich der Intercept geändert. Der Grund dafür: Die Variable Raucher war ursprünglich 0 / 1 kodiert. Der Intercept ist die Odds Krank zu sein für die Personen mit Alter 0, die gleichzeitig Nichtraucher sind. In einer Regression benutzt R die Faktor Variable immer noch als 0 / 1 (obschon im Vordergrund die Levels jetzt 1 und 2 sind). Wenn wir aus der Faktor-Variable eine numerische machen, wird aus 1 und 2 in den Levels auch 1 und 2 in der numerischen Variable.
Wir sehen das hier:
summarytools::freq(data$Raucher_numeric_from_factor)## Frequencies
## data$Raucher_numeric_from_factor
## Type: Numeric
##
## Freq % Valid % Valid Cum. % Total % Total Cum.
## ----------- ------ --------- -------------- --------- --------------
## 1 201 20.10 20.10 20.10 20.10
## 2 799 79.90 100.00 79.90 100.00
## <NA> 0 0.00 100.00
## Total 1000 100.00 100.00 100.00 100.0052.2 Bespiel mit numerischer Variable und Faktor-Variable als abhängige Variable.
Wir sehen auch im folgenden Beispiel, dass R trotz Faktor korrekterweise im Hintergrund 0 benutzt. Das Problem würde erst entstehen, wenn wir die Faktor Variable wieder direkt in eine numerische Variable umwandeln.
Hmisc::describe(data$Krank)## data$Krank
## n missing distinct Info Sum Mean Gmd
## 1000 0 2 0.672 661 0.661 0.4486res<-glm(Krank~1, data=data, family = "binomial")
res$coefficients## (Intercept)
## 0.667753752.2.1 Wir können jetzt aus dem Intercept auch die Probabilität der Krankheit, respektive die Prävalenz der Krankheit ungeachtet der anderen Variablen, berechnen.
# prevalence
exp(res$coefficients)/(1+exp(res$coefficients))## (Intercept)
## 0.661Alles in Ordnung wenn wir das mit der neu erstellten Faktor-Variable *Krank.factor” tun, genau das gleiche Resultat wird erzielt:
data<-data %>%
mutate(Krank.factor=factor(Krank))
Hmisc::describe(data$Krank.factor)## data$Krank.factor
## n missing distinct
## 1000 0 2
##
## Value 0 1
## Frequency 339 661
## Proportion 0.339 0.661res2<-glm(Krank.factor~1,data=data, family = binomial)
res2$coefficients## (Intercept)
## 0.6677537# prevalence
exp(res2$coefficients)/(1+exp(res2$coefficients)) # here, the glm commands does use the 0 ## (Intercept)
## 0.661Das Problem entsteht wirklich erst, wenn wir aus der Faktor-Variable direkt eine numerische Variable berechnen: Die folgende Version wird gar nicht funktionieren, weil die neu erstellte Variable Krank_numeric_from_factor 1 und 2 kodiert ist.
data<-data %>%
mutate(Krank_numeric_from_factor=as.numeric(Krank.factor))
Hmisc::describe(data$Krank_numeric_from_factor)## data$Krank_numeric_from_factor
## n missing distinct Info Mean Gmd
## 1000 0 2 0.672 1.661 0.4486
##
## Value 1 2
## Frequency 339 661
## Proportion 0.339 0.661Da die neue Variable Krank_numeric_from_factor 1 und 2 kodiert ist, wird die logistische Regression nicht laufen und eine Fehlermeldung erscheint. (By the way: Wenn ihr nicht wollt, dass das RMarkdown abbricht bei einem Fehler, müsst ihr in die chunk options {r, error=TRUE} eingeben.)
res3<-glm(Krank_numeric_from_factor~1,data=data, family = binomial)## Error in eval(family$initialize): y values must be 0 <= y <= 1Wie weiter oben gesehen ist die Lösung einfach: Einfach zuerst aus der Faktor-Variable eine Charakter-Variable erstellen und daraus die numerische Variable:
data<-data %>%
mutate(Krank_numeric_from_string_from_factor=as.numeric(as.character(Krank.factor)))
Hmisc::describe(data$Krank_numeric_from_string_from_factor)## data$Krank_numeric_from_string_from_factor
## n missing distinct Info Sum Mean Gmd
## 1000 0 2 0.672 661 0.661 0.4486Alles klar für die logistische Regression:
res4<-glm(Krank_numeric_from_string_from_factor~1,data=data, family = binomial)
res4$coefficients## (Intercept)
## 0.6677537# prevalence
exp(res4$coefficients)/(1+exp(res4$coefficients)) # here, the glm commands does use the 0 ## (Intercept)
## 0.661
Null ist nicht nichts.