我想手动进行泊松回归并定义一个可用于估计任意数量的系数的函数。我有两个问题:
第一:我怎样才能得到一个 beta 矩阵,而不必明确地写出每个 beta。我想以这种方式编写 lambda lambda = exp(t(x)%*%beta) 。我以为我可以做一个 for 循环并为 xa beta 中的每一列创建并将它们 ab 求和到一个矩阵中,但我不知道如何编码。
第二:由于我不知道如何为 i beta 编写代码,因此我尝试编写用于估计 6 个 beta 的函数。我得到了数据集 warpbreaks 的结果,但系数与 glm 中的不同,为什么?我也不知道我必须将哪些值粘贴到 par,也不知道如果我不将 x 和 y 粘贴到函数中,为什么 optim 不起作用。
希望你能帮忙!
daten <- warpbreaks
LogLike <- function(y,x, par) {
beta <- par
# the deterministic part of the model:
lambda <- exp(beta%*%t(x))
# and here comes the negative log-likelihood of the whole dataset, given the
# model:
LL <- -sum(dpois(y, lambda, log = TRUE))
return(LL)
}
PoisMod<-function(formula, data){
# #formula
form <- formula(formula)
#
# # dataFrame
model <- model.frame(formula, data = data)
#
# # Designmatrix
x <- model.matrix(formula,data = data)
#
# # Response Variable
y <- model.response(model)
par <- rep(0,ncol(x))
call <- match.call()
koef <- optim(par=par,fn=LogLike,x=x,y=y)$par
estimation <- return(list("coefficients" = koef,"call"= call))
class(result) <- "PoisMod"
}
print.PoisMod <- function(x, ...) {
# Call
cat("Call:", "\n")
#
print(x$call)
#
cat("\n")
# Coefficients
cat("Coefficents:", "\n")
#
Koef <- (t(x$coefficients))
#
rownames(Koef) <- ""
#
print(round(Koef, 3))
}