Подход к применению логарифмических линейных моделей на категориальных данных с приложениями в экологических исследованиях.
Как специалисты по науке о данных, мы вооружаем себя стандартными инструментами из нашего пояса, которые мы склонны использовать, когда мы формулируем проблемную область как проблему классификации или регрессии (для обучения с учителем) или проблему кластеризации (для обучения без учителя). . Эти инструменты дают довольно хорошие результаты при внедрении решений машинного обучения применительно к структурированным данным.
В этой статье я обсуждаю надежный инструмент моделирования для работы с категориальными данными: Обобщенная линейная модель (GLM).
GLM широко используются в социальных науках, маркетинговых исследованиях, медицинских науках, где преимущественно используются категориальные шкалы. В сфере финансовых услуг страховая отрасль также широко использует GLM. В маркетинговых исследованиях, например, исследование потребительских предпочтений, проведенное компанией, занимающейся исследованием рынка, на их панели клиентов между выбором бренда для шампуня, может попросить их указать свой выбор между брендами, произведенными L'Oreal, P&G и Unilever. Моделирование здесь используется для определения того, как объяснить выбор или предпочтение клиента на основе других переменных (также называемых ковариатами), каждая из которых дополнительно представлена в категориальных шкалах.
Что такое категориальные данные?
Рассмотрим следующие примеры —
Упорядоченные категориальные данные: опрос обслуживания клиентов, показывающий, что клиенты оценивают услугу по выбору {плохо, хорошо, отлично}. В этом списке категорий есть определенный порядок от низшей к высшей.
Неупорядоченные категориальные данные: в исследовании потребительских предпочтений демографическая переменная того, где живет потребитель, может быть переменной. Варианты могут включать {городской, пригородный, сельский и т. д.}. Для этой опции нет понятия порядка.
GLM — отличный инструмент для работы с исследованиями, требующими анализа и моделирования на основе данных этого типа.
Что такое обобщенная линейная модель?
Прежде чем мы перейдем к тому, как, давайте кратко рассмотрим, что такое GLM.
GLM — это широкий набор моделей, включающий модели линейной регрессии и ANOVA для непрерывных данных в качестве переменных отклика, а также модели для дискретных откликов. Одним из таких частных случаев модели дискретного ответа для бинарного ответа является модель логистической регрессии. Этот метод моделирования используется для моделирования кредитных и других рисков в сфере финансовых услуг, поскольку он известен своей эффективностью, простотой и экономичностью применительно к этой области. GLM также позволяют моделировать переменную отклика на основе количества (частоты).
GLM состоит из трех компонентов:
1. Случайная составляющая (Y), в которой мы предполагаем наличие определенного вероятностного распределения.
2. Систематический компонент, который представляет ковариаты или независимые переменные для модели.
3. Функция связи, представляющая функцию ожидаемого значения Y. И эта функция позволяет выразить Y как некоторую линейную комбинацию независимых переменных.
Моделирование и анализ с помощью GLM
Далее давайте рассмотрим анализ и моделирование с использованием GLM вместе с кодом.
В этом анализе мы используем небольшой набор данных из области экологии, где была собрана информация о дневной среде обитания двух видов ящериц — Grahami и Opalinus. Эти данные были получены путем наблюдения конкретной информации о местах обитания и времени их наблюдения (McCullagh and Nelder, 1989).
Давайте посмотрим на данные, чтобы убедиться, что мы понимаем, как выглядит весь набор данных категориального типа данных.
Виды ящериц (L): Grahami или Opalinus, указывающие на то, какой вид наблюдался. Это наша целевая переменная.
Освещение (I): указывает, было ли это Солнце или Тень.
Диаметр насеста (D): в дюймах со следующими значениями (также известными как коэффициенты): D≤ 2 и D›2.
Высота насеста (H): в дюймах (H) со следующими значениями: H‹5 и H≥ 5.
Время суток (T): также факторная переменная с тремя уровнями: Раннее, Полуденное и Позднее.
Как мы видим, непрерывных данных нет, все они являются категориальными данными.
Это моделирование было реализовано в R. Вы можете реализовать это на своем любимом языке, включая Python. Концепция применения категориального анализа данных и моделирования одинакова.
Посмотрим на данные:
> lizard.data <-read.table(file = "lizard_data.txt",header = TRUE) > head(lizard.data) Illumination Diameter Height Time Species Count 1 Sun <=2 <=5 Early Grahami 20 2 Sun <=2 <=5 Early Opalinus 2 3 Sun <=2 <=5 Midday Grahami 8 4 Sun <=2 <=5 Midday Opalinus 1 5 Sun <=2 <=5 Late Grahami 4 6 Sun <=2 <=5 Late Opalinus 4
Данные настраиваются в виде так называемой таблицы непредвиденных обстоятельств более высокого порядка. Мы применяем различные логарифмические линейные модели к этим данным. Мы рассчитываем размер эффекта, осуществляем выбор лучшей модели на основе AIC. Мы также проводим тесты на соответствие и остаточный анализ для оценки адекватности модели.
Шаги, предпринятые для этого анализа, вывода и обсуждения производительности модели, подробно описаны ниже.
Подбор логарифмических моделей
Используя данные о ящерицах, мы сначала подогнали три логарифмически-линейные модели, в том числе:
(i) Модель взаимной независимости,
(ii) Модель гомогенной ассоциации,
(iii) Модель, содержащая все условия трехфакторного взаимодействия.
# Create the features from the dataset
I<-factor(lizard.data[,1])
D<-factor(lizard.data[,2])
H<-factor(lizard.data[,3])
T<-factor(lizard.data[,4])
L<-factor(lizard.data[,5])
Count<-factor(lizard.data[,6])
library(MASS) #need this package for using loglm
##create the three models first
#create saturated model
liz.model.sat <- glm(Count ~ I*D*H*T*L, family = poisson, data=lizard.data)
liz.model.sat$formula
#mutual independence
liz.model.mi <- glm(Count ~I+D+H+T+L, family = poisson, data=lizard.data)
summary(liz.model.mi)
#homogeneous association
liz.model.ha <-glm(Count ~ I*D+I*H+I*T+I*L+D*H+D*T+D*L+H*T+H*L+T*L,
family = poisson, data=lizard.data)
summary(liz.model.ha)
## then, get LRT stat (deviance, or, g-sq), aic, pearson x-sq, p-value for model
#mutual independence:
deviance(liz.model.mi) #152.6136
df.residual(liz.model.mi) # 41
AIC(liz.model.mi) #329.8624
mi.1 <- loglm(Count ~ I + D + H + T + L, data = lizard.data)
x2.mi<-mi.1$pearson
x2.mi #pearson chi sq is 157.8766
anova(liz.model.mi, liz.model.sat, test = "Chisq")
#p-value is 9.186e-15 ***
#homogenous association
deviance(liz.model.ha)#5.05375
df.residual(liz.model.ha) # 27
AIC(liz.model.ha) #230.3026
ha.1 <- loglm( Count ~ I * D + I * H + I * T + I * L + D * H +
D * T + D * L + H * T + H * L + T * L, family = poisson,
data = lizard.data)
x2.ha<-ha.1$pearson
x2.ha #pearson chi sq is 20.93927
Затем мы вычисляем качество подгонки, используя остаточное отклонение, которое совпадает со статистикой квадрата G (хорошесть подгонки), показанной в таблице 1 ниже.
Использование AIC для определения лучшей модели
Чтобы получить хорошее соответствие, нам нужна насыщенная модель, модель с df = 0. В отличие от (биномиальной) логистической регрессии, мы не можем использовать только AIC из-за возможной условной связи между двумя или более терминами в логарифмической линейной модели. Мы начинаем оценку с рассмотрения статистики выбора модели, включая качество соответствия (статистика LRT или отклонение), p-значение, AIC, значение хи-квадрата Пирсона и остаточные степени свободы.
Усложнение моделирования
Мы оцениваем три модели, созданные в предыдущем разделе, то есть однофакторные термины, двухфакторные термины и все трехфакторные термины, и дополнительно настраиваем четырехфакторную модель на основе терминов (код ниже). Глядя на строки с 1 по 4 (таблица 2), модель взаимной независимости подходит очень плохо.
## create a three factor interaction model, and get diagnostics
liz.model.trif <-glm(Count ~
I*D*H+I*D*T+I*D*L+I*H*T+I*H*L+I*T*L+D*H*T+D*H*L+D*T*L+H*T*L,
family = poisson, data=lizard.data)
summary(liz.model.trif)
deviance(liz.model.trif)#13.23126
df.residual(liz.model.trif) # 11
AIC(liz.model.trif) #250.4801
trif.1 <- loglm( Count ~ I * D * H + I * D * T + I * D * L + I *
H * T + I * H * L + I * T * L + D * H * T + D * H * L + D *
T * L + H * T * L, family = poisson,
data = lizard.data)
x2.trif<-trif.1$pearson
x2.trif #pearson chi sq is 11.85212
anova(liz.model.trif, liz.model.sat, test = "Chisq")
#p-value is 0.2785
#create four factor interaction model
liz.model.fourf <- glm(Count ~ I*D*H*T+D*H*T*L,family = poisson, data=lizard.data)
# #Illumination*Species*Height*Time possibly has a zero in cross-class tab, so term removed
summary(liz.model.fourf)
deviance(liz.model.fourf) # 14.63944
df.residual(liz.model.fourf) # 12
AIC(liz.model.fourf) # 249.8883
fourf.1 <- loglm( Count ~ I * D * H * T + D * H * T * L, family = poisson,
data = lizard.data)
x2.fourf<-fourf.1$pearson
x2.fourf #pearson chi sq is NaN.
anova(liz.model.fourf, liz.model.sat, test = "Chisq")
#p-value is 0.2617
Слишком сложное сокращение: использование StepAIC для оценки
Нам нужна модель, более сложная, чем модель однородной ассоциации, но более простая, чем трехфакторные и четырехфакторные модели взаимодействий. Чтобы получить это, мы использовали функцию StepAIC, начав с четырехфакторной оценки моделей с разными AIC. Мы также проделали тот же шаг с трехфакторным термином. модель. Мы выбрали модели в обоих случаях с самым низким AIC в каждом из них, а затем приступили к удалению избыточных членов в каждом (выбор моделей показан в строках 5 и 6 в таблице 2).
Основываясь на этом наблюдении, мы выбираем строку 6 (D,T,IH,IT,DH,DT,DL,TL,IHL) как лучшую модель.
## setup multiple models, evaluate each using AIC and goodness of fit
#use stepAIC to manually select few models that 'appear' good for evaluation.
stepAIC(liz.model.mi) #329.9, residual dev = 152.6
stepAIC(liz.model.ha) #225.7, residual dev = 28.44,
stepAIC(liz.model.trif)#AIC: 227.4, residual dev =26.1.
stepAIC(liz.model.fourf)#AIC: 249.9, residual dev=14.64
##pick the lowest residual models and remove redundant params to fit few models
##and choose one final best model
#model 1
liz.model.besteval.1 <- glm(formula = Count ~ I + I:H + D:H + I:T + D:T +
+ D:L + T:L + I:H:L+ D:H:T:L, family = poisson, data = lizard.data)
#I:D:H:T + (fails to converge...so, term removed)
aic1<-AIC(liz.model.besteval.1)
aic1 #238.254
d1<-deviance(liz.model.besteval.1)
d1 # 15.00517
df.residual(liz.model.besteval.1) #18
anova(liz.model.besteval.1, liz.model.sat, test = "Chisq") #0.6616
besteval11 <- loglm(Count ~ I + I:H + D:H + I:T + D:T +
I:L + D:L + H:L + T:L + I:H:L + D:H:T:L, family = poisson, data = lizard.data)
x2.bestevall11 <-besteval11$pearson
x2.bestevall11 #12.58482
#model 2
liz.model.besteval.2 <- glm(formula = Count ~ D + T + I:H + D:H + I:T + D:T +
D:L + T:L + I:H:L, family = poisson, data = lizard.data)
aic2<-AIC(liz.model.besteval.2)
aic2 #227.3522
d2<-deviance(liz.model.besteval.2) #26.1033
d2
df.residual(liz.model.besteval.2) #29
anova(liz.model.besteval.2, liz.model.sat, test = "Chisq")#p-value is 0.62
besteval21 <- loglm(Count ~ D + T + I:H + D:H + I:T + D:T +
+ D:L + T:L + I:H:L, family = poisson, data = lizard.data)
x2.bestevall21 <-besteval21$pearson
x2.bestevall21#25.28304
liz.best=liz.model.besteval.2 #best model.
deviance(liz.best)
Остаточный анализ лучшей модели
Статистические данные о соответствии, рассчитанные выше, показывают модель, которая, по-видимому, достаточно хорошо соответствует данным. Важно, чтобы остатки оценивались в каждой ячейке, чтобы показать качество соответствия. Остаточный анализ может показать, что некоторые ячейки могут демонстрировать несоответствие модели, которая в противном случае могла бы хорошо соответствовать данным. Мы рассчитываем стандартизированный остаток с помощью функции rstandard, где наблюдаемые и подогнанные числа делятся на их стандартные ошибки, показанные в таблице 3 ниже.
Вывод таблицы 3 был сгенерирован с использованием следующего кода:
#get residuals for best model pred.best<-fitted(liz.best) infl<-influence.measures(liz.best) #residual diagnostics influencePlot(liz.best) influenceIndexPlot(liz.best) std.resid<-rstandard(liz.best, type = "pearson") # standardized resid likeli.resid<-rstudent(liz.best, type = "pearson") #likelihood resid liz.resid<-data.frame(I,D,H,T,L,Count, pred.best, std.resid, likeli.resid) write.csv(file="lizard_residuals_pearson.csv", liz.resid)
Мы выполняем несколько проверок модели, включая influencePlot. На рисунке ниже показано несколько наблюдений в виде выбросов. Стьюдентизированное остаточное значение, большее чем равное 2 или меньшее или равное -2, считается выбросом. Мы видим две клетки, всего 5 ящериц, как выброс. Размер пузырька указывает на меру влияния наблюдения (это квадратный корень из расстояния Кука).
Представить это как проблему логистической регрессии
Мы также можем выполнить логистическую регрессию, где нам потребуется преобразовать набор данных из 48 записей в набор данных из 564 записей, где каждая запись представляет виды GrahamiилиOpalinusвиды . Это позволило бы нам предсказать для заданных условий I, D, H, T, какие виды ящериц мы могли бы найти. Чтобы использовать этот метод, мы преобразуем ящерицу в качестве переменной ответа.
Данные были преобразованы в длинную таблицу с 564 наблюдениями, и каждая переменная была перекодирована 0,1 для всех переменных и 1,2,3 для T. Записи заголовка выглядят следующим образом:
library(car) # for using Anova
lizard.longdata <-read.csv("lizard_final_binomial.csv")
lizard.longdata
head(lizard.longdata)
lizard.longdata$L
liz.logistic <- glm(lizard.longdata$L ~ ., family = binomial, data=lizard.longdata)
#effects
logistic.LRT.stats<-Anova(liz.logistic) #default test.stat is LRT, LRT stats
logistic.LRT.stats
cbind(liz.logistic$coefficients,confint(liz.logistic),exp(confint(liz.logistic)))
anova(liz.logistic)
AIC(liz.logistic)
liz.logistic.sat <-glm(lizard.longdata$L ~ I*P*D*T, family = binomial, data=lizard.longdata)
Мы используем функцию Anova из пакета car, чтобы получить статистику LRT, которая показана ниже.
Эффекты показаны в таблице ниже для 95% доверительного интервала. Последние два столбца показывают экспоненциальный CI.
Лучшее подмножество с использованием AIC
Мы применили функцию stepAIC к модели с полной переменной, в результате чего мы получили ту же модель. Другими словами, это был тот, у которого было самое низкое значение AIC 577,62 и остаточное отклонение 567,6216.
L ~ I + P + D + T
Как мы видим на графике влияния ниже, остатки на лучшей полученной модели намного выше, чем то, что мы получаем на логлинейной модели (выше график влияния для логлинейной модели).
Мы делаем вывод, что эта модель плохо соответствует данным.
Заключительные слова об эквивалентности логистической и лог-линейной моделей
Было бы логично, если бы эти две модели были эквивалентны, мы также знаем, что логистическая модель не знает условий взаимодействия. Однако в этом анализе мы не видим эквивалентности, так как девиация двух моделей совершенно различна. Мы получили 26,1033 для лучшей логарифмической линейной модели и отклонение 555,4047 для лучшей модели логистической регрессии.
Мы пришли к выводу, что логарифмическая линейная модель хорошо выражает отношения между видами и их переменными, которые помогают нам их различать.
Есть много замечательных тем, в которые стоит углубиться при работе с категориальными данными. Я расскажу о них в будущих постах. Я надеюсь, что этот анализ был полезен для понимания того, как работать с GLM.
Хотите продолжить разговор?
Если вы не подписаны на мои электронные письма, сделайте это. Вот пять причин, почему вы можете захотеть: www.vijayreddiar.com/email
И подписывайтесь на меня в Твиттере [@ReddiarVijay]
До скорого!
Использованная литература:
Маккалла П., Нелдер Дж. А., 1989. Обобщенные линейные модели. Лондон, Чепмен и Холл, 511 стр.
Агрести, А. 2007. Введение в категориальный анализ данных. Уайли.
https://stats.stackexchange.com/questions/26930/residuals-for-logistic-regression-and-cooks-distance
http://users.stat.ufl.edu/~presnell/Links/R-links.shtml
