`ddply` n'applique pas la régression logistique (GLM) par groupe à mon jeu de données

Question

Je travaille sur la DL50 (dosage létal) pour plusieurs populations d'expériences différentes en utilisant le paquet MASS. C'est assez simple quand je sous-ensemble les données et faire un à la fois, mais je reçois une erreur quand j'utilise ddply. Essentiellement, j'ai besoin d'une DL50 pour chaque population à chaque température.

Mes données ressemble un peu à ceci:

# dput(d) 
d <- structure(list(Pop = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
3L), .Label = c("a", "b", "c"), class = "factor"), Temp = structure(c(1L, 
1L, 1L, 1L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 1L, 
1L, 1L, 1L, 2L, 2L, 2L, 2L), .Label = c("high", "low"), class = "factor"), 
Dose = c(1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 
1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L), Dead = c(0L, 
11L, 12L, 14L, 2L, 16L, 17L, 7L, 5L, 3L, 17L, 15L, 9L, 20L, 
8L, 19L, 7L, 2L, 20L, 14L, 9L, 15L, 1L, 15L), Alive = c(20L, 
9L, 8L, 6L, 18L, 4L, 3L, 13L, 15L, 17L, 3L, 5L, 11L, 0L, 
12L, 1L, 13L, 18L, 0L, 6L, 11L, 5L, 19L, 5L)), .Names = c("Pop", 
"Temp", "Dose", "Dead", "Alive"), class = "data.frame", row.names = c(NA, 
-24L)) 

Ce qui suit fonctionne très bien:

d$Mortality <- cbind(d$Alive, d$Dead) 
a <- d[d$Pop=="a" & d$Temp=="high",] 
library(MASS) 
dose.p(glm(Mortality ~ Dose, family="binomial", data=a), p=0.5)[1] 

Mais quand je mets cela en ddply je reçois l'erreur suivante:

library(plyr) 
d$index <- paste(d$Pop, d$Temp, sep="_") 
ddply(d, 'index', function(x) dose.p(glm(Mortality~Dose, family="binomial", data=x), p=0.5)[1]) 

Error in eval(expr, envir, enclos) : y values must be 0 <= y <= 1

Je ca n obtenir la bonne DL50 quand j'utilise une proportion mais ne peux pas comprendre où je me suis trompé avec mon approche (et avait déjà écrit cette question).

Answer 1

Peut-être vous étonnera. Mais si vous choisissez d'utiliser la formule

cbind(Alive, Dead) ~ Dose 

au lieu de

Mortality ~ Dose 

le problème aura disparu.

---

library(MASS) 
library(plyr) 

## `d` is as your `dput` result 

## a function to apply 
f <- function(x) { 
    fit <- glm(cbind(Alive, Dead) ~ Dose, family = "binomial", data = x) 
    dose.p(fit, p=0.5)[[1]] 
    } 

## call `ddply` 
ddply(d, .(Pop, Temp), f) 

# Pop Temp  V1 
#1 a high 2.6946257 
#2 a low 2.1834099 
#3 b high 2.5000000 
#4 b low 0.4830998 
#5 c high 2.2899553 
#6 c low 2.5000000 

---

Qu'est-il arrivé avec Mortality ~ Dose? Fixons .inform = TRUE lorsque vous appelez ddply:

## `d` is as your `dput` result 
d$Mortality <- cbind(d$Alive, d$Dead) 

## a function to apply 
g <- function(x) { 
    fit <- glm(Mortality ~ Dose, family = "binomial", data = x) 
    dose.p(fit, p=0.5)[[1]] 
    } 

## call `ddply` 
ddply(d, .(Pop, Temp), g, .inform = TRUE) 

#Error in eval(expr, envir, enclos) : y values must be 0 <= y <= 1 
#Error: with piece 1: 
# Pop Temp Dose Dead Alive Mortality 
#1 a high 1 0 20  20 
#2 a high 2 11  9   9 
#3 a high 3 12  8   8 
#4 a high 4 14  6   6 

Maintenant, nous nous voyons cette variable Mortality a perdu dimension, et seule la première colonne (Alive) est conservée. Pour une réponse glm avec binomial, si la réponse est un vecteur unique, glm attend 0-1 binaire ou un facteur de deux niveaux. Maintenant, nous avons des nombres entiers 20, 9, 8, 6, ..., d'où glm se plaindra

Error in eval(expr, envir, enclos) : y values must be 0 <= y <= 1

Il n'y a vraiment aucun moyen de résoudre ce problème. J'ai essayé d'employer un protecteur:

d$Mortality <- I(cbind(d$Alive, d$Dead)) 

mais il finit toujours avec la même panne.

Answer 2

library(MASS) 
d$Mortality <- do.call("cbind", list(d$Alive, d$Dead)) 

base de R Solution:

ld50 <- sapply(X = unique(d$Pop), 
       FUN = function(x){ 
       df1 = subset(d, Pop == x) 
       vec1 <- sapply(X = unique(df1$Temp), 
           FUN = function(y){ 
            dose.p(glm(Mortality ~ Dose, family="binomial", data = subset(df1, Temp == y)), p=0.5)[1] 
           }) 
       names(vec1) <- unique(df1$Temp) 
       return(vec1) 
      }) 
colnames(ld50) <- unique(d$Pop) 
ld50 
#    a   b  c 
# high 2.694626 2.5000000 2.289955 
# low 2.183410 0.4830998 2.500000 

purrr Solution:

library(purrr)  
ld50 <- map(.x = unique(d$Temp), 
      .f = ~{ 
       df1 = subset(d, Temp == .x) 
       vec1 <- map_dbl(.x = unique(df1$Pop), 
           .f = ~{ 
           dose.p(glm(Mortality ~ Dose, family="binomial", data = subset(df1, Pop == .x)), p=0.5)[1] 
           }) 
       names(vec1) <- unique(df1$Pop) 
       return(vec1) 
      }) 
names(ld50) <- unique(d$Temp) 
do.call('rbind', ld50) 
#    a   b  c 
# high 2.694626 2.5000000 2.289955 
# low 2.183410 0.4830998 2.500000