Understanding where the link function is computed in Poisson regression

Viewed 83

From Poisson Regression by hand this 'manual' Poisson coefficient function is provided:

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

Can any one tell me exactly where the link function is computed here? What would this code look like with an identity link function? I have basic understanding from YouTube videos etc, but no one explains the actual computation.

How would this code look like with offset and weights?

1 Answers

The GLM is specified as:

or equivalently as

In the case of the Poisson model, the canonical link is , so

You can see that in your code where lambda = exp(beta %*% t(x)) If you wanted to estimate a model with an identity link, you would use lambda = beta %*% t(x) because in that case

Addressing the follow-up questions in the comments, with weights, often these are used to weight the likelihood contribution for each observation. So, imagine you had a variable wt that you passed into your function LogLike(), you could modify it by multiplying the likelihood contribution by the weight:

LL <- -sum(dpois(y, lambda, log = TRUE)*wt)

An offset is just a variable whose coefficient is forced to be 1 (see this post for a discussion). Lets say you had a variable called offset that you wanted to include in your model. Again, you could pass that argument to LogLike() and you would need to modify the function by:

  lambda <- exp(beta%*%t(x) + offset)

Following the example in the CrossValidated post I linked, offset = log(time). Which offset is correct is more of a theoretical or substantive matter.

Related