Using an offset variable in a Poisson GeneralizedLinearModelFit

Question

In section "2.6 Exposure: Modeling Over Time, Area, and Space" of his book "Modeling Count Data" Joseph M Hilbe uses an offset (exposure) variable in his Poisson regression analysis of the FASTRAK data set (available at Hilbe's download page).

{{"die", "cases", "anterior", "hcabg", "killip", "kk1", "kk2", "kk3", 
  "kk4"}, {5, 19, "Inferior", 0, 4, 0, 0, 0, 1}, {10, 83, "Inferior", 
  0, 3, 0, 0, 1, 0}, {15, 412, "Inferior", 0, 2, 0, 1, 0, 0}, {28, 
  1864, "Inferior", 0, 1, 1, 0, 0, 0}, {1, 1, "Inferior", 1, 4, 0, 0, 
  0, 1}, {0, 3, "Inferior", 1, 3, 0, 0, 1, 0}, {1, 18, "Inferior", 1, 
  2, 0, 1, 0, 0}, {2, 70, "Inferior", 1, 1, 1, 0, 0, 0}, {10, 28, 
  "Anterior", 0, 4, 0, 0, 0, 1}, {9, 139, "Anterior", 0, 3, 0, 0, 1, 
  0}, {39, 443, "Anterior", 0, 2, 0, 1, 0, 0}, {50, 1374, "Anterior", 
  0, 1, 1, 0, 0, 0}, {1, 6, "Anterior", 1, 3, 0, 0, 1, 0}, {3, 16, 
  "Anterior", 1, 2, 0, 1, 0, 0}, {2, 27, "Anterior", 1, 1, 1, 0, 0, 
  0}}

The source, description, and metadata for this this data set is available in Hilbe's MCD Description Data Files: Stata-R-SAS-Excel document. This document also gives the Stata and R-code to analyze this as a proportional intensity Poisson model.

A simple Poisson regression model of the "die" count response variable is easily accomplished in Mathematica (I have recoded the "anterior" variable into a binary variable "ant"; but it can, as well, be entered with the NominalVariables->Automatic option):

glm6 = GeneralizedLinearModelFit[Transpose[{ant, hcabg, kk2, kk3, kk4, 
     die}], {x1, x2, x3, x4, x5}, {x1, x2, x3, x4, x5}, 
   ExponentialFamily -> "Poisson"];
glm6["ParameterTable"]

However, I am unable to succeed with crafting a rate-parameterized Poisson model, where the log of the "cases" variable would be the exposure variable. I do not know how to introduce the cases variable (with something like

LinearOffsetFunction->Log[cases[[#]]]&

) into each record analyzed using the LinearOffsetFunction option of the GeneralizedLinearModelFit command. I tried, unsuccessfully, to use a modified version of the simple template available in the Help documentation:

offsetVar = N[Log[cases]]
glm6 = GeneralizedLinearModelFit[
  Transpose[{ant, hcabg, kk2, kk3, kk4, offsetVar, die}], {x1, x2, x3,
    x4, x5, x6}, {x1, x2, x3, x4, x5, x6}, 
  ExponentialFamily -> "Poisson", LinearOffsetFunction -> (#6 &)]
glm6["ParameterTable"]

The command runs and produces output; but the coefficients do not match the Stata or R results.

Unfortunately I am hampered by both my incomplete knowledge of offset variables in regression statistics and Mathematica. At a minimum, I am not understanding now to enter the cases variable as an offset variable using the LinearOffsetFunction option without adding it to the variable list. The correct solution should return a regression coefficient of 1 with the other coefficients in the model adjusted to account for the influence of the "cases" variable on the "die" count variable, entered as the natural log of "cases" as an offset into the estimating algorithm.

Can the members please assist me with the appropriate syntax?

Answer 1

You'll need to rearrange the data so that the response is found at the end of each row:

data = {{5, 19, "Inferior", 0, 4, 0, 0, 0, 1}, {10, 83, "Inferior", 0,
     3, 0, 0, 1, 0}, {15, 412, "Inferior", 0, 2, 0, 1, 0, 0}, {28, 
    1864, "Inferior", 0, 1, 1, 0, 0, 0}, {1, 1, "Inferior", 1, 4, 0, 
    0, 0, 1}, {0, 3, "Inferior", 1, 3, 0, 0, 1, 0}, {1, 18, 
    "Inferior", 1, 2, 0, 1, 0, 0}, {2, 70, "Inferior", 1, 1, 1, 0, 0, 
    0}, {10, 28, "Anterior", 0, 4, 0, 0, 0, 1}, {9, 139, "Anterior", 
    0, 3, 0, 0, 1, 0}, {39, 443, "Anterior", 0, 2, 0, 1, 0, 0}, {50, 
    1374, "Anterior", 0, 1, 1, 0, 0, 0}, {1, 6, "Anterior", 1, 3, 0, 
    0, 1, 0}, {3, 16, "Anterior", 1, 2, 0, 1, 0, 0}, {2, 27, 
    "Anterior", 1, 1, 1, 0, 0, 0}};
data = data[[All, {2, 3, 4, 5, 6, 7, 8, 9, 1}]];

The GeneralizedLinearModelFit statement should be

glm6 = GeneralizedLinearModelFit[data, {anterior, hcabg, kk2, kk3, kk4},
  {cases, anterior, hcabg, killip, kk1, kk2, kk3, kk4 },
  ExponentialFamily -> "Poisson", NominalVariables -> anterior, 
  LinearOffsetFunction -> (Log[#1] &)];

The parameter table is given by

glm6["ParameterTable"]

Mathematica gives the parameter estimates while Stata(?) gives the "incidence rate ratio" (IRR) which is just the antilog of the parameter estimates:

TableForm[Exp[glm6["BestFitParameters"]], 
 TableHeadings -> {{"Constant", "Anterior", "hcabg", "kk2", "kk3", "kk4"}, None}]