R Language
Modele liniowe (regresja)

Składnia

• lm (formuła, dane, podzbiór, wagi, na.action, method = "qr", model = TRUE, x = FALSE, y = FALSE, qr = TRUE, liczba pojedyncza .ok = PRAWDA, kontrasty = NULL, przesunięcie, .. .)

Parametry

Regresja liniowa w zestawie danych mtcars

Wbudowana ramka danych mtcars zawiera informacje o 32 samochodach, w tym o masie, zużyciu paliwa (w milach na galon), prędkości itp. (Aby dowiedzieć się więcej o zestawie danych, skorzystaj z help(mtcars) ).

Jeśli interesuje nas związek między efektywnością zużycia paliwa ( mpg ) a wagą ( wt ), możemy zacząć rysować te zmienne za pomocą:

plot(mpg ~ wt, data = mtcars, col=2)

Wykresy pokazują (liniowy) związek! Następnie, jeśli chcemy wykonać regresję liniową w celu ustalenia współczynników modelu liniowego, lm funkcji lm :

fit <- lm(mpg ~ wt, data = mtcars)

~ Oznacza tutaj „wyjaśnione przez”, więc wzór mpg ~ wt oznacza, że przewidujemy mpg, jak wyjaśniono przez wt. Najbardziej pomocnym sposobem przeglądania wyników jest:

summary(fit)

Co daje wynik:

Call:
lm(formula = mpg ~ wt, data = mtcars)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.5432 -2.3647 -0.1252  1.4096  6.8727 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  37.2851     1.8776  19.858  < 2e-16 ***
wt           -5.3445     0.5591  -9.559 1.29e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.046 on 30 degrees of freedom
Multiple R-squared:  0.7528,    Adjusted R-squared:  0.7446 
F-statistic: 91.38 on 1 and 30 DF,  p-value: 1.294e-10

Zawiera informacje o:

• szacowane nachylenie każdego współczynnika ( wt i punkt przecięcia y), co sugeruje, że prognoza najlepszego dopasowania mpg wynosi 37.2851 + (-5.3445) * wt
• Wartość p każdego współczynnika, co sugeruje, że przecięcie i waga prawdopodobnie nie wynikają z przypadku
• Ogólne oszacowania dopasowania, takie jak R ^ 2 i skorygowane R ^ 2, które pokazują, jak wiele zmian w mpg wyjaśnia model

Możemy dodać linię do naszej pierwszej fabuły, aby pokazać przewidywane mpg :

abline(fit,col=3,lwd=2)

Możliwe jest również dodanie równania do tego wykresu. Najpierw uzyskaj współczynniki z coef . Następnie za pomocą paste0 współczynniki odpowiednimi zmiennymi i +/- , aby zbudować równanie. Na koniec dodajemy go do wykresu za pomocą mtext :

bs <- round(coef(fit), 3) 
lmlab <- paste0("mpg = ", bs[1],
             ifelse(sign(bs[2])==1, " + ", " - "), abs(bs[2]), " wt ")
mtext(lmlab, 3, line=-2) 

Wynik to:

Rysowanie regresji (podstawa)

Kontynuując przykład mtcars , oto prosty sposób na utworzenie wykresu regresji liniowej, który jest potencjalnie odpowiedni do publikacji.

Najpierw dopasuj model liniowy i

fit <- lm(mpg ~ wt, data = mtcars)

Następnie wykreśl dwie zmienne będące przedmiotem zainteresowania i dodaj linię regresji w domenie definicji:

plot(mtcars$wt,mtcars$mpg,pch=18, xlab = 'wt',ylab = 'mpg')
lines(c(min(mtcars$wt),max(mtcars$wt)),
as.numeric(predict(fit, data.frame(wt=c(min(mtcars$wt),max(mtcars$wt))))))

Prawie na miejscu! Ostatnim krokiem jest dodanie do wykresu, równania regresji, rsquare, a także współczynnika korelacji. Odbywa się to za pomocą funkcji vector :

rp = vector('expression',3)
rp[1] = substitute(expression(italic(y) == MYOTHERVALUE3 + MYOTHERVALUE4 %*% x), 
          list(MYOTHERVALUE3 = format(fit$coefficients[1], digits = 2),
                        MYOTHERVALUE4 = format(fit$coefficients[2], digits = 2)))[2]
rp[2] = substitute(expression(italic(R)^2 == MYVALUE), 
             list(MYVALUE = format(summary(fit)$adj.r.squared,dig=3)))[2]
rp[3] = substitute(expression(Pearson-R == MYOTHERVALUE2), 
             list(MYOTHERVALUE2 = format(cor(mtcars$wt,mtcars$mpg), digits = 2)))[2]

legend("topright", legend = rp, bty = 'n')

Pamiętaj, że możesz dodać dowolny inny parametr, taki jak RMSE, dostosowując funkcję wektorową. Wyobraź sobie, że chcesz legendę z 10 elementami. Definicja wektora byłaby następująca:

rp = vector('expression',10)

i musisz zdefiniować r[1] .... do r[10]

Oto wynik:

Ważenie

Czasami chcemy, aby model przypisywał większą wagę niektórym punktom danych lub przykładom niż innym. Jest to możliwe poprzez określenie wagi danych wejściowych podczas uczenia się modelu. Istnieją na ogół dwa rodzaje scenariuszy, w których możemy użyć niejednolitych wag na przykładach:

• Wagi analityczne: odzwierciedlają różne poziomy dokładności różnych obserwacji. Na przykład, jeśli analizuje się dane, w których każda obserwacja jest średnią wyników z obszaru geograficznego, waga analityczna jest proporcjonalna do odwrotności szacowanej wariancji. Przydatne, gdy mamy do czynienia ze średnimi w danych, ponieważ zapewniają proporcjonalną wagę, biorąc pod uwagę liczbę obserwacji. Źródło

• Wagi próbkowania (odwrotne wagi prawdopodobieństwa - IPW): technika statystyczna służąca do obliczania statystyk znormalizowanych dla populacji innej niż populacja, w której zebrano dane. Plany badań z odmienną populacją próbkowania i populacją wnioskowania docelowego (populacja docelowa) są powszechne w zastosowaniu. Przydatne w przypadku danych, w których brakuje wartości. Źródło

Funkcja lm() wykonuje ważenie analityczne. W przypadku wag próbkowania pakiet survey służy do zbudowania obiektu projektu ankiety i uruchomienia svyglm() . Domyślnie pakiet survey wykorzystuje wagi próbkowania. (UWAGA: lm() , a svyglm() z rodziny gaussian() .. Czy wszystkie generują te same szacunki punktowe, ponieważ oba rozwiązania dla współczynników minimalizując ważonych najmniejszych kwadratów Różnią się jak standardowe obliczane są błędy)

Dane testowe

data <- structure(list(lexptot = c(9.1595012302023, 9.86330744180814, 
8.92372556833205, 8.58202430280175, 10.1133857229336), progvillm = c(1L, 
1L, 1L, 1L, 0L), sexhead = c(1L, 1L, 0L, 1L, 1L), agehead = c(79L, 
43L, 52L, 48L, 35L), weight = c(1.04273509979248, 1.01139605045319, 
1.01139605045319, 1.01139605045319, 0.76305216550827)), .Names = c("lexptot", 
"progvillm", "sexhead", "agehead", "weight"), class = c("tbl_df", 
"tbl", "data.frame"), row.names = c(NA, -5L))

Wagi analityczne

lm.analytic <- lm(lexptot ~ progvillm + sexhead + agehead, 
                            data = data, weight = weight)
summary(lm.analytic)

Wynik

Call:
lm(formula = lexptot ~ progvillm + sexhead + agehead, data = data, 
    weights = weight)

Weighted Residuals:
         1          2          3          4          5 
 9.249e-02  5.823e-01  0.000e+00 -6.762e-01 -1.527e-16

Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept) 10.016054   1.744293   5.742    0.110
progvillm   -0.781204   1.344974  -0.581    0.665
sexhead      0.306742   1.040625   0.295    0.818
agehead     -0.005983   0.032024  -0.187    0.882

Residual standard error: 0.8971 on 1 degrees of freedom
Multiple R-squared:  0.467, Adjusted R-squared:  -1.132 
F-statistic: 0.2921 on 3 and 1 DF,  p-value: 0.8386

Odważniki do próbkowania (IPW)

library(survey)
data$X <- 1:nrow(data)             # Create unique id
        
# Build survey design object with unique id, ipw, and data.frame
des1 <- svydesign(id = ~X,  weights = ~weight, data = data)
        
# Run glm with survey design object
prog.lm <- svyglm(lexptot ~ progvillm + sexhead + agehead, design=des1)

Wynik

Call:
svyglm(formula = lexptot ~ progvillm + sexhead + agehead, design = des1)

Survey design:
svydesign(id = ~X, weights = ~weight, data = data)

Coefficients:
             Estimate Std. Error t value Pr(>|t|)  
(Intercept) 10.016054   0.183942  54.452   0.0117 *
progvillm   -0.781204   0.640372  -1.220   0.4371  
sexhead      0.306742   0.397089   0.772   0.5813  
agehead     -0.005983   0.014747  -0.406   0.7546  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 0.2078647)

Number of Fisher Scoring iterations: 2

Sprawdzanie nieliniowości za pomocą regresji wielomianowej

Czasami podczas pracy z regresją liniową musimy sprawdzić nieliniowość danych. Jednym ze sposobów jest dopasowanie modelu wielomianowego i sprawdzenie, czy pasuje on do danych lepiej niż model liniowy. Istnieją inne powody, takie jak teoretyczne, które wskazują na dopasowanie modelu kwadratowego lub wyższego rzędu, ponieważ uważa się, że związek zmiennych ma z natury charakter wielomianowy.

mtcars model kwadratowy do zestawu danych mtcars . Aby uzyskać model liniowy, zobacz Regresja liniowa w zestawie danych mtcars .

Najpierw tworzymy wykres rozproszenia zmiennych mpg (Mil / galon), disp (Przemieszczenie (cu.in.)) i wt (Waga (1000 funtów)). Zależność między mpg i disp wydaje się nieliniowa.

plot(mtcars[,c("mpg","disp","wt")])

Dopasowanie liniowe pokaże, że disp nie jest znaczący.

fit0 = lm(mpg ~ wt+disp, mtcars)
summary(fit0)

# Coefficients:
#            Estimate Std. Error t value Pr(>|t|)    
#(Intercept) 34.96055    2.16454  16.151 4.91e-16 ***
#wt          -3.35082    1.16413  -2.878  0.00743 ** 
#disp        -0.01773    0.00919  -1.929  0.06362 .  
#---
#Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Residual standard error: 2.917 on 29 degrees of freedom
#Multiple R-squared:  0.7809,    Adjusted R-squared:  0.7658

Następnie, aby uzyskać wynik modelu kwadratowego, dodaliśmy I(disp^2) . Nowy model wygląda lepiej, gdy patrzy się na R^2 a wszystkie zmienne są znaczące.

fit1 = lm(mpg ~ wt+disp+I(disp^2), mtcars)
summary(fit1)

# Coefficients:
#              Estimate Std. Error t value Pr(>|t|)    
#(Intercept) 41.4019837  2.4266906  17.061  2.5e-16 ***
#wt          -3.4179165  0.9545642  -3.581 0.001278 ** 
#disp        -0.0823950  0.0182460  -4.516 0.000104 ***
#I(disp^2)    0.0001277  0.0000328   3.892 0.000561 ***
#---
#Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Residual standard error: 2.391 on 28 degrees of freedom
#Multiple R-squared:  0.8578,    Adjusted R-squared:  0.8426 

Ponieważ mamy trzy zmienne, dopasowany model jest powierzchnią reprezentowaną przez:

mpg = 41.4020-3.4179*wt-0.0824*disp+0.0001277*disp^2

Innym sposobem określenia regresji wielomianowej jest użycie parametru poly z parametrem raw=TRUE , w przeciwnym razie uwzględnione zostaną wielomiany ortogonalne help(ploy) więcej informacji można znaleźć w help(ploy) ). Ten sam wynik uzyskujemy, używając:

summary(lm(mpg ~ wt+poly(disp, 2, raw=TRUE),mtcars))

Wreszcie, co jeśli musimy pokazać wykres szacowanej powierzchni? Istnieje wiele opcji tworzenia wykresów 3D w R Tutaj używamy Fit3d z pakietu p3d .

library(p3d)
Init3d(family="serif", cex = 1)
Plot3d(mpg ~ disp+wt, mtcars)
Axes3d()
Fit3d(fit1)

Ocena jakości

Po zbudowaniu modelu regresji ważne jest, aby sprawdzić wynik i zdecydować, czy model jest odpowiedni i dobrze współpracuje z dostępnymi danymi. Można tego dokonać, badając wykres resztkowy, a także inne wykresy diagnostyczne.

# fit the model
fit <- lm(mpg ~ wt, data = mtcars)
# 
par(mfrow=c(2,1))
# plot model object
plot(fit, which =1:2)

Te wykresy sprawdzają dwa założenia przyjęte podczas budowy modelu:

1. Że oczekiwana wartość przewidywanej zmiennej (w tym przypadku mpg ) jest dana przez liniową kombinację predyktorów (w tym przypadku wt ). Oczekujemy, że ta ocena będzie obiektywna. Zatem reszty powinny być wyśrodkowane wokół średniej dla wszystkich wartości predyktorów. W tym przypadku widzimy, że reszty są zwykle dodatnie na końcach i ujemne na środku, co sugeruje nieliniowy związek między zmiennymi.
2. Faktyczna przewidywana zmienna jest zwykle rozkładana wokół jej oszacowania. Tak więc reszty powinny być normalnie rozłożone. W przypadku normalnie rozłożonych danych punkty na normalnym wykresie QQ powinny leżeć na przekątnej lub blisko niej. Na końcach jest trochę pochylenia.

Korzystanie z funkcji „przewidywania”

Po zbudowaniu modelu predict jest główną funkcją do testowania na nowych danych. W naszym przykładzie wykorzystamy wbudowany zestaw danych mtcars do regresji mil na galon przed przemieszczeniem:

my_mdl <- lm(mpg ~ disp, data=mtcars)
my_mdl

Call:
lm(formula = mpg ~ disp, data = mtcars)

Coefficients:
(Intercept)         disp  
   29.59985     -0.04122

Gdybym miał nowe źródło danych z przesunięciem, widziałbym szacunkowe mile na galon.

set.seed(1234)
newdata <- sample(mtcars$disp, 5)
newdata
[1] 258.0  71.1  75.7 145.0 400.0

newdf <- data.frame(disp=newdata)
predict(my_mdl, newdf)
       1        2        3        4        5 
18.96635 26.66946 26.47987 23.62366 13.11381

Najważniejszą częścią tego procesu jest utworzenie nowej ramki danych o takich samych nazwach kolumn jak oryginalne dane. W tym przypadku oryginalne dane miały kolumnę z etykietą disp , z pewnością nazwałem nowe dane tą samą nazwą.

Uwaga

Spójrzmy na kilka typowych pułapek:

1. nieużywając data.frame w nowym obiekcie:

   predict(my_mdl, newdata)
   Error in eval(predvars, data, env) : 
      numeric 'envir' arg not of length one
   

2. niestosowanie tych samych nazw w nowej ramce danych:

   newdf2 <- data.frame(newdata)
   predict(my_mdl, newdf2)
   Error in eval(expr, envir, enclos) : object 'disp' not found
   

Precyzja

Aby sprawdzić dokładność prognozy, potrzebujesz rzeczywistych wartości y nowych danych. W tym przykładzie newdf będzie potrzebował kolumny dla „mpg” i „disp”.

newdf <- data.frame(mpg=mtcars$mpg[1:10], disp=mtcars$disp[1:10])
#     mpg  disp
# 1  21.0 160.0
# 2  21.0 160.0
# 3  22.8 108.0
# 4  21.4 258.0
# 5  18.7 360.0
# 6  18.1 225.0
# 7  14.3 360.0
# 8  24.4 146.7
# 9  22.8 140.8
# 10 19.2 167.6

p <- predict(my_mdl, newdf)

#root mean square error
sqrt(mean((p - newdf$mpg)^2, na.rm=TRUE))
[1] 2.325148