`ddply`はグループ別に私のデータセットにロジスティック回帰（GLM）を適用できません

Question

私はパッケージを使用して、異なる実験の複数の集団についてLD50（致死量）を計算しています。私はデータをサブセット化して一度に1つずつ行うだけで十分ですが、ddplyを使用するとエラーが発生します。基本的には、各温度で各集団にLD50が必要です。

私のデータはいくぶん次のようになります。

# 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)) 

次正常に動作します：

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] 

しかし、私はddplyにこれを入れたときに、私は次のエラーを取得する：

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

I ca私は割合を使用して私は私のアプローチ（と、この質問を書いていた）と間違っている場所を把握することはできません右LD50を取得します。

Answer 1

おそらくこれはあなたを驚かせるでしょう。しかし、あなたが式

cbind(Alive, Dead) ~ Dose 

代わりの

Mortality ~ Dose 

を使用することを選択した場合、問題が消えてしまいます。

---

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 

---

何がMortality ~ Doseで起こったのか？ ddplyを呼び出すときのが.inform = TRUEを設定してみましょう：

## `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 

は、今、私たちは、私たちは、変数Mortalityが次元を失っている、とだけ最初の列（Alive）が保持されていることがわかります。 glmがbinomialの応答の場合、応答が単一のベクトルである場合、glmは0-1バイナリまたは2レベルの係数を期待します。今、私たちは、整数20、9、8、6を、持っている...、それゆえglmは

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

この問題を解決する方法は本当にありません文句を言うでしょう。私はプロテクターを使用しようとしました：

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

しかし、それでもやはり同じ失敗で終わります。

Answer 2

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

ベースRソリューション：

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ソリューション：

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