Глава 10. Метод наименьших квадратов
В этой и следующей главах представлена концепция подгонки модели к данным. В данном контексте модель состоит из математического описания взаимосвязи между переменными (например, в виде прямой линии) и описания случайной изменчивости (например, в виде нормального распределения).
Когда говорят, что модель соответствует данным, обычно подразумевают, что она минимизирует ошибки, которые представляют собой расстояния между модельной аппроксимацией и данными. Начнем с одного из самых распространенных способов подгонки модели, минимизирующего сумму квадратов ошибок, — метода наименьших квадратов (МНК).
Сначала вы познакомитесь с моделями, которые работают только с двумя переменными. В следующей главе представлены модели, которые способны описывать большее количество переменных.
Линейная регрессия методом наименьших квадратов
Для примера вернемся к сценарию, описанному в разделе «Вес пингвинов» главы 8 на с. 144. Предположим, вы исследователь в Антарктиде, изучающий местные популяции пингвинов. В процессе сбора данных вы отлавливаете выборку пингвинов, измеряете и взвешиваете их, а затем отпускаете целыми и невредимыми.
Как вскоре становится известно, не так просто заставить пингвина стоять на весах нужное время, чтобы провести точные измерения. Представьте, что для некоторых пингвинов у нас есть такие показатели, как размеры крыльев и клюва, но нет показаний веса. Посмотрим, сможем ли мы использовать другие измерения для получения недостающих данных. Этот процесс называют восстановлением пропущенных данных или импутацией (imputation).
Начнем с изучения взаимосвязи между весом и размерами, используя данные, собранные в период с 2007 по 2010 год исследователями на станции Палмер в Антарктиде. Собранные ими данные находятся в свободном доступе: инструкции по их загрузке приведены в Jupyter-блокноте этой главы.
Для чтения данных можно использовать функцию Pandas read_csv:
penguins = pd.read_csv("penguins_raw.csv").dropna(subset=["Body Mass (g)"])
penguins.shape
(342, 17)
В датасете имеются измерения для 151 пингвина Адели (Adélie penguin). Воспользуемся методом query для выбора строк, содержащих эти данные:
adelie = penguins.query('Species.str.startswith("Adelie")')
len(adelie)
151
Теперь посмотрим, насколько точно можно предсказать вес пингвина Адели, если известна длина его крыльев (flipper length). Сначала извлечем эти данные из объекта DataFrame:
xvar = "Flipper Length (mm)"
yvar = "Body Mass (g)"
flipper_length = adelie[xvar]
body_mass = adelie[yvar]
Вот диаграмма рассеяния, показывающая взаимосвязь между этими величинами:
plt.scatter(flipper_length, body_mass, marker=".", alpha=0.5)
decorate(xlabel="Длина крыла (мм)", ylabel="Масса тела (г)")
Похоже, что они связаны, — оценим количественно силу этой взаимосвязи, вычислив коэффициент корреляции:
np.corrcoef(flipper_length, body_mass)[0, 1]
0.4682016942179394
Корреляция составляет около 0.47, поэтому пингвины с более длинными крыльями, как правило, тяжелее. Знание такого факта полезно, поскольку предполагает, что мы можем более точно определить вес пингвина, если знаем длину его крыльев. Но корреляция сама по себе не говорит нам, как получать такие приблизительные оценки. Для этого нам нужно найти линию наилучшего соответствия.
Существует много способов определить «наилучшую» линию, но для таких данных обычно используется линейная регрессия методом наименьших квадратов, которая дает прямую линию, минимизирующую средний квадрат ошибки (MSE).
В библиотеке SciPy есть функция linregress, которая выполняет подгонку по методу наименьших квадратов. Название функции — это сокращение от линейной регрессии (linear regression), еще одного названия описанной модели. Аргументами функции linregress являются значения x и y в указанном порядке:
from scipy.stats import linregress
result = linregress(flipper_length, body_mass)
result
LinregressResult(slope=32.83168975115009, intercept=-2535.8368022002514,
rvalue=0.46820169421793933, pvalue=1.3432645947790051e-09,
stderr=5.076138407990821, intercept_stderr=964.7984274994059)
Возвращаемое значение — объект LinregressResult, содержащий значения углового коэффициента, или наклона (slope), и свободного члена, или пересечения оси ординат (intercept), подогнанной линии, а также другую информацию, которую мы скоро разберем. Наклон составляет около 32.8, это означает, что на каждый дополнительный миллиметр длины крыла приходятся дополнительные 32.8 г веса тела.
Пересечение оси ординат происходит в точке –2535 г, что может показаться абсурдным, поскольку измеренный вес не может быть отрицательным. Возможно, эта цифра приобретет больший смысл, если использовать значения углового коэффициента и свободного члена для вычисления ординаты подогнанной линии (линии регрессии) для средней длины крыльев:
x = flipper_length.mean()
y = result.intercept + result.slope * x
x, y
(189.95364238410596, 3700.662251655629)
Для пингвина со средней длиной крыльев около 190 мм ожидаемый вес тела составляет около 3700 г.
Следующая функция извлекает параметры из объекта, возвращаемого функцией linregress, и для каждого значения x из последовательности xs находит ординату точки на подогнанной линии:
def predict(result, xs):
ys = result.intercept + result.slope * xs
return ys
Название функции, predict, здесь может показаться странным: в человеческом языке предсказание (prediction) обычно относится к происходящему в будущем, но в контексте регрессии точки на подогнанной линии также называются предсказаниями.
Воспользуемся функцией predict, чтобы найти точки на подогнанной линии для некоторого диапазона размеров плавников:
fit_xs = np.linspace(np.min(flipper_length), np.max(flipper_length))
fit_ys = predict(result, fit_xs)
Вот линия, совмещенная с точечной диаграммой данных:
plt.scatter(flipper_length, body_mass, marker=".", alpha=0.5)
plt.plot(fit_xs, fit_ys, color="C1")
decorate(xlabel="Длина крыла (мм)", ylabel="Масса тела (г)")
Как и ожидалось, линия проходит через центр облака точек данных и отражает тенденцию. Некоторые из полученных прогнозов точны, но множество точек данных находятся далеко от линии. Чтобы понять, насколько хорошо (или плохо) предсказываются данные, можно вычислить ошибку предсказания (prediction error), равную расстоянию по вертикали от каждой точки до линии. Следующая функция вычисляет эти ошибки, которые также называются остатками (residuals):
def compute_residuals(result, xs, ys):
fit_ys = predict(result, xs)
return ys - fit_ys
Ниже приведены значения остатков для массы тела как функции длины крыла:
residuals = compute_residuals(result, flipper_length, body_mass)
В качестве примера посмотрим на результаты для первого пингвина из датасета:
x = flipper_length[0]
y = predict(result, x)
x, y
(181.0, 3406.699042757914)
Длина крыльев выбранного пингвина составляет 181 мм, а предсказанная масса тела — 3407 г. Теперь узнаем, какова его реальная масса:
body_mass[0], residuals[0]
(3750.0, 343.30095724208604)
Фактическая масса этого пингвина составляет 3750 г, а остаток после вычитания предсказанного значения — 343 г.
Среднее значение квадратов остатков — это средний квадрат ошибки (MSE) предсказаний:
mse = np.mean(residuals**2)
mse
163098.85902884745
Само по себе это число не имеет особого значения. Оно станет более осмысленным, когда мы вычислим коэффициент детерминации.
Коэффициент детерминации
Представьте, что вы хотите угадать вес пингвина. Если вы знаете длину его крыльев, то можете использовать метод наименьших квадратов, чтобы предсказать вес, а значение MSE — для количественной оценки средней точности предсказаний.
Но что, если вы не знаете длину крыльев? Какое предположение вы сделаете в этом случае? Оказывается, предсказывать среднее значение — это лучшая стратегия с точки зрения минимизации MSE. Если мы всегда предсказываем среднее значение, то ошибки предсказания — это отклонения (deviations) от среднего значения:
deviations = body_mass - np.mean(body_mass)
А MSE — это средний квадрат отклонения (mean squared deviation):
np.mean(deviations**2)
208890.28989956583
Возможно, вы помните, что средний квадрат отклонения — это дисперсия:
np.var(body_mass)
208890.28989956583
Таким образом, если мы всегда угадываем среднее значение, то в качестве MSE можно рассматривать дисперсию масс, а если используем линию регрессии — то дисперсию остатков. Если вычислить отношение этих величин и вычесть его из 1, то результат покажет, насколько уменьшается MSE при использовании длины крыльев для обоснования наших предположений.
Следующая функция вычисляет это значение, которое формально называется коэффициентом детерминации, но, поскольку оно обозначается как R2, многие называют его «R квадрат».
def coefficient_of_determination(ys, residuals):
return 1 - np.var(residuals) / np.var(ys)
В приведенном примере значение R2 равно примерно 0.22, это означает, что подогнанная линия уменьшает MSE на 22 %:
R2 = coefficient_of_determination(body_mass, residuals)
R2
0.21921282646854912
Оказывается, существует взаимосвязь между коэффициентом детерминации R2 и коэффициентом корреляции r. Как вы могли бы догадаться, глядя на эти условные обозначения, r2 = R2. В этом можно убедиться, вычислив квадратный корень из R2:
r = np.sqrt(R2)
r
0.4682016942179397
И сравнив его с коэффициентом корреляции, который мы считали ранее:
corr = np.corrcoef(flipper_length, body_mass)[0, 1]
corr
0.4682016942179394
Они одинаковы, если не считать небольшой разницы, обусловленной погрешностью операций с плавающей точкой.
Функция linregress также вычисляет это значение и возвращает его в качестве атрибута объекта RegressionResult:
result.rvalue
0.46820169421793933
Коэффициенты детерминации и корреляции содержат главным образом одну и ту же информацию, но интерпретируются по-разному:
• корреляция количественно выражает силу взаимосвязи на шкале от –1 до 1;
• R2 измеряет способность подогнанной линии уменьшить MSE.
Кроме того, величина R2 всегда положительна, поэтому по ней нельзя понять, является корреляция положительной или отрицательной.
Минимизация MSE
Ранее я упоминал, что результатом подгонки методом наименьших квадратов является прямая линия, которая минимизирует средний квадрат ошибки (MSE). Мы не будем это доказывать, но можем проверить, если добавим небольшие случайные значения к свободному члену (intercept) и коэффициенту наклона (slope) и посмотрим, не ухудшает ли это MSE:
intercept = result.intercept + np.random.normal(0, 1)
slope = result.slope + np.random.normal(0, 1)
Чтобы прогнать этот тест, нужно создать объект с атрибутами intercept и slope. Воспользуемся для этого объектом SimpleNamespace, определенным в модуле types:
from types import SimpleNamespace
fake_result = SimpleNamespace(intercept=intercept, slope=slope)
fake_result
namespace(intercept=-2535.738911382989, slope=34.24509022936497)
Передадим этот объект в функцию compute_residuals и используем полученные остатки для вычисления MSE:
fake_residuals = compute_residuals(fake_result, flipper_length, body_mass)
fake_mse = np.mean(fake_residuals**2)
Если сравнить полученный результат с MSE для линии наименьших квадратов, то первый будет всегда хуже:
mse, fake_mse, fake_mse > mse
(163098.85902884745, 235318.11301937344, True)
Достичь минимума MSE — хорошо, но это не единственный критерий «наилучшей» модели. Альтернативный вариант — минимизировать абсолютные значения ошибок. Еще один — минимизировать кратчайшее расстояние от каждой точки до подогнанной линии, которое называется суммарной ошибкой (total error). В одних случаях предполагаемое значение лучше завысить, в других — занизить. В таком случае может понадобиться вычислить функцию потерь (cost function) для каждого остатка и минимизировать общую величину потерь (cost).
Впрочем, метод наименьших квадратов используется гораздо чаще, чем упомянутые альтернативы, в первую очередь потому, что он вычислительно эффективен. Следующая функция демонстрирует это:
def least_squares(xs, ys):
xbar = np.mean(xs)
ybar = np.mean(ys)
xdev = xs — xbar
ydev = ys - ybar
slope = np.sum(xdev * ydev) / np.sum(xdev**2)
intercept = ybar - slope * xbar
return intercept, slope
Для проверки этой функции снова используем длину крыльев и массу тела:
intercept, slope = least_squares(flipper_length, body_mass)
intercept, slope
(-2535.8368022002524, 32.831689751150094)
И можно убедиться, что получены те же результаты, что и при использовании функции linregress:
np.allclose([intercept, slope], [result.intercept, result.slope])
True
Минимизация MSE имела смысл, когда эффективность вычислений превалировала над выбором метода, наиболее подходящего для конкретной задачи. Но это более не актуально, поэтому стоит поразмышлять о том, насколько правильно минимизировать именно квадраты остатков.
Оценка
Параметры slope и intercept — это оценки, полученные на основе некоторой выборки. Подобно другим оценкам, на них сказываются нерепрезентативность выборки, погрешности измерений и изменчивость вследствие случайности выборки. Как правило, влияние нерепрезентативной выборки и ошибок измерений численно оценить трудно. Гораздо проще измерить влияние процесса составления случайной выборки.
Один из способов это сделать — использовать разновидность ресемплинга, называемую бутстрепингом (bootstrapping). Будем считать, что выборка представляет собой всю популяцию, и генерировать из наблюдаемых данных новые выборки с возвращением (replacement). Следующая функция принимает на вход объект DataFrame, использует метод sample для повторной выборки строк и возвращает новый объект DataFrame:
def resample(df):
n = len(df)
return df.sample(n, replace=True)
А функция ниже принимает на вход объект DataFrame, находит линию наилучшего соответствия по методу наименьших квадратов и возвращает наклон найденной подогнанной линии:
def estimate_slope(df):
xs, ys = df["Flipper Length (mm)"], df["Body Mass (g)"]
result = linregress(xs, ys)
return result.slope
Воспользуемся этими функциями для создания множества моделей датасетов и вычисления коэффициента наклона в случае каждого из них:
resampled_slopes = [estimate_slope(resample(adelie)) for i in range(1001)]
В результате получаем выборку из выборочного распределения коэффициента наклона. Вот как она выглядит:
from thinkstats import plot_kde
plot_kde(resampled_slopes)
decorate(xlabel="Коэффициент наклона подогнанной линии (г/мм)",
ylabel="Плотность вероятности")
С помощью функции percentile можно вычислить 90%-ный доверительный интервал (ДИ):
ci90 = np.percentile(resampled_slopes, [5, 95])
print(result.slope, ci90)
32.83168975115009 [25.39604591 40.21054526]
Таким образом, в своем отчете мы могли бы указать, что оценка углового коэффициента составляет 33 г/мм с 90%-ным ДИ [25, 40] г/мм.
Стандартная ошибка оценки — это стандартное отклонение выборочного распределения:
stderr = np.std(resampled_slopes)
stderr
4.570238986584832
Объект RegressionResult, возвращаемый функцией linregress, содержит аппроксимацию стандартной ошибки, исходя из некоторых предположений о форме распределения:
result.stderr
5.076138407990821
Стандартная ошибка, которую мы вычислили путем ресемплинга, была немного меньше, но на практике разница, вероятно, роли не играет.
Визуализация неопределенности
Каждый раз, когда мы проводим ресемплинг датасета, мы получаем новую подогнанную линию. Чтобы увидеть, насколько велик разброс линий, можно, например, пройтись по ним циклом и отобразить их все. Функция ниже принимает на вход результат ресемплинга в виде объекта DataFrame, находит линию наилучшего соответствия по методу наименьших квадратов и для последовательности xs генерирует предсказанные значения:
def fit_line(df, fit_xs):
xs, ys = df["Flipper Length (mm)"], df["Body Mass (g)"]
result = linregress(xs, ys)
fit_ys = predict(result, fit_xs)
return fit_ys
Вот последовательность xs, которую мы будем использовать:
xs = adelie["Flipper Length (mm)"]
fit_xs = np.linspace(np.min(xs), np.max(xs))
А вот как выглядят подогнанные линии, совмещенные с диаграммой рассеяния исходных данных:
plt.scatter(flipper_length, body_mass, marker=".", alpha=0.5)
for i in range(101):
fit_ys = fit_line(resample(adelie), fit_xs)
plt.plot(fit_xs, fit_ys, color="C1", alpha=0.05)
decorate(xlabel="Длина крыла (мм)", ylabel="Масса тела (г)")
В середине диаграммы линии расположены близко друг к другу, а в крайних точках — расходятся.
Другой способ представить изменчивость подогнанных линий — построить 90%-ный доверительный интервал для каждого предсказанного значения. Начнем с того, что соберем подогнанные линии в списке массивов:
fitted_ys = [fit_line(resample(adelie), fit_xs) for i in range(1001)]
Этот список массивов можно считать двумерным массивом, где каждая строка соответствует подогнанной линии, а столбец — значению из xs.
Чтобы найти 5-й, 50-й и 95-й процентили ys, соответствующие каждому из значений в xs, можно вызвать функцию percentile с аргументом axis=0:
low, median, high = np.percentile(fitted_ys, [5, 50, 95], axis=0)
Теперь воспользуемся Matplotlib-функцией fill_between для заполнения области между 5-м и 95-м процентилями, представляющей 90 % ДИ, отображения медианных значений всех столбцов и отрисовки точечной диаграммы исходных данных:
plt.scatter(flipper_length, body_mass, marker=".", alpha=0.5)
plt.fill_between(fit_xs, low, high, color="C1", lw=0, alpha=0.2)
plt.plot(fit_xs, median, color="C1")
decorate(xlabel="Длина крыла (мм)", ylabel="Масса тела (г)")
Это мой любимый способ визуализации изменчивости подогнанной линии, обусловленной случайностью выборки.
Преобразование переменных
Перед тем как приступать к подгонке линии регрессии данных, иногда полезно преобразовать одну или обе переменные — например, возведя значения в квадрат либо рассчитав их квадратный корень или логарифм. В качестве примера рассмотрим данные о росте и весе из Системы наблюдения за поведенческими факторами риска (BRFSS), описанные в разделе «Логнормальное распределение» главы 5 на с. 94.
Сначала загрузим данные BRFSS:
from thinkstats import read_brfss
brfss = read_brfss()
Далее отберем строки с корректными данными и выберем столбцы, содержащие значения роста и веса:
valid = brfss.dropna(subset=["htm3", "wtkg2"])
heights, weights = valid["htm3"], valid["wtkg2"]
С помощью функции linregress можно по методу наименьших квадратов получить значения коэффициента наклона и точки пересечения ординаты линией регрессии:
result_brfss = linregress(heights, weights)
result_brfss.intercept, result_brfss.slope
(-82.65926054409877, 0.957074585033226)
Коэффициент наклона составляет около 0.96, что можно интерпретировать как прибавление в среднем почти 1 кг веса на каждый добавочный 1 см роста. Снова воспользуемся функцией predict для предсказания значений зависимой переменной в диапазоне значений независимой xs:
fit_xs = np.linspace(heights.min(), heights.max())
fit_ys = predict(result_brfss, fit_xs)
Прежде чем строить диаграмму рассеяния данных, полезно применить джиттеринг к значениям роста и веса:
from thinkstats import jitter
jittered_heights = jitter(heights, 2)
jittered_weights = jitter(weights, 1.5)
Также используем среднее значение и стандартное отклонение величины роста, чтобы подобрать граничные значения на оси x:
m, s = heights.mean(), heights.std()
xlim = m - 4 * s, m + 4 * s
ylim = 0, 200
Вот диаграмма рассеяния данных после добавления шума вместе с подогнанной линией:
plt.scatter(jittered_heights, jittered_weights, alpha=0.01, s=0.1)
plt.plot(fit_xs, fit_ys, color="C1")
decorate(xlabel="Рост (см)", ylabel="Вес(кг)", xlim=xlim, ylim=ylim)
Линия регрессии на диаграмме проходит несколько выше области наибольшего скопления точек данных. Это объясняется тем, что вес не подчиняется нормальному распределению. Как вы узнали из раздела «Логнормальное распределение» главы 5 на с. 94, вес взрослых людей обычно описывается логнормальным распределением, которое скошено в сторону бо́льших значений. Именно этот перекос смещает подогнанную линию вверх.
Еще один момент, вызывающий беспокойство, — распределение остатков, которое выглядит так:
residuals = compute_residuals(result_brfss, heights, weights)
from thinkstats import make_pmf
pmf_kde = make_pmf(residuals, -60, 120)
pmf_kde.plot()
decorate(xlabel="Остаток регрессии (кг)", ylabel="Плотность вероятности")
Распределение остатков скошено вправо. Сам по себе этот факт может не представлять сложности, но он подразумевает, что методом наименьших квадратов не удалось должным образом описать взаимосвязь между этими переменными.
Если вес описывается логнормальным распределением, то его логарифм должен соответствовать нормальному распределению. Поэтому посмотрим, как изменится результат, если линия регрессии будет выражать логарифм веса как функцию показаний роста:
log_weights = np.log10(weights)
result_brfss2 = linregress(heights, log_weights)
result_brfss2.intercept, result_brfss2.slope
(0.9930804163932876, 0.005281454169417777)
Поскольку мы преобразовали одну из переменных, коэффициент наклона и свободный член теперь сложнее интерпретировать. Но при помощи predict можно найти точки на линии регрессии:
fit_xs = np.linspace(heights.min(), heights.max())
fit_ys = predict(result_brfss2, fit_xs)
А затем изобразить ее на одном графике с диаграммой рассеяния преобразованных данных:
jittered_log_weights = jitter(log_weights, 1.5)
plt.scatter(jittered_heights, jittered_log_weights, alpha=0.01, s=0.1)
plt.plot(fit_xs, fit_ys, color="C1")
decorate(xlabel="Рост (см)", ylabel="Вес (log10 кг)", xlim=xlim)
Теперь линия регрессии проходит через область с наибольшей плотностью точек, а сами точки примерно равноудалены по обе стороны от линии. Поэтому распределение остатков приблизительно симметрично:
residuals = compute_residuals(result_brfss2, heights, log_weights)
pmf_kde = make_pmf(residuals, -0.6, 0.6)
pmf_kde.plot()
decorate(xlabel="Остаток регрессии (кг)", ylabel="Плотность вероятности")
Внешний вид точечной диаграммы и распределение остатков свидетельствуют, что линия регрессии хорошо описывает связь роста и логарифма веса.
Если сравнить значения r обеих регрессий, то окажется, что коэффициент корреляции роста с логарифмом веса немного больше:
result_brfss.rvalue, result_brfss2.rvalue
(0.5087364789734582, 0.5317282605983435)
Это означает, что значение R2 тоже слегка больше:
result_brfss.rvalue**2, result_brfss2.rvalue**2
(0.2588128050383119, 0.28273494311893993)
При расчете предполагаемых значений веса с помощью показателей роста результаты несколько улучшатся, если работать с логарифмом веса.
Однако преобразование данных усложняет интерпретацию параметров модели — улучшить представление результатов может помочь обратное преобразование. Например, обратным логарифму по основанию 10 является возведение в степень по основанию 10. Вот как выглядит линия регрессии после обратного преобразования вместе с исходными данными:
plt.scatter(jittered_heights, jittered_weights, alpha=0.01, s=0.1)
plt.plot(fit_xs, 10**fit_ys, color="C1")
decorate(xlabel="Рост (см)", ylabel="Вес(кг)", xlim=xlim, ylim=ylim)
Линия регрессии, выглядящая как прямая в отношении к преобразованным данным, становится изогнутой применительно к исходным данным.
Глоссарий
Модель (model)
Применительно к задаче регрессии модель — это математическое описание взаимосвязи между переменными, например прямая линия, а также описание случайной изменчивости, например нормальное распределение.
Импутация (imputation)
Процесс оценки и заполнения пропущенных значений в датасете.
Линия наилучшего соответствия (line of best fit)
Прямая линия (или кривая), которая наилучшим образом описывает взаимосвязь между переменными в соответствии с некоторым определением «наилучшего».
Линейная регрессия (linear regression)
Метод поиска прямой линии наилучшего соответствия.
Предсказание (prediction)
Точка на линии наилучшего соответствия. В контексте регрессии необязательно относится к будущему.
Остаток (residual)
Разница между наблюдаемым значением и значением, предсказанным с помощью линии наилучшего соответствия.
Линейная регрессия методом наименьших квадратов (linear least squares fit)
Поиск линии, которая минимизирует сумму квадратов остатков.
Коэффициент детерминации (coefficient of determination)
Статистический показатель, обозначаемый как R2 и часто называемый «R квадрат», который измеряет, насколько хорошо модель соответствует данным.
Бутстрепинг (bootstrap resampling)
Метод ресемплинга (повторной выборки), в котором на основе исходной выборки, рассматриваемой как генеральная совокупность, осуществляется новая выборка того же размера, что и исходная, с возвращением.
Упражнения
Упражнение 10.1
В этой главе мы методом наименьших квадратов находили зависимость веса пингвинов от длины их крыльев. В использованном датасете есть еще два измерения, которые также можно принять во внимание, — длина (culmen length) и глубина (culmen depth) верхней части клюва пингвина.
С помощью метода наименьших квадратов найдите линию регрессии, выражающую зависимость веса от длины верхней части клюва. Постройте диаграмму рассеяния для этих переменных и нанесите на нее подогнанную линию.
Ориентируясь на значение атрибута rvalue объекта RegressionResult, ответьте на вопросы: какова величина корреляции этих переменных? каков коэффициент детерминации? что лучше предсказывает вес — длина верхней части клюва или длина крыла?
Упражнение 10.2
В этой главе мы использовали ресемплинг, чтобы аппроксимировать выборочное распределение коэффициента наклона линии регрессии. Точно так же можно аппроксимировать выборочное распределение для свободного члена (пересечения оси ординат).
1. Создайте функцию estimate_intercept, которая принимает в качестве аргумента объект DataFrame — результат ресемплинга, по методу наименьших квадратов находит линейную аппроксимацию веса как функцию длины крыла и возвращает величину свободного члена (intercept).
2. Пройдите в цикле по большому множеству повторных выборок датафрейма adelie, вызывая данную функцию, и соберите все значения свободного члена.
3. Используйте функцию plot_kde для построения графика выборочного распределения свободного члена.
4. Вычислите стандартную ошибку и 90%-ный доверительный интервал.
5. Убедитесь, что стандартная ошибка, которую вы получаете при повторной выборке, согласуется со значением атрибута intercept_stderr объекта RegressionResult: она может быть несколько меньше.
Упражнение 10.3
Индекс массы тела, ИМТ (Body Mass Index, BMI), человека равен весу в килограммах, деленному на квадрат роста в метрах. В датасете BRFSS мы можем рассчитать ИМТ после перевода роста из сантиметров в метры:
heights_m = heights / 100
bmis = weights / heights_m**2
В этом определении рост возводится во вторую, а не какую-либо другую, степень, поскольку на заре истории статистики было обнаружено, что средний вес увеличивается ориентировочно пропорционально квадрату роста.
Чтобы проверить, так ли это, воспользуемся данными из датасета BRFSS, методом наименьших квадратов и математическими приемами. Предположим, что вес пропорционален росту, возведенному в степень с неизвестным показателем a. В этом случае можно записать:
w = bha,
где w — вес, h — рост, а b — неизвестная константа.
Логарифмируя обе части равенства, получаем:
log w = log b + a log h.
Поэтому, если методом наименьших квадратов вычислить линейную регрессию для логарифма веса как функции логарифма роста, то коэффициент наклона линии регрессии будет давать оценку неизвестного показателя степени a.
Вычислите логарифмы роста и веса. Можно использовать любое основание логарифма, если только оно одинаково для обоих преобразований. Найдите линию регрессии методом наименьших квадратов. Насколько коэффициент наклона близок к 2?
