Предпосылки: я подбираю смешанную логистическую модель с двумя фиксированными эффектами и двумя перекрещенными случайными эффектами. Две переменные, составляющие фиксированные эффекты, являются двоичными категориальными переменными (т. Е. Числовых предикторов нет).
Данные взяты с айтрекера, и каждый объект имеет тысячи нулей и единиц с течением времени. Одна из моделей, которую я хотел рассмотреть, - это логистическая смешанная модель, учитывающая корреляцию нулей и единиц внутри человека и внутри изображения, на которое человек смотрит в данный момент времени (в течение всего эксперимента они смотрят на шесть разных изображений, из 12 возможных). (Кстати, я планирую также проанализировать эти данные с помощью анализа кривой роста, но мне также было любопытно попробовать это, временно игнорируя динамику времени.)
Проблема: я запустил модель, используя фиктивное кодирование обоих двоичных предикторов с уровнями базовых показателей по умолчанию для каждого. Модель сошлась. Это также имело смысл на основе графиков данных, которые я уже сделал: оценки были в ожидаемых направлениях, значимость фиксированных эффектов была такой, как ожидалось, и т. Д.
Затем я хотел подогнать модель, параметризованную с одним из предикторов на его другом уровне базовой линии, чтобы получить связанное значение p (см., Например: https://stats.stackexchange.com/questions/449723/calculating-simple-effects-from-standard-regression-output). Я использовал relevel() для этого, а затем снова установил ту же самую модель. Однако на этот раз было предупреждение о сближении. Я попробовал еще раз (почему бы и нет), и снова ничего не вышло. Я не могу придумать причину, по которой базовый уровень должен влиять на конвергенцию (или ее отсутствие).
Структура данных: вот небольшой фиктивный набор данных, показывающий общую структуру данных (значительно уменьшенную по сравнению с фактическим набором данных, который имеет ~ 5000 наблюдений (нулей и единиц) на человека. на изображение, а общая длина составляет более 1 миллиона строк (в фактическом наборе данных ~ 60 объектов)). Обратите внимание: Time_stamp не входит в модель, так как я не рассматриваю временной ход в этом анализе. AOI в AOI_look обозначает область интереса, общий термин в данных отслеживания взгляда.
structure(list(Subject = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), Image = c("A", "A",
"A", "B", "B", "B", "C", "C", "C", "D", "D", "D", "E", "E", "E",
"F", "F", "F", "G", "G", "G", "E", "E", "E", "D", "D", "D", "A",
"A", "A", "J", "J", "J", "L", "L", "L"), Time_stamp = c(0L, 1L,
10L, 0L, 8L, 16L, 0L, 7L, 16L, 0L, 1L, 10L, 0L, 8L, 16L, 0L,
7L, 16L, 0L, 8L, 16L, 0L, 7L, 16L, 0L, 1L, 9L, 0L, 8L, 16L, 0L,
8L, 16L, 0L, 1L, 9L), Intervention = c("Pre", "Pre", "Pre", "Pre",
"Pre", "Pre", "Pre", "Pre", "Pre", "Post", "Post", "Post", "Post",
"Post", "Post", "Post", "Post", "Post", "Pre", "Pre", "Pre",
"Pre", "Pre", "Pre", "Pre", "Pre", "Pre", "Post", "Post", "Post",
"Post", "Post", "Post", "Post", "Post", "Post"), Sex = c("Female",
"Female", "Female", "Female", "Female", "Female", "Female", "Female",
"Female", "Female", "Female", "Female", "Female", "Female", "Female",
"Female", "Female", "Female", "Male", "Male", "Male", "Male",
"Male", "Male", "Male", "Male", "Male", "Male", "Male", "Male",
"Male", "Male", "Male", "Male", "Male", "Male"), AOI_look = c(1L,
1L, 1L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 0L, 0L, 1L, 1L, 1L, 0L, 0L,
1L, 1L, 0L, 0L, 0L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 1L, 1L, 0L, 1L,
1L, 1L, 1L)), class = "data.frame", row.names = c(NA, -36L))
Модель 1: Вмешательство = Публикация в качестве базового показателя: предупреждение о конвергенции отсутствует
logit_model01 <- glmer(AOI_look ~ factor(Sex)*factor(Intervention) +
(1 | Image) +
(1 | Subject),
data = data01,
family="binomial")
Выход модели 1 (конвергентный):
summary(logit_model01)
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) [glmerMod]
Family: binomial ( logit )
Formula: AOI_look ~ factor(Sex) * factor(Intervention) + (1 | Image) +
(1 | Subject)
Data: data01
AIC BIC logLik deviance df.resid
1519023.6 1519095.5 -759505.8 1519011.6 1182007
Scaled residuals:
Min 1Q Median 3Q Max
-2.6200 -0.9733 0.5504 0.8444 2.2317
Random effects:
Groups Name Variance Std.Dev.
Subject (Intercept) 0.2320 0.4816
Image (Intercept) 0.1318 0.3631
Number of obs: 1182013, groups: Subject, 59; Image, 9
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.052133 0.036172 1.441 0.15
factor(Sex)Male 0.187736 0.039967 4.697 2.64e-06 ***
factor(Intervention)Pre 0.257262 0.006136 41.924 < 2e-16 ***
factor(Sex)Male:factor(Intervention)Pre -0.220240 0.007932 -27.766 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) fc(S)M fc(I)P
fctr(Sex)M -0.006
factr(Int)Pr 0.011 0.012
f(S)M:(I)P -0.013 -0.018 -0.747
Модель 2: то же, что и Модель 1, но Вмешательство = Предварительно было установлено на базовое значение перед запуском: получить предупреждение о конвергенции:
logit_model02 <- glmer(AOI_look ~ factor(Sex)*factor(Intervention) +
(1 | Image) +
(1 | Subject),
data = data01,
family="binomial")
Предупреждение о сходимости:
Warning message:
In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
Model failed to converge with max|grad| = 0.00210935 (tol = 0.002, component 1)
Хотя было предупреждение, я все еще могу получить результат для Модели 2 с summary():
summary(logit_model02)
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) [glmerMod]
Family: binomial ( logit )
Formula: AOI_look ~ factor(Sex) * factor(Intervention) + (1 | Image) +
(1 | Subject)
Data: data01
AIC BIC logLik deviance df.resid
1519023.6 1519095.5 -759505.8 1519011.6 1182007
Scaled residuals:
Min 1Q Median 3Q Max
-2.6200 -0.9733 0.5504 0.8444 2.2317
Random effects:
Groups Name Variance Std.Dev.
Subject (Intercept) 0.2320 0.4816
Image (Intercept) 0.1318 0.3631
Number of obs: 1182013, groups: Subject, 59; Image, 9
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.309400 0.038106 8.120 4.68e-16 ***
factor(Sex)Male -0.032522 0.042520 -0.765 0.444
factor(Intervention)Post -0.257262 0.006158 -41.774 < 2e-16 ***
factor(Sex)Male:factor(Intervention)Post 0.220234 0.007926 27.787 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) fc(S)M fc(I)P
fctr(Sex)M -0.107
fctr(Int)Pst 0.009 0.024
f(S)M:(I)P -0.013 -0.025 -0.748
convergence code: 0
Model failed to converge with max|grad| = 0.00210935 (tol = 0.002, component 1)
Для ясности / полноты добавлена информация, найденная в комментариях под ответом Бена Болкера:
Модель 3: то же, что и Модель 1, но для агрегированных данных (нули и единицы суммированы в успехи и неудачи по Тема и Изображение < / em> комбинация):
Результаты не такие же, как у Модели 1 (соответствие неагрегированным данным 0/1), хотя в Модели 1 не было предупреждения о конвергенции. А именно, все точечные оценки аналогичны, но не для всех SE: SE для взаимодействия и вмешательства (внутри субъектов) аналогичны, но для пола (между -subjects) и перехват заметно отличаются:
logit_model03 <- glmer(cbind(Successes, Failures) ~ factor(Sex)*factor(Intervention) +
(1 | Image) +
(1 | Subject),
data = data01_agg,
family="binomial")
summary(logit_model03)
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.052151 0.150564 0.346 0.729
factor(Sex)Male 0.187734 0.125563 1.495 0.135
factor(Intervention)Pre 0.257261 0.006227 41.314 <2e-16 ***
factor(Sex)Male:factor(Intervention)Pre -0.220239 0.008061 -27.322 <2e-16 ***
Чтобы увидеть, может ли явный размер необработанного набора данных (более 1,67 миллиона строк) вызывать проблему, я взял случайную выборку из 15000 наблюдений из моего исходного набора данных и обновил модель, используя как неагрегированные (0/1), так и агрегированные (успех, неудача) данные. (Я сделал это для обоих уровней базового уровня для Вмешательства, но я просто показываю это для Вмешательства = Опубликовать здесь.):
Модель 4: случайная выборка из 15000 наблюдений. из фактического набора данных, неагрегированный (0/1) результат, вмешательство = Опубликовать как базовый:
logit_model04 <- glmer(AOI_look ~ factor(Sex)*factor(Intervention) +
(1 | Image) +
(1 | Subject),
data = smaller_raw_data,
family="binomial")
summary(logit_model04)
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.03212 0.15564 -0.206 0.836505
factor(Sex)Male 0.28384 0.14121 2.010 0.044429 *
factor(Intervention)Pre 0.22727 0.06449 3.524 0.000425 ***
factor(Sex)Male:factor(Intervention)Pre -0.25061 0.08410 -2.980 0.002884 **
Модель 5: случайная выборка из 15000 наблюдений. из фактического набора данных, агрегированный результат (успех, неудача), вмешательство = Опубликовать как базовый:
logit_model05 <- glmer(cbind(Successes, Failures) ~ factor(Sex)*factor(Intervention) +
(1 | Image) +
(1 | Subject),
data = smaller_agg_data,
family="binomial")
summary(logit_model05)
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.03211 0.15569 -0.206 0.836582
factor(Sex)Male 0.28383 0.14126 2.009 0.044500 *
factor(Intervention)Pre 0.22727 0.06450 3.524 0.000425 ***
factor(Sex)Male:factor(Intervention)Pre -0.25061 0.08412 -2.979 0.002890 **
Как видите, точечные оценки и SE для фиксированных эффектов для моделей 4 и 5 одинаковы, по крайней мере, до пятого знака после запятой. Я также повторил это для Intervention = Pre as baseline, и результаты эквивалентны в отношении их сходства в неагрегированных и агрегированных моделях (как я и ожидал).
Таким образом, использование неагрегированных или агрегированных данных не должно иметь значения в теории, по крайней мере, для набора данных разумного размера (т. Е. Не слишком большого). Это заставляет меня задаться вопросом, есть ли нестабильность в моделях, когда набор данных становится слишком большим. Согласно ?lme4::convergence, есть основания полагать, что это могло быть так: «Исследовательский анализ показывает, что (1) наивная оценка Гессе может потерпеть неудачу для больших наборов данных (количество наблюдений больше, чем приблизительно 1e5)». Так что, возможно, это корень того, что я наблюдал в моделях 1-3 выше.
В моем случае у меня есть возможность агрегировать данные, поскольку оба моих предиктора являются бинарными категориальными, но для числовых предикторов, которые не позволяют агрегирование, это не было бы вариантом для уменьшения вычислительной нагрузки до подбора модели.
(1|Subject:Trial)в качестве случайного эффекта и посмотрю, решит ли это проблему. Обратите внимание, что этот тот же сценарий, что и stackoverflow.com/questions/63360751/, я только что упомянул, что это глаз данные отслеживания, поэтому результат называется AOI_look вместо Binary_outcome. В очередной раз благодарим за помощь. (Кроме того, следует ли перенести этот вопрос и на CrossValidated, или здесь все в порядке? Для меня этот вопрос более четко касался кодирования.) - person Meg schedule 12.08.2020Postв качестве базовой линии. ), но значение p для Sex сильно отличается (0,135 против 2,64e-06, приведенного выше), как и для члена перехвата. Это заставляет меня задаться вопросом, не истекает ли тайм-аут lmer () (или, в последней модели выше, проблема сходимости) из-за размера необработанных данных (~ 1,6 миллиона строк) и, таким образом, дает неточный вывод. Работаем над этим ... Отчитаюсь. - person Meg schedule 13.08.2020