手工泊松回归 [英] Poisson Regression by hand
本文介绍了手工泊松回归的处理方法,对大家解决问题具有一定的参考价值,需要的朋友们下面随着小编来一起学习吧!
问题描述
我想手工进行泊松回归,并定义一个函数,该函数可用于估计任意数量的系数。我有两个问题:
首先:如何才能获得Beta矩阵,而不必编写每个Beta显式。我想这样编写lambda=exp(t(X)%*%beta)。我以为我可以做一个for循环,为x、a、beta中的每一列创建一个for,然后在矩阵中将它们相加ab,但我不知道如何编码。
第二个: 因为我不知道怎么写,所以我试着写了估算6个贝塔的函数。我得到了数据集翘曲破缺的结果,但系数与GLM中的不同,为什么?我也不知道必须粘贴哪些值才能进行解析,也不知道如果不将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:", "
")
#
print(x$call)
#
cat("
")
# Coefficients
cat("Coefficents:", "
")
#
Koef <- (t(x$coefficients))
#
rownames(Koef) <- ""
#
print(round(Koef, 3))
}
推荐答案
这里是一个基于您的代码的工作示例。但不包括解释变量的平方:
LogLike <- function(y,x, par) {
beta0 <- par[1]
beta1 <- par[2]
beta2 <- par[3]
beta3 <- par[4]
# the deterministic part of the model:
lambda <- exp(beta0*x[,1] + beta1 * x[,2] +beta2*x[,3]+beta3*x[,4])
# 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){
# # definiere Regressionsformel
form <- formula(formula)
#
# # dataFrame wird erzeugt
model <- model.frame(formula, data = data)
#
# # Designmatrix erzeugt
x <- model.matrix(formula,data = data)
#
# # Response Variable erzeugt
y <- model.response(model)
par <- c(0,0,0,0)
erg <- list(optim(par=par,fn=LogLike,x=x,y=y)$par)
return(erg)
}
PoisMod(breaks~wool+tension, as.data.frame(daten))
您可以与GLM进行比较:
glm(breaks~wool+tension, family = "poisson", data = as.data.frame(daten))
编辑:任意数量的解释变量
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){
# # definiere Regressionsformel
form <- formula(formula)
#
# # dataFrame wird erzeugt
model <- model.frame(formula, data = data)
#
# # Designmatrix erzeugt
x <- model.matrix(formula,data = data)
#
# # Response Variable erzeugt
y <- model.response(model)
par <- rep(0,ncol(x))
erg <- list(optim(par=par,fn=LogLike,x=x,y=y)$par)
return(erg)
}
PoisMod(breaks~wool+tension, as.data.frame(daten))
glm(breaks~wool+tension, family = "poisson", data = as.data.frame(daten))
这篇关于手工泊松回归的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!
查看全文