Czasami interesuje nas efekt wpływu kilku zmiennych wyjaśniających \(X_i\) (predyktorów) na zmienną wynikową \(Y\). Jeżeli uważamy, że zmienne te mogą być wewnętrznie powiązane, nie powinniśmy osobno obserwować wpływu na \(Y\) zmieniających się wartości pojedynczego \(X\) lecz jednocześnie wziąć pod uwagę wartości pozostałych \(X\)-ów. Funkcja lm z pakietu R działa również w takim przypadku.
# Przykład 8.1
n <- 100
x1 <- rnorm(n, 175, 7)
x2 <- rnorm(n, 30, 8)
x3 <- abs(rnorm(n, 60, 30))
y1 <- 0.2*x1 - 0.3*x2 - 0.4*x3 + 10 + rnorm(n, 0, 5)
y2 <- -0.3*x2 + 0.2*x3 + rnorm(n, 10, 2)
data <- data.frame(x1, x2, x3, y1, y2)
y.lm <- lm(cbind(y1, y2) ~ x1 + x2 + x3, data=data)
summary(y.lm)## Response y1 :
##
## Call:
## lm(formula = y1 ~ x1 + x2 + x3, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.0325 -3.4906 -0.4393 3.5103 12.2894
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.52101 13.25460 0.190 0.849555
## x1 0.26061 0.07285 3.578 0.000546 ***
## x2 -0.34687 0.06375 -5.441 4.05e-07 ***
## x3 -0.43144 0.01917 -22.502 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.026 on 96 degrees of freedom
## Multiple R-squared: 0.8594, Adjusted R-squared: 0.855
## F-statistic: 195.6 on 3 and 96 DF, p-value: < 2.2e-16
##
##
## Response y2 :
##
## Call:
## lm(formula = y2 ~ x1 + x2 + x3, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.4987 -1.5225 0.3341 1.5403 5.3096
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.171548 5.597555 1.460 0.148
## x1 0.007990 0.030763 0.260 0.796
## x2 -0.260577 0.026923 -9.679 7.3e-16 ***
## x3 0.193215 0.008097 23.862 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.123 on 96 degrees of freedom
## Multiple R-squared: 0.8707, Adjusted R-squared: 0.8667
## F-statistic: 215.5 on 3 and 96 DF, p-value: < 2.2e-16Regresja logistyczna jest podobna do regresji liniowej, jednakże w tym przypadku zmienna wynikowa \(Y\) ma postać binarną. Ze względu na fakt, że regresja liniowa zakłada normalność zmiennej wynikowej, musimy dokonać pewnej transformacji zwanej logistyczną (lub logitową). Zamiast przewidywać tylko jedną z dwóch kategorii (np. "chory"=1, "zdrowy"=0), rozważamy prawdopodobieństwo \(p\) tego, że osobnik został zakwalifikowany do jednej z dwóch kategorii. Po takiej transformacji zmienną wynikową \(Z\) jest tzw. logit, czyli logarytm naturalny szansy wystąpienia wydarzenia:
\(Z = \ln(\frac{p}{1-p}) = a + b_1x_1 + b_2x_2+\dots+b_nx_n\).
W pakiecie R regresję logistyczną można wykonać przy pomocy ogólniejszej funkcji glm (Generalized Linear Models). Przy wywołaniu funkcji należy podać parametr family = binomial a zmienną wynikową można przekazać na jeden z trzech sposobów:
factor określająca sukces,weights całkowitą liczbę przypadków),Rozpatrzmy zbiór danych Titanic z pakietu datasets. W zbiorze tym znajduje się 4-wymiarowa tablica zawierająca informacje o 2201 pasażerach feralnego rejsu.
# Przykład 8.2
titanic <- Titanic; titanic## , , Age = Child, Survived = No
##
## Sex
## Class Male Female
## 1st 0 0
## 2nd 0 0
## 3rd 35 17
## Crew 0 0
##
## , , Age = Adult, Survived = No
##
## Sex
## Class Male Female
## 1st 118 4
## 2nd 154 13
## 3rd 387 89
## Crew 670 3
##
## , , Age = Child, Survived = Yes
##
## Sex
## Class Male Female
## 1st 5 1
## 2nd 11 13
## 3rd 13 14
## Crew 0 0
##
## , , Age = Adult, Survived = Yes
##
## Sex
## Class Male Female
## 1st 57 140
## 2nd 14 80
## 3rd 75 76
## Crew 192 20# Liczba osób, którzy przeżyli z podziałem na klasę, płeć i wiek
titanic.df <- as.data.frame(titanic[,,,"Yes"])
titanic.df## Class Sex Age Freq
## 1 1st Male Child 5
## 2 2nd Male Child 11
## 3 3rd Male Child 13
## 4 Crew Male Child 0
## 5 1st Female Child 1
## 6 2nd Female Child 13
## 7 3rd Female Child 14
## 8 Crew Female Child 0
## 9 1st Male Adult 57
## 10 2nd Male Adult 14
## 11 3rd Male Adult 75
## 12 Crew Male Adult 192
## 13 1st Female Adult 140
## 14 2nd Female Adult 80
## 15 3rd Female Adult 76
## 16 Crew Female Adult 20# Dodanie kolumny, z liczbą osób, które nie przeżyły i zmiana nazwy czwartej kolumny
titanic.df$Died <- as.data.frame(titanic[,,,"No"])$Freq
names(titanic.df)[4] <- "Survived"
titanic.df## Class Sex Age Survived Died
## 1 1st Male Child 5 0
## 2 2nd Male Child 11 0
## 3 3rd Male Child 13 35
## 4 Crew Male Child 0 0
## 5 1st Female Child 1 0
## 6 2nd Female Child 13 0
## 7 3rd Female Child 14 17
## 8 Crew Female Child 0 0
## 9 1st Male Adult 57 118
## 10 2nd Male Adult 14 154
## 11 3rd Male Adult 75 387
## 12 Crew Male Adult 192 670
## 13 1st Female Adult 140 4
## 14 2nd Female Adult 80 13
## 15 3rd Female Adult 76 89
## 16 Crew Female Adult 20 3# Wykonanie regresji logistycznej
titanic.lr <- glm(cbind(Survived,Died)~Class+Sex+Age, data=titanic.df, family=binomial(link = "logit"))
summary(titanic.lr)##
## Call:
## glm(formula = cbind(Survived, Died) ~ Class + Sex + Age, family = binomial(link = "logit"),
## data = titanic.df)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -4.1356 -0.7004 0.3039 2.2401 4.3833
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.6853 0.2730 2.510 0.0121 *
## Class2nd -1.0181 0.1960 -5.194 2.05e-07 ***
## Class3rd -1.7778 0.1716 -10.362 < 2e-16 ***
## ClassCrew -0.8577 0.1573 -5.451 5.00e-08 ***
## SexFemale 2.4201 0.1404 17.236 < 2e-16 ***
## AgeAdult -1.0615 0.2440 -4.350 1.36e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 671.96 on 13 degrees of freedom
## Residual deviance: 112.57 on 8 degrees of freedom
## AIC: 171.19
##
## Number of Fisher Scoring iterations: 5Z analizy wynika, że wszystkie zmienne są statystycznie istotne. Wyniki wskazują, że kobiety na statku miały 11.25 raza większą szansę (exp(2.4201)=11.25), a pasażerowie trzeciej klasy 5.92 razy mniejszą szansę na przeżycie. Co więcej, w modelu regresji logistycznej następuje efekt multiplikacji, co oznacza, że mężczyźni z trzeciej klasy mieli 66.5 raza mniejszą szansę na przeżycie od pozostałych pasażerów.
Aby wykorzystać otrzymany model do przewidywania należy użyć funkcji predict
# Przykład 8.3
test <- data.frame(Class=c("1st","2nd","3rd"), Age=rep("Adult",3), Sex=rep("Male",3))
predict(titanic.lr, newdata = test, type = "link")## 1 2 3
## -0.3762229 -1.3943179 -2.1539851predict(titanic.lr, newdata = test, type = "response")## 1 2 3
## 0.4070382 0.1987193 0.1039594W przypadku gdy zmienna wynikowa jest częstością wystąpienia jakiegoś zdarzenia (lub po prostu liczbą wystąpień w jakimś okresie) możemy posłużyć się regresją Poissona w celu znalezienia jej zależności od zmiennych wyjaśniających \(x_i\). Przykładem może być częstotliwość napadów padaczkowych w zależności od różnych czynników (przejście określonych chorób, dieta, przyjmowanie określonych leków). Model regresji Poissona przybiera podobną formę do modelu regresji logistycznej:
\(\ln(r) = a + b_1x_1 + b_2x_2+\dots+b_nx_n\),
gdzie \(r\) jest oczekiwaną częstością lub liczbą zliczeń wystąpień dla osobników z określonym zestawem wartości \(x_1,\dots,x_n\), a \(b_1,\dots,b_n\) są szukanymi współczynnikami regresji. Eksponenta poszczególnych współczynników są oszacowanymi częstościami względnymi powiązanymi z odpowiednimi zmiennymi.
W pakiecie R regresję Poissona można wykonać przy pomocy funkcji glm (Generalized Linear Models) z parametrem family = poisson. Jako przykład rozważmy zbiór danych warpbreaks zawierający informacje o liczbie pęknięć osnowy na krosno w zależności od rodzaju wełny i napięcia.
# Przykład 8.4
warp <- warpbreaks
warp## breaks wool tension
## 1 26 A L
## 2 30 A L
## 3 54 A L
## 4 25 A L
## 5 70 A L
## 6 52 A L
## 7 51 A L
## 8 26 A L
## 9 67 A L
## 10 18 A M
## 11 21 A M
## 12 29 A M
## 13 17 A M
## 14 12 A M
## 15 18 A M
## 16 35 A M
## 17 30 A M
## 18 36 A M
## 19 36 A H
## 20 21 A H
## 21 24 A H
## 22 18 A H
## 23 10 A H
## 24 43 A H
## 25 28 A H
## 26 15 A H
## 27 26 A H
## 28 27 B L
## 29 14 B L
## 30 29 B L
## 31 19 B L
## 32 29 B L
## 33 31 B L
## 34 41 B L
## 35 20 B L
## 36 44 B L
## 37 42 B M
## 38 26 B M
## 39 19 B M
## 40 16 B M
## 41 39 B M
## 42 28 B M
## 43 21 B M
## 44 39 B M
## 45 29 B M
## 46 20 B H
## 47 21 B H
## 48 24 B H
## 49 17 B H
## 50 13 B H
## 51 15 B H
## 52 15 B H
## 53 16 B H
## 54 28 B H# Regresja Poissona
warp.pr <- glm(breaks~wool+tension, data=warp, family = poisson(link = "log"))
summary(warp.pr)##
## Call:
## glm(formula = breaks ~ wool + tension, family = poisson(link = "log"),
## data = warp)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.6871 -1.6503 -0.4269 1.1902 4.2616
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.69196 0.04541 81.302 < 2e-16 ***
## woolB -0.20599 0.05157 -3.994 6.49e-05 ***
## tensionM -0.32132 0.06027 -5.332 9.73e-08 ***
## tensionH -0.51849 0.06396 -8.107 5.21e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 297.37 on 53 degrees of freedom
## Residual deviance: 210.39 on 50 degrees of freedom
## AIC: 493.06
##
## Number of Fisher Scoring iterations: 4Powyższe wyniki wskazują na to, że wełna typu B powoduje względnie 0.814 (exp(-0.206)=0.814) raza mniej uszkodzeń osnowy niż wełna typu A. Innymi słowy zmiana wełny z A na B sprawi, że liczba uszkodzeń spadnie o 18.6%. Sprawdźmy jaką skuteczność predykcyjną ma zbudowany przez nas model. W tym celu podzielimy dane na dwa zbiory: zbiór treningowy i testowy
# Przykład 8.5
# Predykcja dla nowych danych
warp.tr <- warpbreaks[1:44,]
warp.te <- warpbreaks[45:54,]
warp.pr <- glm(breaks~wool+tension, data=warp.tr, family = poisson(link = "log"))
warp.predict <- predict(warp.pr, newdata=warp.te[,-1], type="response")
warp.error <- mean(abs(warp.te[,1]-warp.predict))
warp.error## [1] 4.305155Pima.tr i Pima.te z biblioteki MASS i zapoznaj się z nim. Dla zbioru Pima.tr wykonaj model regresji logistycznej i określ które cechy są statystycznie istotne przy określaniu cukrzycy. Następnie sprawdź skuteczność predykcyjną modelu przy użyciu zbioru Pima.te. Niech metoda predict zwraca prawdopodobieństwo bycia chorym na cukrzyce, które następnie przy pomocy funkcji ifelse należy zamienić na odpowiedź binarną (chory=1, zdrowy=0). Jako próg prawdopodobieństwa powyżej którego uznajemy osobnika za chorego przyjąć 0.5. Oblicz skuteczność modelu według poniższego wzoru.
pima.accuracy <- 1-mean(abs(pima.true-pima.result))
pima.accuracy## [1] 0.8012048