在R中手动构建逻辑回归模型以进行预测

Question

我正在尝试测试一个 logistic 回归模型（例如，3 个预测变量 X1、X2、X3 的 3 个系数），并使用数据集进行测试。我知道如何在创建模型对象后测试模型，例如：

mymodel <- glm( Outcome ~  X1 + X2 + X3 , family = binomial,data=trainDat)

然后测试数据

prob <- predict(mymodel,type="response",newdata=test)

但是我想现在使用我已经拥有的系数和截距创建一个物流模型，然后在数据上测试这个模型。

基本上，我不清楚如何在不运行glm的情况下创建“mymodel”。

问题的背景：我使用偏移量运行了逻辑回归。

mymodel <- glm(Outcome ~ offset(C1 * X1) + offset(C2 * X2) + X3, 
               family = binomial, data = trainDat)

因此，mymodel对象生成了一个只有截距（I）和C3系数（针对特征X3）的模型。
现在我需要在测试数据集上测试完整模型（即I + C1*X1 + C2*X2 + C3*X3），但是我不知道如何获得完整模型，因为mymodel的输出仅具有截距和C3。因此，我认为我的更一般性问题是：“如何手动构建逻辑回归模型对象？”
感谢您的耐心等待。

- nerdlyfe

---

那么你只是想保存 mymodel 对象并在以后重新加载它，这样你就不必再次运行 glm 了吗？它运行起来需要很长时间或者有什么问题吗？ - MrFlick

他们还包括了方差/协方差矩阵吗？以及所有因子水平的编码？ - MrFlick

不，不需要协方差矩阵。是的，我知道因子水平编码。协方差矩阵不应该是不必要的吗？我知道你可以在Excel中手动使用系数和因子水平规则制作逻辑函数。这里，我会添加一些内容到我的问题，以允许多种答案类型。 - nerdlyfe

包含测试数据的问题将会产生更多兴趣。 - IRTFM

Answer 1

我找不到一个简单的函数来完成这个任务。在predict函数中有一些代码依赖于拟合模型（例如确定模型的等级）。但是，我们可以创建一个函数来制作一个虚假的glm对象，您可以使用它与predict一起使用。下面是我尝试创建这样一个函数的第一步：

makeglm <- function(formula, family, data=NULL, ...) {
    dots <- list(...)
    out<-list()
    tt <- terms(formula, data=data)
    if(!is.null(data)) {
        mf <- model.frame(tt, data)
        vn <- sapply(attr(tt, "variables")[-1], deparse)

        if((yvar <- attr(tt, "response"))>0)
            vn <- vn[-yvar]
            xlvl <- lapply(data[vn], function(x) if (is.factor(x))
           levels(x)
        else if (is.character(x))
           levels(as.factor(x))
        else
            NULL)
        attr(out, "xlevels") <- xlvl[!vapply(xlvl,is.null,NA)]
        attr(tt, "dataClasses") <- sapply(data[vn], stats:::.MFclass)
    }
    out$terms <- tt
    coef <- numeric(0)
    stopifnot(length(dots)>1 & !is.null(names(dots)))
    for(i in seq_along(dots)) {
        if((n<-names(dots)[i]) != "") {
            v <- dots[[i]]
            if(!is.null(names(v))) {
                coef[paste0(n, names(v))] <- v
            } else {
                stopifnot(length(v)==1)
                coef[n] <- v
            }
        } else {
            coef["(Intercept)"] <- dots[[i]]
        }   
    }
    out$coefficients <- coef
    out$rank <- length(coef)
    out$qr <- list(pivot=seq_len(out$rank))
    out$family <- if (class(family) == "family") {
        family
    } else if (class(family) == "function") {
        family()
    } else {
        stop(paste("invalid family class:", class(family)))
    }
    out$deviance <- 1
    out$null.deviance <- 1
    out$aic <- 1
    class(out) <- c("glm","lm")
    out
}

因此，该函数创建一个对象并传递所有预期在该对象上找到的predict和print值。现在我们可以测试它。首先，这是一些测试数据

set.seed(15)
dd <- data.frame(
    X1=runif(50),
    X2=factor(sample(letters[1:4], 50, replace=T)),
    X3=rpois(50, 5),
    Outcome = sample(0:1, 50, replace=T)
)

我们可以使用标准的二项式模型进行拟合

mymodel<-glm(Outcome~X1+X2+X3, data=dd, family=binomial)

这提供了

Call:  glm(formula = Outcome ~ X1 + X2 + X3, family = binomial, data = dd)

Coefficients:
(Intercept)           X1          X2b          X2c          X2d           X3  
    -0.4411       0.8853       1.8384       0.9455       1.5059      -0.1818  

Degrees of Freedom: 49 Total (i.e. Null);  44 Residual
Null Deviance:      68.03 
Residual Deviance: 62.67    AIC: 74.67

现在假设我们想在相同的数据上尝试一下我们在出版物中读到的模型。这里是如何使用makeglm函数。

newmodel <- makeglm(Outcome~X1+X2+X3, binomial, data=dd, 
    -.5, X1=1, X2=c(b=1.5, c=1, d=1.5), X3=-.15)

第一个参数是模型的公式。这定义了响应和协变量，就像在运行 glm 时一样。接下来，您需要像使用 glm() 一样指定family（分布类型）。并且您需要传递一个数据框，以便R可以嗅探涉及的每个变量的正确数据类型。这也将使用data.frame标识所有因子变量及其级别。因此，这可以是编码方式与拟合数据框相同的新数据，也可以是原始数据。

现在，我们开始指定要在模型中使用的系数。系数将使用参数名称填充。未命名参数将用作截距。如果您有一个因子，则需要通过命名参数为所有水平指定系数。在这里，我决定将拟合估计值舍入为“好看”的数字。

现在我可以使用我们的newmodel进行预测。

predict(mymodel, type="response")
#         1         2         3         4         5
# 0.4866398 0.3553439 0.6564668 0.7819917 0.3008108

predict(newmodel, newdata=dd, type="response")

#         1         2         3         4         5
# 0.5503572 0.4121811 0.7143200 0.7942776 0.3245525

在这里，我使用原始模型以及使用旧数据和指定系数的新模型进行预测。我们可以看到概率估计结果略有变化。

现在，我并没有彻底测试这个函数，因此使用时请自行承担风险。我并没有进行足够的错误检查，可能还有更好的方法，希望有其他人能提供。

- MrFlick

---

那是一次惊人的努力，我欠你一个人情。我非常感激这个，也知道其他人也会感激。我会使用并审查它，并尝试在下面添加一些进一步的见解。 - nerdlyfe

我还将代码发布为gist。我稍后会在那里发布我对代码所做的任何更新。 - MrFlick