`ddply` 无法按组将逻辑回归 (GLM) 应用于我的数据集

Question

我正在使用 MASS 包计算来自不同实验的多个人群的 LD50（致死剂量）。当我对数据进行子集化并一次做一个时，这很简单，但是当我使用 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

当我使用比例时，我可以获得正确的 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）。对于具有 binomial 响应的 glm，如果响应是单个向量，则 glm 需要 0-1 二进制或两个级别的因子。现在，我们有整数 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))

但还是以同样的失败告终。

`ddply` 无法按组将逻辑回归 (GLM) 应用于我的数据集

`ddply` fails to apply logistic regression (GLM) by group to my dataset

regression

r

plyr

glm

logistic-regression