Глава 11. Множественная регрессия

Линейный метод наименьших квадратов, описанный в предыдущей главе, является примером регрессии. В более общем плане регрессия представляет собой задачу моделирования взаимосвязи между одним набором переменных, называемых зависимыми переменными (dependent variable) или откликом (response variable), и другим набором переменных, именуемыми объясняющими переменными (explanatory variable) или независимыми переменными (independent variable).

В примерах, приведенных в предыдущей главе, есть только одна зависимая переменная и одна объясняющая переменная. Такой случай называется простой регрессией (simple regression). В этой главе вы познакомитесь с множественной регрессией (multiple regression) с несколькими объясняющими переменными, но по-прежнему только одной переменной-откликом. Случай, когда отклика более одного, представляет собой многомерную регрессию (multivariate regression). Ее мы не будем рассматривать в этой книге.

Модуль StatsModels

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

В первом примере этой главы продолжим изучение датасета о пингвинах. При загрузке данных зададим столбцам переменных имена без пробелов — это упростит их использование со StatsModels. Для этой цели пригодится следующий словарь:

columns = {

    "Body Mass (g)": "mass",

    "Flipper Length (mm)": "flipper_length",

    "Culmen Length (mm)": "culmen_length",

    "Culmen Depth (mm)": "culmen_depth",

}

Теперь загрузим данные, удалим строки с отсутствующими значениями массы тела и переименуем столбцы:

penguins = (

    pd.read_csv("penguins_raw.csv")

    .dropna(subset=["Body Mass (g)"])

    .rename(columns=columns)

)

penguins.shape

(342, 17)

Наш датасет содержит данные по трем видам пингвинов. Мы будем работать только с пингвинами Адели:

adelie = penguins.query('Species.str.startswith("Adelie")').copy()

len(adelie)

151

В предыдущей главе мы находили линию зависимости веса пингвина от длины его крыльев методом наименьших квадратов:

flipper_length = adelie["flipper_length"]

body_mass = adelie["mass"]

Напомню, как мы это делали с помощью функции linregress:

from scipy.stats import linregress

 

result_linregress = linregress(flipper_length, body_mass)

result_linregress.intercept, result_linregress.slope

(-2535.8368022002514, 32.83168975115009)

Модуль StatsModels имеет два типа программных интерфейсов (API). Мы будем использовать «формульный» API, в котором зависимую и объясняющие переменные можно задать с помощью языка формул пакета Patsy. Следующая формульная запись означает, что зависимая переменная mass является линейной функцией одной независимой переменной flipper_length:

formula = "mass ~ flipper_length"

Теперь эту формулу вместе с данными можно передать в функцию ols модуля StatsModels:

import statsmodels.formula.api as smf

 

model = smf.ols(formula, data=adelie)

type(model)

statsmodels.regression.linear_model.OLS

Название ols расшифровывается как «обычный метод наименьших квадратов» (ordinary least squares). Слово «обычный» указывает на то, что эта функция реализует метод наименьших квадратов в условиях наиболее распространенного, или «обычного», набора допущений.

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

Метод fit подгоняет модель под данные и возвращает объект типа Regres­sionResults, содержащий результаты:

result_ols = model.fit()

Объект RegressionResults содержит много информации, поэтому модуль thinkstats предоставляет функцию, отображающую только ту информацию, которая нам сейчас понадобится:

from thinkstats import display_summary

 

display_summary(result_ols)

	coef	std err	t	P>| t |	[0.025	0.975]
Intercept	–2535.8368	–964.798	–2.628	0.009	–4442.291	–629.382
flipper_length	32.8317	5.076	6.468	0.000	22.801	42.862

R-squared: 0.2192

В первом столбце указаны коэффициенты модели — свободный член (intercept) и наклон (slope). Можно убедиться, что они совпадают с коэффициентами, полученными из функции linregress:

result_linregress.intercept, result_linregress.slope

(-2535.8368022002514, 32.83168975115009)

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

result_linregress.intercept_stderr, result_linregress.stderr,

(964.7984274994059, 5.076138407990821)

В следующем столбце представлена t-статистика, которая используется для вычисления p-значений. Мы можем пропустить его, поскольку p-значения приведены в следующем столбце с заголовком P>| t |. p-значение для коэффициента при переменной flipper_length округлено до 0, но можно отобразить его так:

result_ols.pvalues["flipper_length"]

1.3432645947789321е-09

А затем убедиться, что linregress выдала тот же результат:

result_linregress.pvalue

1.3432645947790051е-09

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

В последних двух столбцах, обозначенных как [0.025 и 0.975], указан 95%-ный доверительный интервал для свободного члена и наклона. Таким образом, 95 %-ный ДИ для коэффициента наклона составляет [22.8, 42.9].

В последней строке выведено значение R2 модели — около 0.22. Оно подразумевает, что, если вместо среднего значения веса предсказывать вес с использованием длины крыла, величину MSE можно уменьшить примерно на 22 %.

Получаемое при расчете простой корреляции значение R2 равно квадрату коэффициента корреляции r. Поэтому можно сравнить rsquared, вычисленный с помощью функции ols, с квадратом rvalue, вычисленным функцией linregress:

result_ols.rsquared, result_linregress.rvalue**2

(0.21921282646854878, 0.21921282646854875)

Если пренебречь небольшой разницей из-за погрешности операций с плавающей точкой, они одинаковы.

Прежде чем переходить к множественной регрессии, рассчитаем еще одну простую регрессию, взяв culmen_length (длина верхнего края клюва пингвина) в качестве объясняющей переменной:

formula = "mass ~ culmen_length"

result = smf.ols(formula, data=adelie).fit()

display_summary(result)

	coef	std err	t	P>| t |	[0.025	0.975]
Intercept	34.8830	458.439	0.076	0.939	–870.998	940.764
culmen_length	94.4998	11.790	8.015	0.000	71.202	117.798

R-squared: 0.3013

И опять полученное p-значение коэффициента наклона очень мало — случайно наблюдать такую величину наклона в отсутствие связи между массой и длиной клюва крайне маловероятно. Как вы могли заметить, p-значение, связанное со свободным членом, велико. Но этот факт нас не должен волновать, поскольку мы не проверяем гипотезу о равенстве свободного члена нулю. У данной модели свободный член оказался близким к нулю только по стечению обстоятельств — на проблему с моделью он не указывает.

Коэффициент R2 у этой модели составляет около 0.30. Поэтому если в качестве объясняющей переменной вместо длины крыла использовать длину верхнего края клюва, величина MSE сократится немного больше (значение R2 для длины крыла равно 0.22). Теперь посмотрим, что получится, если их объединить.

Знакомство с множественной регрессией

Вот формула Patsy для модели множественной регрессии, в которой масса является линейной комбинацией длины плавника и длины верхнего края клюва:

formula = "mass ~ flipper_length + culmen_length"

А вот результат подгонки этой модели к данным:

result = smf.ols(formula, data=adelie).fit()

display_summary(result)

	coef	std err	t	P>| t |	[0.025	0.975]
Intercept	–3573.0817	866.739	–4.122	0.000	–5285.864	–1860.299
flipper_length	22.7024	4.742	4.787	0.000	13.331	32.074
culmen_length	76.3402	11.644	6.556	0.000	53.331	99.350

R-squared: 0.3949

У этой модели три коэффициента — свободный член и два коэффициента наклона. Наклон для длины крыла равен 22.7. Итак, можно ожидать, что пингвин с крылом, которое длиннее на один миллиметр, будет тяжелее на 22.7 г при условии, что длина верхнего края клюва такая же. Точно так же мы ожидаем, что пингвин с клювом, который длиннее на один миллиметр, будет тяжелее на 76.3 г при условии, что длина крыльев одинакова.

У обоих коэффициентов наклона p-значения малы, это означает, что вклад обеих объясняющих переменных вряд ли случаен.

А коэффициент R2 равен 0.39 — это выше, чем у модели только с длиной верхней части клюва (0.30) и модели с одной длиной крыла (0.22). Таким образом, предсказания с использованием обеих объясняющих переменных точнее, чем предсказания с помощью только одной из них.

Но не настолько, как можно было ожидать. Если длина крыла уменьшает показатель MSE на 22 %, а длина клюва — на 30 %, то почему бы им обеим вместе не уменьшить его в общей сложности на 52 %? Причина в том, что объясняющие переменные коррелируют друг с другом:

from thinkstats import corrcoef

 

corrcoef(adelie, "flipper_length", "culmen_length")

0.32578471516515944

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

Мы увидим ту же картину, если добавим глубину верхнего края клюва (culmen depth) в качестве третьей объясняющей переменной:

formula = "mass ~ flipper_length + culmen_length + culmen_depth"

result = smf.ols(formula, data=adelie).fit()

display_summary(result)

	coef	std err	t	P > | t |	[0.025	0.975]
Intercept	–4341.3019	795.117	–5.460	0.000	–5912.639	–2769.964
flipper_length	17.4215	4.385	3.973	0.000	8.756	26.087
culmen_length	55.3676	11.133	4.973	0.000	33.366	77.369
culmen_depth	140.8946	24.216	5.818	0.000	93.037	188.752

R-squared: 0.5082

У этой модели четыре коэффициента. Все p-значения малы, а значит, вклад какой-либо независимой переменной вряд ли можно объяснить случайностью. А коэффициент R2 равен примерно 0.51, что несколько лучше, чем у предыдущей модели с двумя независимыми переменными (0.39), и заметно лучше, чем у любой из моделей с одной переменной (0.22 и 0.30).

Но и в этот раз его увеличение меньше, чем мы могли бы надеяться, потому что величина глубины клюва коррелирует с двумя другими измерениями:

[

    corrcoef(adelie, "culmen_depth", "flipper_length"),

    corrcoef(adelie, "culmen_depth", "culmen_length"),

]

[0.30762017939668534, 0.39149169183587634]

Этот пример показывает, как обычно используется множественная регрессия, объединяющая несколько независимых переменных для получения более точных прогнозов. Еще одно ее распространенное применение — количественная оценка вклада одного набора переменных при одновременном контроле вклада другого набора.

Контрольная переменная

Из раздела «Сравнение ИФР» главы 4 на с. 67 вы узнали, что первенцы в среднем легче, чем другие новорожденные. А в разделе «Проверка значимости корреляции» главы 9 на с. 170 убедились, что вес новорожденного коррелирует с возрастом матери: чем старше мать, тем тяжелее в среднем ее дети.

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

Инструкции по загрузке данных NSFG приведены в Jupyter-блокноте этой главы.

Функция get_nsfg_groups поможет со считыванием данных, выбором случаев рождения живых детей и разбивкой последних на первенцев и остальных детей:

from nsfg import get_nsfg_groups

 

live, firsts, others = get_nsfg_groups()

Используем метод dropna для отбора строк с корректными значениями веса при рождении, порядка рождения и возраста матери:

valid = live.dropna(subset=["agepreg", "birthord", "totalwgt_lb"]).copy()

Теперь задействуем функционал модуля StatsModels, чтобы подтвердить, что вес новорожденного коррелирует с возрастом матери, и оценить коэффициент наклона в предположении, что эта зависимость линейна:

formula = "totalwgt_lb ~ agepreg"

result_age = smf.ols(formula, data=valid).fit()

display_summary(result_age)

	coef	std err	t	P>| t |	[0.025	0.975]
Intercept	6.8304	0.068	100.470	0.000	6.697	6.964
agepreg	0.0175	0.003	6.559	0.000	0.012	0.023

R-squared: 0.004738

Коэффициент наклона мал — всего 0.0175 фунта на год. Поэтому, если две матери различаются в возрасте на десять лет, можно ожидать, что вес их детей будет отличаться на 0.175 фунта. Но p-значение тоже мало, так что этот наклон, каким бы малым он ни был, был бы маловероятен, если бы на самом деле не было никакой взаимосвязи.

Значение R2 также мало, а значит, возраст матери не очень полезен в качестве предсказывающей переменной. Если мы узнаем возраст матери, наша способность предсказать вес ребенка практически не улучшится.

Такое сочетание малых величин p-значения и коэффициента R2 часто сбивает с толку, поскольку кажется противоречивым: если связь статистически значима, то кажется, что она должна обладать предсказательной силой. Но данный пример показывает, что противоречия нет: взаимосвязь может быть статистически значимой, но малополезной для целей предсказания. Если мы визуализируем результаты, то увидим почему. Сначала выберем соответствующие столбцы:

totalwgt = valid["totalwgt_lb"]

agepreg = valid["agepreg"]

Чтобы посчитать координаты линии регрессии, можно было бы извлечь величины свободного члена и коэффициента наклона из result_age, но в этом нет необходимости. Объект RegressionResults предоставляет метод predict, который можно использовать взамен. Сначала сгенерируем массив значений для переменной agepreg:

agepreg_range = np.linspace(agepreg.min(), agepreg.max())

Чтобы воспользоваться методом predict, нужно поместить значения объясняющих переменных в объект DataFrame:

df = pd.DataFrame({"agepreg": agepreg_range})

Столбцы в датафрейме должны иметь те же имена, что и объясняющие переменные. Теперь его можно передать методу predict:

fit_ys = result_age.predict(df)

Метод возвращает объект Series с предсказанными значениями. Вот как они выглядят вместе с диаграммой рассеяния данных:

plt.scatter(agepreg, totalwgt, marker=".", alpha=0.1, s=5)

plt.plot(agepreg_range, fit_ys, color="C1", label="линейная модель")

 

decorate(xlabel="Возраст матери", ylabel="Вес при рождении (фунты)")

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

Теперь с помощью модуля StatsModels подтвердим, что первенцы при рождении весят меньше, чем последующие дети. Для этого сначала создадим объект Series с логическими значениями True для первых детей и False для всех остальных и добавим в него новый столбец с именем is_first:

valid["is_first"] = valid["birthord"] == 1

Далее идет формула для модели с весом новорожденного в качестве зависимой переменной и is_first в качестве объясняющей переменной. В формульном языке Patsy буква C с круглыми скобками вокруг имени переменной указывает на то, что переменная является категориальной (categorical), то есть представляющей такие категории, как «первый ребенок», а не измерения, такие как вес новорожденного:

formula = "totalwgt_lb ~ C(is_first)"

Теперь выполним подгонку (fit) модели и выведем результаты как обычно:

result_first = smf.ols(formula, data=valid).fit()

display_summary(result_first)

	coef	std err	t	P > | t |	[0.025	0.975]
Intercept	7.3259	0.021	356.007	0.000	7.286	7.366
C(is_first)[T.True]	–0.1248	0.030	–4.212	0.000	–0.183	–0.067

R-squared: 0.00196

В таблице результатов заголовок C(is_first)[T.True] указывает на то, что is_first является категориальной переменной, а коэффициент связан со значением True. Буква Т перед True расшифровывается как «воздействие» (treatment): на языке контролируемого эксперимента первые дети считаются экспериментальной группой (treatment group), а остальные дети — контрольной группой (reference group). Такое соотнесение произвольно: мы могли бы считать, что первенцы являются контрольной группой, а остальные дети — экспериментальной группой. Но для интерпретации результатов нам нужно знать, что есть что.

Свободный член равен примерно 7.3, это означает, что средний вес контрольной группы составляет 7.3 фунта. Коэффициент переменной is_first равен –0.12, что означает, что средний вес детей экспериментальной группы, первенцев, на 0.12 фунта меньше. Мы можем проверить оба этих результата, вычислив их напрямую:

others["totalwgt_lb"].mean()

7.325855614973262

diff_weight = firsts["totalwgt_lb"].mean() - others["totalwgt_lb"].mean()

diff_weight

-0.12476118453549034

В дополнение к этим коэффициентам StatsModels вычисляет p-значения, доверительные интервалы и коэффициент R2. Величина p-значения, связанного с первенцами, мала, а значит, разница между двумя группами статистически значима. А малость R2 в свою очередь показывает, что при попытке угадать вес ребенка знание того, первый ли это ребенок, не очень-то помогает.

Теперь посмотрим, насколько правдоподобно, что разница в весе новорожденных обусловлена разницей в возрасте матерей. В среднем матери первенцев примерно на 3.6 года моложе матерей других детей:

diff_age = firsts["agepreg"].mean() - others["agepreg"].mean()

diff_age

-3.5864347661500275

А угол наклона линии регрессии для веса новорожденного как функции возраста матери составляет 0.0175 фунта на год:

slope = result_age.params["agepreg"]

slope

0.017453851471802638

Если мы умножим величину наклона на разницу в возрасте, то получим ожидаемое различие в весе при рождении первенцев и других детей, обусловленное возрастом матери:

slope * diff_age

-0.0625970997216918

Результат равен 0.063 фунта, что составляет примерно половину наблюдаемой разницы. Таким образом, кажется, что наблюдаемую разницу в весе новорожденных можно только частично объяснить различием в возрасте матерей.

Используя множественную регрессию, можно одновременно оценить коэффициенты для возраста матери и категориальной переменной «первый ребенок»:

formula = "totalwgt_lb ~ agepreg + C(is_first)"

result = smf.ols(formula, data=valid).fit()

display_summary(result)

	coef	std err	t	P > | t |	[0.025	0.975]
Intercept	6.9142	0.078	89.073	0.000	6.762	7.066
C(is_first)[T.True]	–0.0698	0.031	–2.236	0.025	–0.131	–0.009
agepreg	0.0154	0.003	5.499	0.000	0.010	0.021

R-squared: 0.005289

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

А p-значение равно 0.025, что по-прежнему считается статистически значимым, но находится в пограничном диапазоне, где нельзя исключить вероятность того, что разница такого размера могла возникнуть случайно.

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

Нелинейная зависимость

Чтобы проверить, может ли вклад переменной agepreg быть нелинейным, ­добавим в датасет новый столбец со значениями agepreg, возведенными в квадрат:

valid["agepreg2"] = valid["agepreg"] ** 2

Теперь мы можем определить модель, включающую в себя и линейную, и квад­ратичную зависимости:

formula = "totalwgt_lb ~ agepreg + agepreg2"

Выполним подгонку модели как обычно:

result_age2 = smf.ols(formula, data=valid).fit()

display_summary(result_age2)

	coef	std err	t	P > | t |	[0.025	0.975]
Intercept	5.5720	0.275	20.226	0.000	5.032	6.112
agepreg	0.1186	0.022	5.485	0.000	0.076	0.161
agepreg2	–0.0019	0.000	–4.714	0.000	–0.003	–0.001

R-squared: 0.00718

Связанное с членом второй степени, agepreg2, p-значение очень мало, что свидетельствует о том, что его вклад в информацию о весе новорожденного больше, чем могла бы давать чистая случайность. А коэффициент R2 этой модели, равный 0.0072, выше, чем у линейной модели (0.0047).

Оценивая коэффициенты для agepreg и agepreg2, мы фактически подгоняем параболу к данным. Чтобы в этом убедиться, можно воспользоваться объектом RegressionResults для генерации прогнозов для диапазона значений возраста матери.

Сначала создадим временный объект DataFrame, столбцы agepreg и agepreg2 которого заполнены с использованием массива диапазона возраста agepreg_range:

df = pd.DataFrame({"agepreg": agepreg_range})

df["agepreg2"] = df["agepreg"] ** 2

Теперь вызовем метод predict, передав ему этот датафрейм в качестве аргумента и получая на выходе объект Series с предсказаниями:

fit_ys = result_age2.predict(df)

Вот как выглядит подогнанная парабола вместе с диаграммой рассеяния данных:

plt.scatter(agepreg, totalwgt, marker=".", alpha=0.1, s=5)

plt.plot(agepreg_range, fit_ys, color="C1", label="квадратичная модель")

 

decorate(xlabel="Возраст матери", ylabel="Вес при рождении (фунты)")

Хотя кривизна едва различима, ее наличие говорит о том, что вес новорожденного ниже у самых юных и самых возрастных матерей и выше у рожениц из середины возрастного диапазона.

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

formula = "totalwgt_lb ~ agepreg + agepreg2 + C(is_first)"

result = smf.ols(formula, data=valid).fit()

display_summary(result)

	coef	std err	t	P > | t |	[0.025	0.975]
Intercept	5.6923	0.286	19.937	0.000	5.133	6.252
C(is_first)[T.True]	–0.0504	0.031	–1.602	0.109	–0.112	0.011
agepreg	0.1124	0.022	5.113	0.000	0.069	0.155
agepreg2	–0.0018	0.000	–4.447	0.000	–0.003	–0.001

R-squared: 0.007462

При повышении эффективности контроля возраста матери оценочная разница между первыми и остальными детьми стала равна 0.0504 фунта — меньше, чем при использовании только линейной модели (0.0698 фунта). А p-значение, связанное с is_first, стало равно 0.109, то есть остаток разницы между этими группами может быть обусловлен случайностью.

Мы приходим к выводу, что разница в весе новорожденных объясняется — по крайней мере частично, а возможно и полностью — различием в возрасте матери.

Логистическая регрессия

Линейная регрессия по своей сути представляет модель, в которой ожидаемое значение зависимой переменной равно взвешенной сумме объясняющих переменных и свободного члена. Такая модель уместна, когда зависимая переменная является непрерывной величиной, такой как вес новорожденного или масса пингвина, но не когда она представляет собой дискретную величину, такую как натуральное число или категория.

Для такого типа зависимых переменных подходят обобщенные линейные модели, ОЛМ (generalized linear model, GLM). Например:

• Если зависимая переменная — натуральное число, можно использовать регрессию Пуассона.

• Если зависимая переменная категориальная и категорий только две, можно применить логистическую регрессию.

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

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

В книге мы не будем разбирать все эти виды — только логистическую регрессию, которая используется наиболее широко. Для примера снова воспользуемся датасетом о пингвинах и посмотрим, можно ли на основе веса и размеров пингвина определить его пол.

В модуле StatsModels реализована функция, которая считает логистическую регрессию. Она называется logit, потому что математическая функция с таким названием фигурирует в определении логистической регрессии. Прежде чем мы сможем использовать функцию logit, необходимо преобразовать зависимую переменную так, чтобы ее значения были равны 0 и 1:

adelie["y"] = (adelie["Sex"] == "MALE").astype(int)

Начнем с простой модели с y в качестве зависимой переменной и mass в качестве объясняющей. Определим и подгоним ее под данные (аргумент disp=False подавляет сообщения о процессе подгонки):

model = smf.logit("y ~ mass", data=adelie)

result = model.fit(disp=False)

Вот что получается:

display_summary(result)

	coef	std err	z	P > | z |	[0.025	0.975]
Intercept	–25.9871	4.221	–6.156	0.000	–34.261	–17.713
mass	0.0070	0.001	6.138	0.000	0.005	0.009

Pseudo R-squared: 0.5264

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

Коэффициент при переменной mass положителен, это означает, что более тяжелые пингвины, вероятнее всего, самцы. В остальном эти коэффициенты нелегко интерпретировать. Мы сможем лучше понять эту модель, если построим графики предсказанных значений. Для этого создадим объект DataFrame с диапазоном значений переменной mass и воспользуемся методом predict, чтобы получить объект Series с предсказаниями модели:

mass = adelie["mass"]

mass_range = np.linspace(mass.min(), mass.max())

df = pd.DataFrame({"mass": mass_range})

fit_ys = result.predict(df)

Каждое предсказанное значение как функция веса пингвина — вероятность того, что он самец. Вот график предсказанных значений:

plt.plot(mass_range, fit_ys)

 

decorate(xlabel="Масса (г)", ylabel="P(самец)")

Самые легкие пингвины почти наверняка самки, а самые тяжелые, скорее всего, самцы. Пингвин со средним весом 3750 г примерно равновероятно может быть самцом или самкой.

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

formula = "y ~ mass + flipper_length + culmen_length + culmen_depth"

model = smf.logit(formula, data=adelie)

result = model.fit(disp=False)

display_summary(result)

	coef	std err	z	P > | z |	[0.025	0.975]
Intercept	–60.6075	13.793	–4.394	0.000	–87.642	–33.573
mass	0.0059	0.001	4.153	0.000	0.003	0.009
flipper_length	–0.0209	0.052	–0.403	0.687	–0.123	0.081
culmen_length	0.6208	0.176	3.536	0.000	0.277	0.965
culmen_depth	1.0111	0.349	2.896	0.004	0.327	1.695

Pseudo R-squared: 0.6622

Значение псевдо-R2 у этой модели равно 0.662 — выше, чем у предыдущей модели (0.526). Значит, в добавленных измерениях содержится дополнительная информация, помогающая отличить самцов от самок пингвинов.

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

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

def plot_predictions(mass_range, culmen_length, **options):

    """Построить график предсказанной вероятности в зависимости от массы."""

    df = pd.DataFrame({"mass": mass_range})

    df["flipper_length"] = adelie["flipper_length"].mean()

    df["culmen_length"] = culmen_length

    df["culmen_depth"] = adelie["culmen_depth"].mean()

    fit_ys = result.predict(df)

    plt.plot(mass_range, fit_ys, **options)

Вот как выглядят результаты для трех значений culmen_length — одного стандартного отклонения выше среднего, среднего значения и одного стандартного отклонения ниже среднего:

culmen_length = adelie["culmen_length"]

m, s = culmen_length.mean(), culmen_length.std()

plot_predictions(mass_range, m + s, ls="--", label="Длина клюва выше среднего")

plot_predictions(mass_range, m, alpha=0.5, label="Средняя длина клюва")

plot_predictions(mass_range, m - s, ls=":", label="Длина клюва ниже среднего")

 

decorate(xlabel="Масса (г)", ylabel="P(самец)")

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

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

Чтобы опробовать эту методику, применим эту же модель к пингвинам другого вида. Помимо пингвина Адели, с которым мы работали до сих пор, датасет содержит измерения для 123 особей субантарктического пингвина (Gentoo penguin). Выделим их, воспользовавшись следующей функцией:

def get_species(penguins, species):

    df = penguins.query(f'Species.str.startswith("{species}")').copy()

    df["y"] = (df["Sex"] == "MALE").astype(int)

    return df

gentoo = get_species(penguins, "Gentoo")

len(gentoo)

123

Вот результаты работы модели логистической регрессии:

formula = "y ~ mass + flipper_length + culmen_length + culmen_depth"

model = smf.logit(formula, data=gentoo)

result = model.fit(disp=False)

display_summary(result)

	coef	std err	z	P > | z |	[0.025	0.975]
Intercept	–173.9123	62.326	–2.790	0.005	–296.069	–51.756
mass	0.0105	0.004	2.948	0.003	0.004	0.017
flipper_length	0.2839	0.183	1.549	0.121	–0.075	0.643
culmen_length	0.2734	0.285	0.958	0.338	–0.286	0.833
culmen_depth	3.0843	1.291	2.389	0.017	0.554	5.614

Pseudo R-squared: 0.848

Коэффициент псевдо-R2 составляет 0.848 — больше, чем у пингвинов Адели (0.662). А значит, классификация пингвинов вида субантарктический пингвин с помощью физических измерений даст более верные результаты, чем для пингвинов Адели. Это в свою очередь говорит, что субантарктический пингвин обладает более выраженным половым диморфизмом.

Глоссарий

Регрессия (regression)

Метод оценки коэффициентов, которые подгоняют модель к данным.

Зависимая переменная (response variable)

Переменная, также известная как отклик, которую регрессионная модель пытается предсказать.

Объясняющая переменная (explanatory variable)

Переменная, также известная как независимая переменная, которую модель использует для предсказания зависимых переменных.

Простая регрессия (simple regression)

Регрессия с одной зависимой и одной объясняющей переменной.

Множественная регрессия (multiple regression)

Регрессия с несколькими объясняющими переменными, но только одной зависимой переменной.

Коэффициент (coefficient)

В регрессионной модели коэффициенты — это свободный член (intercept) и коэффициенты наклона (slope) при объясняющих переменных.

Категориальная переменная (categorical variable)

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

Контрольная переменная (control variable)

Переменная, которая добавляется в модель регрессии, чтобы отделить не­посредственный эффект объясняющей переменной от косвенного воздействия.

Обобщенная линейная модель (generalized linear model)

Семейство регрессионных моделей с различными математическими зависимостями между объясняющими переменными и зависимой переменной.

Логистическая регрессия (logistic regression)

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

Упражнения

Упражнение 11.1

Весят ли новорожденные мальчики больше девочек? Чтобы ответить на этот вопрос, снова воспользуемся данными NSFG.

Постройте модель линейной регрессии с totalwgt_lb в качестве зависимой переменной и babysex в качестве категориальной объясняющей переменной: используйте 1 для мальчиков и 2 для девочек. Какова оценочная разница в весе? Является ли она статистически значимой? Если проконтролировать возраст матери, объясняется ли наблюдаемая разница в весе полностью или частично возрастом матери?

Упражнение 11.2

Гипотеза Триверса — Уилларда предполагает, что у многих млекопитающих соотношение полов зависит от «состояния материнской особи», то есть от таких факторов, как возраст, физические параметры, состояние здоровья и социальный статус. Ряд исследований показал наличие этого эффекта у людей, но результаты неоднозначны.

Посмотрим, есть ли связь между возрастом матери и вероятностью рождения мальчика. Постройте модель логистической регрессии, указав пол ребенка в качестве зависимой переменной, а возраст матери — в качестве объясняющей. Шансы матери более старшего возраста родить мальчика повышаются или понижаются? Что получится, если воспользоваться квадратичной моделью возраста матери?

Упражнение 11.3

Для пингвинов Адели постройте модель линейной регрессии, которая предсказывает вес пингвина как функцию переменных flipper_length, culmen_depth, а также Sex в качестве категориальной переменной. Если контролируются длина крыльев и глубина клюва, насколько тяжелее будут самцы пингвинов? Сгенерируйте и постройте диаграмму предсказаний для диапазона длины крыльев у самцов и самок пингвинов, задав значение culmen_depth равным среднему значению.

Упражнение 11.4

Проверим, более или менее диморфен антарктический пингвин (Chinstrap penguin) по сравнению с другими видами, представленными в датасете. Об этом можно судить по коэффициенту псевдо-R2 модели. Используйте функцию get_species, чтобы выбрать пингвинов вида антарктический пингвин, затем постройте модель логистической регрессии с полом в качестве зависимой переменной и всеми четырьмя измерениями в качестве объясняющих переменных. Как значение псевдо-R2 соотносится с другими моделями?

---