我想手工进行泊松回归，并定义一个函数，该函数可用于估计任意数量的系数。我有两个问题：
首先：如何才能获得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))