2013-11-08 4 views
2

Я пытаюсь использовать модели регрессионной модели смешанных эффектов, но я обеспокоен размерами переменных выборок в каждом кластере/группе, а также очень низким числом «успехов», в некоторых моделях.Переменный размер выборки для каждого кластера/группы в смешанных эффектах Логистическая регрессия

У меня есть ~ 700 деревьев, распределенных по 163 полям (т. Е. Группа/группа), ежегодно посещаемых с 2004 по 11. Я устанавливаю отдельные модели логистической регрессии смешанных эффектов (далее GLMM) для каждого года исследования, чтобы сравнить этот результат с выводом из общей модели уязвимости (т. Е. Анализ выживаемости со случайным эффектом).

Количество деревьев на участок варьируется от 1-22. Кроме того, в течение нескольких лет очень мало «успехов» (т. Е. Больных деревьев). Например, в 2011 году было всего 4 успеха из 694 «неудач» (т. Е. Здоровых деревьев).

Мои вопросы: (1) существует общее правило для идеального числа образцов | группа, когда фокус вывода - только для оценки фиксированных эффектов в GLMM, и (2) являются стабильными GLMM, когда существует такое крайняя разница в соотношении успехов: неудачи.

Благодарим за любые советы или предложения источников.

Сара

+1

Потому что это (хорошо) вопрос так очень мало содержания программ, было бы лучше спросил более на HTTP://stats.stackexchange.com/. Просто не забудьте удалить его здесь, прежде чем переместить его, чтобы он не попал в перекрестные сообщения. –

+0

Сердечно согласен с Джошем. Я думаю, что предлагаю некоторый предварительный R-код, когда вы повторно разместите, все равно сделайте это более «конкретным». (4/694) не является особенно низкой частотой событий для модели Пуассона. –

ответ

0

Как сказал Джош, это лучше вопросы для CrossValidated.

Не существует жестких правил логистической регрессии, но одно эмпирическое правило - 10 успехов и 10 отказов необходимы для каждой ячейки в дизайне (кластер в этом случае), умноженное на число непрерывных переменных в модели.

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

1

(Привет, Сара, жаль, что я не ответил ранее по электронной почте ...)

Трудно ответить на эти вопросы в общем - вы застряли с вашими данными, не так ли? Так что дело не в анализе власти. Если вы хотите, чтобы ваши результаты были разумно надежным, вероятно, самое лучшее, что нужно сделать, это запустить некоторые симуляции. Я собираюсь показать довольно недавнюю функцию lme4 (в версии разработки 1.1-1 на Github), которая должна имитировать данные из GLMM с учетом формулы и набора параметров.

Сначала я должен имитировать предикторов (вы не бы сделать это, так как у вас уже есть данные - хотя вы можете попробовать варьируя диапазон числа участков, деревьев на участке , и т.д.).

set.seed(101) 
## simulate number of trees per plot 
## want mean of 700/163=4.3 trees, range=1-22 
## by trial and error this is about right 
r1 <- rnbinom(163,mu=3.3,size=2)+1 
## generate plots and trees within plots 
d <- data.frame(plot=factor(rep(1:163,r1)), 
      tree=factor(unlist(lapply(r1,seq)))) 
## expand by year 
library(plyr) 
d2 <- ddply(d,c("plot","tree"), 
     transform,year=factor(2004:2011)) 

Теперь настроить параметры: Я буду считать год фиксированный эффект и что в целом заболеваемость составляет plogis(-2)=0.12 кроме в 2011 году, когда он plogis(-2-3)=0.0067.Среди-участка стандартного отклонения равно 1 (по шкале логита), как это имеет место среди-дерев внутри-участка стандартного отклонения:

beta <- c(-2,0,0,0,0,0,0,-3) 
theta <- c(1,1) ## sd by plot and plot:tree 

Теперь имитации: год, как фиксированный эффект, участок и дерево-в пределах -plot как случайные эффекты

library(lme4) 
s1 <- simulate(~year+(1|plot/tree),family=binomial, 
    newdata=d2,newparams=list(beta=beta,theta=theta)) 
d2$diseased <- s1[[1]] 

Обобщить/проверка:

d2sum <- ddply(d2,c("year","plot"), 
      summarise, 
      n=length(tree), 
      nDis=sum(diseased), 
      propDis=nDis/n) 
library(ggplot2) 
library(Hmisc) ## for mean_cl_boot 
theme_set(theme_bw()) 
ggplot(d2sum,aes(x=year,y=propDis))+geom_point(aes(size=n),alpha=0.3)+ 
    stat_summary(fun.data=mean_cl_boot,colour="red") 

Теперь установите модель:

g1 <- glmer(diseased~year+(1|plot/tree),family=binomial, 
     data=d2) 
fixef(g1) 

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