手动泊松回归
Poisson Regression by hand
我想手动进行泊松回归并定义一个可用于估计任意数量系数的函数。我有 2 个问题:
首先:我怎样才能得到一个 beta 矩阵,而不必明确地写出每个 beta。我想以这种方式编写 lambda = exp(t(x)%*%beta) 。我以为我可以做一个 for 循环并为 x 中的每一列创建一个 beta,然后将它们求和到一个矩阵中,但我不知道如何编码。
第二个:
因为我不知道如何为 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))
}
这是一个工作示例,基于您的代码..但没有解释变量的平方:
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))
我想手动进行泊松回归并定义一个可用于估计任意数量系数的函数。我有 2 个问题:
首先:我怎样才能得到一个 beta 矩阵,而不必明确地写出每个 beta。我想以这种方式编写 lambda = exp(t(x)%*%beta) 。我以为我可以做一个 for 循环并为 x 中的每一列创建一个 beta,然后将它们求和到一个矩阵中,但我不知道如何编码。
第二个: 因为我不知道如何为 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))
}
这是一个工作示例,基于您的代码..但没有解释变量的平方:
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))