小议非二元响应变量的EC50

其实解决这个问题的R包还是很多的，MASS里就有一个dose.p来计算LD50，我们来仔细看看其中的代码：

dose.p <- function(obj, cf = 1:2, p = 0.5) {
    eta <- family(obj)$link(p)
    b <- coef(obj)[cf]
    x.p <- (eta - b[1])/b[2]
    names(x.p) <- paste("p = ", format(p), ":", sep = "")
    pd <- -cbind(1, x.p)/b[2]
    SE <- sqrt(((pd %*% vcov(obj)[cf, cf]) * pd) %*% c(1, 1))
    structure(x.p, SE = SE, p = p, class = "glm.dose")
}
print.glm.dose <- function(x, ...) {
    M <- cbind(x, attr(x, "SE"))
    dimnames(M) <- list(names(x), c("Dose", "SE"))
    x <- M
    NextMethod("print")
}

前面一个eta求解的是广义线性模型在生存概率p下的响应值，b取的是模型参数，说白了就是给定响应求浓度。回想一下EC50的定义，其实就是先建模后预测一个数值。至于说SE，也是通过模型参数来计算的。所以看来看去整个问题的核心就在于广义线性模型的建立。一般而言LD50因为响应变量为二元的，所以要进行logistic变换，之后如果自变量浓度跨度较大，要进行log变换。

但在EC50的计算过程中，响应变量通常不是二元的，如果检测器给出的是连续变量，那实际上这个问题就略复杂一点，也就是把响应变量缩小到[0,1]范围内,直接应用线性模型对剂量进行线性回归，而EC50就是线性模型截距跟剂量参数比值的负数，但要是求SE就还得是MASS里面的方法，也就是通过模型参数的SE去做估计。这里用的变换其实就是常见的四参数求解EC50的方法，实际情况中用这个的较多，但本质上就是一个广义线性模型中的logistic回归或probit回归。

说到底，在ligistic回归中对响应变量是有明确要求为二元服从二项分布的，但实际应用中的响应值却是转化（多为极值范围限定到[0,1]空间）过去的，这就是很多疑惑的起源。其实对于S型曲线还有其他种类的子模型，有的更有针对性，如Hill equation。另外需要注意的是该模型的适用范围，最好先做图看看，如果混合了其他过程如生长，胁迫，模型就需要修改一下。此外，也可以尝试局部回归与非参的方法，这对于估计EC50这一推断值的标准误或方差非常有效，而这却是经验所不能达的。