Я работаю над докторской диссертацией и пытаюсь подобрать обобщенную линейную модель смешанных эффектов, используя пакет «MCMCglmm» в R. Я неоднократно читал полезные учебные материалы Джаррода.Тем не менее, в спецификации приоров все еще есть некоторые проблемы, которые я не смог выяснить.
В моем исследовании зависимой переменной является количество людей, участвующих в программе энергосбережения в домашних хозяйствах (результат подсчета).Он был неоднократно измерен для каждого дня в течение приблизительно трех лет для каждого из 360 сообществ (таким образом, данные достаточно велики и n = 371, 520).Кроме того, эти сообщества расположены в разных районах (всего 90 районов).Таким образом, данные продольного суточного счета вложены в сообщества внутри районов.Мое исследование направлено на выяснение того, какие факторы будут влиять на (ежедневное) количество участников такой программы.Базовая модель (сверхдисперсная) модель Пуассона:
#give the priors
prior.poi <- list(R = list(V = diag(1), nu = 0.002, n=0, fix=1),
G = list(
G1=list(V = diag(3)*0.02, nu =4),
G2=list(V=diag(3)*0.02, nu=4)
))
#fit the model
model.poi <- MCMCglmm(y ~ 1 + t + x + x:t + t2 + t3 + t4 + c1 + c2 + c3 + d1 + d2 + d3,
random = ~ us(1 + t + x):no_c + us(1 + t + x):no_d,
rcov = ~idh(1):units,
family = "poisson",
data = dat.big,
prior = prior.poi,
burnin = 15000, nitt = 65000, thin = 50)
В части с фиксированными эффектами, y - результат подсчета;'t' измеряет время в истекших днях с начала программы;«х» - еще одно поведенческое вмешательство, реализованное для некоторых сообществ«t2 ~ t4» - другие факторы времени (т. е. манекены, измеряющие выходные и праздничные дни, а также логарифм среднесуточной температуры);«c1 ~ c3» и «d1 ~ d3» измеряют характеристики уровня сообщества и района соответственно, которые являются переменными, не зависящими от времени (например, общая численность населения, размер района).В части случайных эффектов «no_c» и «no_d» - это номер записи каждого сообщества и района.
Поскольку в моих данных много избыточных нулей, поэтому я в дальнейшем использую препятствие (чрезмерно рассредоточенное)Модель Пуассона выглядит следующим образом.
#give the priors
prior.hp <- list(R = list(V = diag(2), nu = 0.002, n=0, fix=1),
G = list(
G1=list(V = diag(6)*0.02, nu =7),
G2=list(V=diag(6)*0.02, nu=7)
)
)
#fit the model
model.hp <- MCMCglmm(y ~ -1 + trait + trait:t + trait:x + trait:x:t + trait:t2 + trait:t3 + trait:t4 + trait:c1 + trait:c2 + trait:c3 + trait:d1 + trait:d2 + trait:d3,
random = ~ us(trait + trait:t + trait:x):no_c + us(trait + trait:t + trait:x):no_d,
rcov = ~idh(trait):units,
family = "hupoisson",
data = dat.big,
prior = prior.hp,
burnin = 15000, nitt = 65000, thin = 50,
pr = T, pl = T)
Обе модели и модели Пуассона с препятствиями могут работать хорошо только тогда, когда 'fix = 1' было добавлено в R-структуру предыдущей спецификации.Когда он был удален из априорных таблиц, обе модели возвращали сообщение об ошибке:
«Уравнения смешанной модели единственного числа: используйте (более сильный) предшествующий»
и останавливайтеБег.Эта ошибка не исчезнет, независимо от того, были ли использованы расширения параметров в G-структуре (то есть alpha.mu = rep (0, 3), alpha.V = diag (3) * 25 ^ 2 для модели Пуассона OD иalpha.mu = rep (0, 6), alpha.V = diag (6) * 25 ^ 2 для модели препятствий) или нет, независимо от того, были ли удалены / отрегулированы другие элементы в R-структуре.
В модели Пуассона с препятствиями, поскольку ковариационная матрица для процесса с нулевым изменением не может быть оценена, в R-структуре следует использовать «fix = 2», а не «fix = 1».Однако модель не может работать хорошо, если остаточная дисперсия для усеченного до нуля процесса Пуассона не будет установлена на 1, как описано выше.Мой вопрос заключается в том, уместно ли фиксировать остаточную дисперсию для процессов с нулевым изменением и для процессов Пуассона в 1 в R-структуре?Это слишком «информативно» для моей оценки модели?Есть ли другие приоры, которые я могу использовать, чтобы модель работала хорошо?Любая идея по этим вопросам приветствуется.