Установите GEE-модель типа «сменный» с гаммой - PullRequest
Новая
       122

Установите GEE-модель типа «сменный» с гаммой

0 голосов
/ 23 сентября 2019

Я хотел бы оценить плавный эффект некоторого ковариата N в маргинальной модели типа "обмениваемой" в R, где переменная кластеризации равна S.Из того, что я мог найти, это должно быть возможно с:

geeglm(..., id = S, corstr = "exchangeable")

, а также:

gamm(..., correlation = corCompSymm(form = ~1|S))

Ниже вы можете найтипример, где результаты выглядят хорошо в том смысле, что две оценки довольно близки.Однако, если я использую реальные данные, о которых говорит наш проект, предполагаемые эффекты сглаживания, как правило, сильно отличаются.Я не могу опубликовать это здесь, но, возможно, кто-то все еще может заметить некоторую проблему в коде.Например (см. Ниже), gamm -объект говорит Number of Groups: 1, что меня беспокоит, поскольку там явно более одного кластера ...

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

########
## Packages
########
library(ggplot2)
library(mgcv)
library(dplyr)
library(geepack)
library(splines)

########
## Data Simulation
########
f        <- function(N) {return((-200+(N-25)^2)/100)}

N        <- sort(sample(1:50, 10, replace = T))
S        <- as.character(1:10)
S_Effect <- rnorm(length(S),0,1)
S_Effect <- rep(S_Effect,N)
S        <- rep(S,N)
N        <- rep(N,N)
E        <- runif(length(N))

data     <- data.frame(O        = rep(0,length(N)),
                       E        = E,
                       N        = N,
                       S        = as.factor(S),
                       S_Effect = S_Effect)

for (i in 1:length(N)) {
  data$O[i] <- rbinom(1, 1, plogis(f(N[i]) + qlogis(E[i]) + S_Effect[i]))}

data <- data %>% mutate(E = qlogis(E))



########
## Fitting
########
formula_gamm   <- as.formula("O ~ 1 + offset(E) + s(N, bs = 'bs')")
model_gamm     <- gamm(formula_gamm, family = binomial(), correlation = corCompSymm(form=~1|S), data = data)
model_gamm

formula_geeglm <- as.formula("O ~ 1 + offset(E) + bs(N)")
model_geeglm   <- geeglm(formula_geeglm, family = binomial(), corstr = "exchangeable", id = S, data = data)



########
## Plot
########
pred_gamm      <- plot.gam(model_gamm$gam, select = 1)
x <- pred_gamm[[1]]$x
pred_geeglm  <- predict(model_geeglm, type = "terms", newdata = data.frame(E = rep(0,length(x)), N = x))

z                 <- qnorm(0.9)

tmp               <- data.frame(x = x,
                               y = pred_gamm[[1]]$fit,
                               group = rep("estimate gamm",length(x)))
tmp2               <- data.frame(x = x,
                                y = as.numeric(pred_geeglm),
                                group = rep("estimate geeglm",length(x)))
tmp3              <- data.frame(x = x,
                               y = f(x),
                               group = rep("actual function",length(x)))

data_pred = bind_rows(tmp,tmp2,tmp3) %>% mutate(group = as.factor(group))

p <- ggplot(data = data_pred, aes(x = x, y = y, color = group)) +
     geom_line(size = 2) +
     xlab("N") +
     ylab("f(N)")
p

Дополнительный вопрос: gamm -объект содержит достаточно информации, чтобы построить доверительный интервал-полосы вокруг оценочной функции, но как я могу сделать это для geeglm -метки?Вы получаете что-то, что выглядит разумно, если вы simulate(model_geeglm, ...) и берете точечное mean и так далее, но это меня не очень устраивает, поскольку (1) документация по simulate не упоминает маргинальные модели и (2)это очень примитивно ...

1 Ответ

0 голосов
/ 24 сентября 2019

GAMM использует штрафованные сплайны, так что степени свободы, используемые результирующим сплайном (сглаживающим), вероятно, будут несколько меньше, чем запрошенный базисный размер, равный 10. GEE подходит для непенализованной модели.При прочих равных, модель без наказания будет более волнистой, чем модель, оштрафованная.

Чтобы сравнить эти подходы на общих основаниях, необходимо убедиться, что оба параметра bs() и s(x, bs = 'bs') дают одинаковое количествобазисные функции (версия s() может выдавать на единицу меньше, так как устраняет отсутствие идентификации с термином перехвата, тогда как в версии bs() вы пропускаете перехват).

Убедившись, что вы получитетот же базисный размер, тогда вы можете настроить GAMM для непенализованного сплайна, добавив fx = TRUE к s(...) члену в формуле.

Сделав это, обе модели должны оценивать одинаковые сглаживающие эффекты.

Однако я бы посоветовал вам использовать наказание;Для модели GAMM используйте fx = FALSE, а затем после оценки модели запустите gam.check(model$gam) (заменив model на ваш подобранный объект модели) и посмотрите, пройдет ли проверка базового размера для сглаживания.

...