Coefficients estimate std error t value pr t

Оценка результатов линейной регрессии

Введение

Сегодня уже все, кто хоть немного интересуется дата майнингом, наверняка слышали про простую линейную регрессию. Про нее уже писали на хабре, а также подробно рассказывал Эндрю Нг в своем известном курсе машинного обучения. Линейная регрессия является одним из базовых и самых простых методов машинного обучения, однако очень редко упоминаются методы оценки качества построенной модели. В этой статье я постараюсь немного исправить это досадное упущение на примере разбора результатов функции summary.lm() в языке R. При этом я постараюсь предоставить необходимые формулы, таким образом все вычисления можно легко запрограммировать на любом другом языке. Эта статья предназначена для тех, кто слышал о том, что можно строить линейную регрессию, но не сталкивался со статистическими процедурами для оценки ее качества.

Модель линейной регрессии

Итак, пусть есть несколько независимых случайных величин X1, X2, …, Xn (предикторов) и зависящая от них величина Y (предполагается, что все необходимые преобразования предикторов уже сделаны). Более того, мы предполагаем, что зависимость линейная, а ошибки рапределены нормально, то есть

где I — единичная квадратная матрица размера n x n.

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

где b с крышкой — оценка вектора коэффициентов, y — вектор значений зависимой величины, а X — матрица размера k x n+1 (n — количество предикторов, k — количество наблюдений), у которой первый столбец состоит из единиц, второй — значения первого предиктора, третий — второго и так далее, а строки соответствуют имеющимся наблюдениям.

Функция summary.lm() и оценка получившихся результатов

Теперь рассмотрим пример построения модели линейной регрессии в языке R:

> library(faraway)
> lm1<-lm(Species~Area+Elevation+Nearest+Scruz+Adjacent, data=gala)
> summary(lm1)

Call:
lm(formula = Species ~ Area + Elevation + Nearest + Scruz + Adjacent, 
    data = gala)

Residuals:
     Min       1Q   Median       3Q      Max 
-111.679  -34.898   -7.862   33.460  182.584 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  7.068221  19.154198   0.369 0.715351    
Area        -0.023938   0.022422  -1.068 0.296318    
Elevation    0.319465   0.053663   5.953 3.82e-06 ***
Nearest      0.009144   1.054136   0.009 0.993151    
Scruz       -0.240524   0.215402  -1.117 0.275208    
Adjacent    -0.074805   0.017700  -4.226 0.000297 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 60.98 on 24 degrees of freedom
Multiple R-squared:  0.7658,	Adjusted R-squared:  0.7171 
F-statistic:  15.7 on 5 and 24 DF,  p-value: 6.838e-07

Таблица gala содержит некоторые данные о 30 Галапагосских островах. Мы будем рассматривать модель, где Species — количество разных видов растений на острове линейно зависит от нескольких других переменных.

Рассмотрим вывод функции summary.lm().
Сначала идет строка, которая напоминает, как строилась модель.
Затем идет информация о распределении остатков: минимум, первая квартиль, медиана, третья квартиль, максимум. В этом месте было бы полезно не только посмотреть на некоторые квантили остатков, но и проверить их на нормальность, например тестом Шапиро-Уилка.
Далее — самое интересное — информация о коэффициентах. Здесь потребуется немного теории.
Сначала выпишем следующий результат:

при этом сигма в квадрате с крышкой является несмещенной оценкой для реальной сигмы в квадрате. Здесь b — реальный вектор коэффициентов, а эпсилон с крышкой — вектор остатков, если в качестве коэффициентов взять оценки, полученные методом наименьших квадратов. То есть при предположении, что ошибки распределены нормально, вектор коэффициентов тоже будет распределен нормально вокруг реального значения, а его дисперсию можно несмещенно оценить. Это значит, что можно проверять гипотезу на равенство коэффициентов нулю, а следовательно проверять значимость предикторов, то есть действительно ли величина Xi сильно влияет на качество построенной модели.
Для проверки этой гипотезы нам понадобится следующая статистика, имеющая распределение Стьюдента в том случае, если реальное значение коэффициента bi равно 0:

где
— стандартная ошибка оценки коэффициента, а t(k-n-1) — распределение Стьюдента с k-n-1 степенями свободы.

Теперь все готово для продолжения разбора вывода функции summary.lm().
Итак, далее идут оценки коэффициентов, полученные методом наименьших квадратов, их стандартные ошибки, значения t-статистики и p-значения для нее. Обычно p-значение сравнивается с каким-нибудь достаточно малым заранее выбранным порогом, например 0.05 или 0.01. И если значение p-статистики оказывается меньше порога, то гипотеза отвергается, если же больше, ничего конкретного, к сожалению, сказать нельзя. Напомню, что в данном случае, так как распределение Стьюдента симметричное относительно 0, то p-значение будет равно 1-F(|t|)+F(-|t|), где F — функция распределения Стьюдента с k-n-1 степенями свободы. Также, R любезно обозначает звездочками значимые коэффициенты, для которых p-значение достаточно мало. То есть, те коэффициенты, которые с очень малой вероятностью равны 0. В строке Signif. codes как раз содержится расшифровка звездочек: если их три, то p-значение от 0 до 0.001, если две, то оно от 0.001 до 0.01 и так далее. Если никаких значков нет, то р-значение больше 0.1.

В нашем примере можно с большой уверенностью сказать, что предикторы Elevation и Adjacent действительно с большой вероятностью влияют на величину Species, а вот про остальные предикторы ничего определенного сказать нельзя. Обычно, в таких случаях предикторы убирают по одному и смотрят, насколько изменяются другие показатели модели, например BIC или Adjusted R-squared, который будет разобран далее.

Значение Residual standart error соответствует просто оценке сигмы с крышкой, а степени свободы вычисляются как k-n-1.

А теперь самая важные статистики, на которые в первую очередь стоит смотреть: R-squared и Adjusted R-squared:

где Yi — реальные значения Y в каждом наблюдении, Yi с крышкой — значения, предсказанные моделью, Y с чертой — среднее по всем реальным значениям Yi.

Начнем со статистики R-квадрат или, как ее иногда называют, коэффициента детерминации. Она показывает, насколько условная дисперсия модели отличается от дисперсии реальных значений Y. Если этот коэффициент близок к 1, то условная дисперсия модели достаточно мала и весьма вероятно, что модель неплохо описывает данные. Если же коэффициент R-квадрат сильно меньше, например, меньше 0.5, то, с большой долей уверенности модель не отражает реальное положение вещей.

Однако, у статистики R-квадрат есть один серьезный недостаток: при увеличении числа предикторов эта статистика может только возрастать. Поэтому, может показаться, что модель с большим количеством предикторов лучше, чем модель с меньшим, даже если все новые предикторы никак не влияют на зависимую переменную. Тут можно вспомнить про принцип бритвы Оккама. Следуя ему, по возможности, стоит избавляться от лишних предикторов в модели, поскольку она становится более простой и понятной. Для этих целей была придумана статистика скорректированный R-квадрат. Она представляет собой обычный R-квадрат, но со штрафом за большое количество предикторов. Основная идея: если новые независимые переменные дают большой вклад в качество модели, значение этой статистики растет, если нет — то наоборот уменьшается.

Для примера рассмотрим ту же модель, что и раньше, но теперь вместо пяти предикторов оставим два:

> lm2<-lm(Species~Elevation+Adjacent, data=gala)
> summary(lm2)

Call:
lm(formula = Species ~ Elevation + Adjacent, data = gala)

Residuals:
    Min      1Q  Median      3Q     Max 
-103.41  -34.33  -11.43   22.57  203.65 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.43287   15.02469   0.095 0.924727    
Elevation    0.27657    0.03176   8.707 2.53e-09 ***
Adjacent    -0.06889    0.01549  -4.447 0.000134 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 60.86 on 27 degrees of freedom
Multiple R-squared:  0.7376,	Adjusted R-squared:  0.7181 
F-statistic: 37.94 on 2 and 27 DF,  p-value: 1.434e-08

Как можно увидеть, значение статистики R-квадрат снизилось, однако значение скорректированного R-квадрат даже немного возросло.

Теперь проверим гипотезу о равенстве нулю всех коэффициентов при предикторах. То есть, гипотезу о том, зависит ли вообще величина Y от величин Xi линейно. Для этого можно использовать следующую статистику, которая, если гипотеза о равенстве нулю всех коэффициентов верна, имеет распределение Фишера c n и k-n-1 степенями свободы:

Значение F-статистики и p-значение для нее находятся в последней строке вывода функции summary.lm().

Заключение

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

Синтаксис

• 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