我找不到一个简单的函数来做到这一点。里面有一些代码predict
函数取决于拟合模型(例如确定模型的排名)。但是,我们可以创建一个函数来制作一个可以与预测一起使用的假 glm 对象。这是我第一次尝试这样的功能
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
and 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()
。您需要传递一个数据帧,以便 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
在这里,我使用旧数据和我指定的系数对原始模型和新模型进行预测。我们可以看到概率的估计发生了一些变化。
现在我还没有彻底测试这个功能,所以使用时需要您自担风险。我没有做尽可能多的错误检查。也许其他人确实知道更好的方法。