Как рассчитать LD50 от блеска?

Question

Я анализирую набор данных, в котором зарегистрировано ~ 10 особей, подвергнутых заданному лечению (Время) и смертности (Alive, Dead). glmer использовался для моделирования данных, так как обработка была заблокирована (пробная версия). Из следующей модели я хочу предсказать Время, в течение которого умирает 50% людей.

Trial <- c(1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3) 
Time <- c(2, 6, 9, 12, 15, 18, 21, 24, 1, 2, 3, 4, 5, 6, 1.5, 3, 4.5, 6, 39) 
Alive <- c(10, 0, 0, 0, 0, 0, 0, 0, 6, 2, 8, 1, 0, 0, 4, 6, 1, 2, 0) 
Dead <- c(0, 10, 6, 10, 10, 10, 7, 10, 0, 8, 1, 9, 10, 10, 5, 0, 8, 6, 10) 
ostrinaA.glmm<- glmer(cbind(Alive, Dead)~Time+(1|Trial), family = binomial(link="logit")) 
summary(ostrinaA.glmm) 

Если бы я был просто моделирования с использованием glm функции dose.p из MASS может быть использован. С другого форума я нашел обобщенный код для дозы.p.glmm от Bill Pikounis. Он выглядит следующим образом:

 dose.p.glmm <- function(obj, cf = 1:2, p = 0.5) { 
     eta <- obj$family$linkfun(p) 
     b <- fixef(obj)[cf] 
     x.p <- (eta - b[1L])/b[2L] 
     names(x.p) <- paste("p = ", format(p), ":", sep = "") 
     pd <- -cbind(1, x.p)/b[2L] 
     SE <- sqrt(((pd %*% vcov(obj)[cf, cf]) * pd) %*% c(1, 1)) 
     res <- structure(x.p, SE = SE, p = p) 
     class(res) <- "glm.dose" 
     res 
    } 

Я новичок в кодировании и мне нужна помощь в настройке этого кода для моей модели. Моя попытка заключается в следующем:

dose.p.glmm <- function(ostrinaA.glmm, cf = 1:2, p = 0.5) { 
    eta <- ostrinaA.glmm$family$linkfun(p) 
    b <- fixef(ostrinaA.glmm)[cf] 
    x.p <- (eta - b[1L])/b[2L] 
    names(x.p) <- paste("p = ", format(p), ":", sep = "") 
    pd <- -cbind(1, x.p)/b[2L] 
    SE <- sqrt(((pd %*% vcov(obj)[cf, cf]) * pd) %*% c(1, 1)) 
    res <- structure(x.p, SE = SE, p = p) 
    class(res) <- "glm.dose" 
    res 
} 
dose.p.glmm(ostrinaA.glmm, cf=1:2, p=0.5) 
Error in ostrinaA.glmm$family : $ operator not defined for this S4 class 

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

Answer 1

На первый взгляд я думаю, что замена

eta <- obj$family$linkfun(p) 

с

f <- family(obj) 
eta <- f$linkfun(p) 

должен сделать трюк.

Кроме того, необходимо заменить res <- ... линии

res <- structure(x.p, SE = matrix(SE), p = p) 

Это довольно туманно, но это необходимо, потому что print.dose.glm метод (от MASS пакета) автоматически пытается cbind() некоторые вещи вместе. Это не удается, если SE является причудливой матрицей из пакета Matrix, а не ванильной матрицей из базы R: matrix() делает преобразование.

Если вы очень новым для кодирования, вы не могли бы понять, что вы не должны изменить имя переменной obj в коде вы скопировали на ostrina.glmm. Другими словами, код Пикуниса должен отлично работать с только двумя предложенными мной модификациями.