手动实现泊松回归

Question

手动实现泊松回归

rregressionestimationpoissonlog-likelihood

3

我想手动进行泊松回归，并定义一个函数，用于估计任意数量的系数。我有两个问题：

第一：如何获得beta矩阵，而不必显式地编写每个beta。我想以以下方式编写lambda 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))
}

- Dima Ku

1

“但它不起作用”..你能更具体一些吗？你是否收到了错误信息？ - MrSmithGoesToWashington

抱歉，我的错。我用更好的代码和两个明确的问题编辑了我的问题。 - Dima Ku

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- MrSmithGoesToWashington · Accepted Answer

以下是一个基于您的代码的工作示例，但不包括解释变量的平方：

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))