Как исправить перехват в биномиальной модели GLM? - PullRequest
Новая
       102

Как исправить перехват в биномиальной модели GLM?

0 голосов
/ 19 декабря 2018

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

intercept <- 1.0
fit <- lm(I(x - intercept) ~ 0 + y, lin)
summary(fit)

Или указав смещение: 

m(x ~ y +0 + offset(rep(1, nrow(lin))),

Но я не могу понять, как это сделать с помощью биномиальной регрессии.Пример кода ниже: 

library(ggplot2)
library(lme4)
library(dplyr)
library(gather)

# simulate dummy dataframe

x <- data.frame(time = c(1, 1, 1, 1, 1, 1,1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2,3, 3, 3, 3, 3, 3, 3, 3, 3,4, 4, 4, 4, 4, 4, 4, 4, 4),
                type = c('a', 'a', 'a', 'b', 'b', 'b','c','c','c','a', 'a', 'a', 'b', 'b', 'b','c','c','c','a', 'a', 'a', 'b', 'b', 'b','c','c','c','a', 'a', 'a', 'b', 'b', 'b','c','c','c'), randef=c('aa','ab','ac','ba','bb','bc','ca','cb','cc','aa','ab','ac','ba','bb','bc','ca','cb','cc','aa','ab','ac','ba','bb','bc','ca','cb','cc','aa','ab','ac','ba','bb','bc','ca','cb','cc'), 
                surv = sample(x = 1:200, size = 36, replace = TRUE), 
                nonsurv= sample(x = 1:200, size = 36, replace = TRUE))


# convert to survival and non survival into individuals following 
# https://stackoverflow.com/questions/51519690/convert-cbind-format-for- binomial-glm-in-r-to-a-dataframe-with-individual-rows

x_long <- x %>%
  gather(code, count, surv, nonsurv)

# function to repeat a data.frame
x_df <- function(x, n){
    do.call('rbind', replicate(n, x, simplify = FALSE))
    }

# loop through each row and then rbind together
x_full <- do.call('rbind', 
            lapply(1:nrow(x_long), 
                   FUN = function(i) 
                     x_df(x_long[i,], x_long[i,]$count)))

# create binary_code
x_full$binary <- as.numeric(x_full$code == 'surv')


### binomial glm with interaction between time and type:
summary(fm2<-glm(binary ~ time*type, data = x_full, family = "binomial"))

### plot glm in ggplot2
ggplot(x_full, aes(x = time, y = as.numeric(x_full$binary), fill= x_full$type)) +
 geom_smooth(method="glm", aes(color = factor(x_full$type)), method.args = list(family = "binomial"))

То, что я хотел бы сделать, это форсировать перехват через 1 (т.е. все выжившие в первой временной точке), почти как кривая выживания.Задание смещения, по-видимому, не работает в этом примере, так как это биноминальная вероятность, поэтому указание 1 или 0 не работает, и число выживших против не выживших людей в первый момент времени варьируется среди факторов).

Переполнение стека, я узнал от вас больше, чем кто-либо из моих студентов и аспирантов.Есть предложения по этому поводу?

...