R Language
Линейные модели (регрессия)

Синтаксис

• lm (формула, данные, подмножество, весы, na.action, method = "qr", model = TRUE, x = FALSE, y = FALSE, qr = TRUE, singular.ok = TRUE, контрасты = NULL, смещение, .. .)

параметры

Линейная регрессия по набору данных mtcars

Встроенный кадр данных mtcars содержит информацию о 32 машинах, включая их вес, топливную экономичность (в мили на галлон), скорость и т. Д. (Чтобы узнать больше о наборе данных, используйте help(mtcars) ).

Если нас интересует взаимосвязь между топливной экономичностью ( mpg ) и весом ( wt ), Мы можем начать строить эти переменные с:

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

Графики показывают (линейное) соотношение !. Тогда, если мы хотим выполнить линейную регрессию для определения коэффициентов линейной модели, мы будем использовать lm функцию:

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

Здесь ~ означает «объяснено», поэтому формула mpg ~ wt означает, что мы прогнозируем mpg, как объясняется wt. Самый полезный способ просмотра результатов:

summary(fit)

Что дает результат:

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

Это дает информацию о:

• оцененный наклон каждого коэффициента ( wt и y-перехват), который предполагает наилучшее предсказание mpg, составляет 37.2851 + (-5.3445) * wt
• P-значение каждого коэффициента, что предполагает, что перехват и вес, вероятно, не из-за случайности
• Общие оценки подгонки, такие как R ^ 2 и скорректированные R ^ 2, которые показывают, какая часть изменения в mpg объясняется моделью

Мы могли бы добавить строку к нашему первому сюжету, чтобы показать предсказанный mpg :

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

Также можно добавить уравнение к этому графику. Сначала получим коэффициенты с коэффициентом coef . Затем, используя paste0 мы paste0 коэффициенты с соответствующими переменными и +/- , чтобы построить уравнение. Наконец, мы добавляем его к сюжету с помощью 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) 

Результат:

Построение регрессии (базы)

Продолжая пример mtcars , вот простой способ создать график вашей линейной регрессии, который потенциально подходит для публикации.

Сначала применим линейную модель и

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

Затем постройте две интересующие переменные и добавьте линию регрессии в пределах области определения:

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))))))

Почти готово! Последний шаг - добавить к графику, уравнение регрессии, rsquare, а также коэффициент корреляции. Это делается с использованием 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')

Обратите внимание, что вы можете добавить любые другие параметры, такие как RMSE, путем адаптации векторной функции. Представьте, что вам нужна легенда с 10 элементами. Векторное определение будет следующим:

rp = vector('expression',10)

и вам нужно будет определить r[1] .... to r[10]

Вот результат:

утяжеление

Иногда мы хотим, чтобы модель придавала больший вес некоторым точкам данных или примерам, чем другие. Это возможно, указав вес для входных данных, изучая модель. Обычно существуют два типа сценариев, в которых мы можем использовать неравномерные веса над примерами:

• Аналитические весы: отражают различные уровни точности различных наблюдений. Например, если анализировать данные, где каждое наблюдение является средним результатом географической области, аналитический вес пропорционален обратному оценочной дисперсии. Полезно при работе со средними данными, предоставляя пропорциональный вес с учетом количества наблюдений. Источник

• Весы выборки (обратные весовые коэффициенты вероятности - IPW): статистический метод для расчета статистики, стандартизованной для населения, отличного от того, в котором собирались данные. В заявке рассматриваются проекты исследований с разрозненной популяцией выборки и населением целевого вывода (целевое население). Полезно при работе с данными, у которых отсутствуют значения. Источник

Функция lm() выполняет аналитическое взвешивание. Для выборочных весов пакет survey используется для создания объекта дизайна съемки и запуска svyglm() . По умолчанию в пакете survey используются весы для выборки. (ПРИМЕЧАНИЕ: lm() и svyglm() с семейством gaussian() будут производить одни и те же точечные оценки, потому что они оба решаются для коэффициентов, сводя к минимуму взвешенные наименьшие квадраты. Они отличаются тем, как вычисляются стандартные ошибки.)

Данные испытаний

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))

Аналитические весы

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

Выход

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

Вес пробы (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)

Выход

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

Проверка нелинейности с полиномиальной регрессией

Иногда при работе с линейной регрессией нам нужно проверить нелинейность данных. Один из способов сделать это - установить полиномиальную модель и проверить, соответствует ли она данным лучше, чем линейная модель. Существуют и другие причины, такие как теоретические, которые указывают на то, что они соответствуют квадратичной модели или модели более высокого порядка, поскольку считается, что отношение переменных по своей природе носит полиномиальный характер.

Построим квадратичную модель для набора данных mtcars . Для линейной модели см. Линейную регрессию по набору данных mtcars .

Сначала мы создаем график рассеяния переменных mpg (Miles / gallon), disp (смещение (cu.in.)) и wt (вес (1000 фунтов)). Связь между mpg и disp выглядит нелинейной.

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

Линейная посадка покажет, что disp не имеет значения.

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

Затем, чтобы получить результат квадратичной модели, мы добавили I(disp^2) . Новая модель выглядит лучше, если смотреть на R^2 и все переменные значительны.

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 

Поскольку у нас есть три переменные, подходящая модель представляет собой поверхность, представленную:

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

Другой способ указать полиномиальную регрессию - использовать poly с параметром raw=TRUE , в противном случае будут рассмотрены ортогональные многочлены (дополнительную информацию см. В help(ploy) ). Мы получаем тот же результат, используя:

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

Наконец, что, если нам нужно показать график предполагаемой поверхности? Ну, есть много вариантов сделать 3D графики в R Здесь мы используем Fit3d из пакета p3d .

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

Оценка качества

После построения регрессионной модели важно проверить результат и решить, подходит ли модель и хорошо работает с данными. Это можно сделать, изучив график остатков, а также другие диагностические графики.

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

Эти графики проверяют два предположения, которые были сделаны при построении модели:

1. То, что ожидаемое значение предсказанной переменной (в данном случае mpg ) задается линейной комбинацией предсказателей (в данном случае wt ). Мы ожидаем, что эта оценка будет беспристрастной. Таким образом, остатки должны быть сосредоточены вокруг среднего значения для всех значений предикторов. В этом случае мы видим, что остатки имеют тенденцию быть положительными на концах и отрицательными в середине, что указывает на нелинейную зависимость между переменными.
2. То, что фактическая предсказанная переменная обычно распределяется вокруг ее оценки. Таким образом, остатки должны нормально распределяться. Для нормально распределенных данных точки в нормальном графике QQ должны лежать на диагонали или близко к ней. На концах есть несколько перекосов.

Использование функции «предсказывать»

Как только модель построена, predict является основная функция тестирования с новыми данными. В нашем примере будет использоваться встроенный набор данных mtcars для регрессии миль на галлон против перемещения:

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

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

Coefficients:
(Intercept)         disp  
   29.59985     -0.04122

Если бы у меня был новый источник данных со смещением, я мог бы видеть приблизительные мили за галлон.

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

Наиболее важной частью процесса является создание нового фрейма данных с теми же именами столбцов, что и исходные данные. В этом случае исходные данные имели столбец с меткой disp , я был уверен, что вы вызовите новые данные с таким же именем.

предосторожность

Давайте рассмотрим несколько распространенных ошибок:

1. не используя data.frame в новом объекте:

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

2. не используя одинаковые имена в новом кадре данных:

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

точность

Чтобы проверить точность прогноза, вам понадобятся фактические значения y новых данных. В этом примере newdf потребуется столбец для «mpg» и «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