3
$\begingroup$

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.

enter image description here

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?

$\endgroup$
5
$\begingroup$

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"]

Parameter estimates

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}]

incidence risk rate

$\endgroup$
  • $\begingroup$ Thank you for clarifying this syntax. I was unaware that one could include variables in the variable list that do not enter into the model. Now that you have shown me a practical example, I understand the help bullet explanation for LinearOffsetfunction->h and how the exposure variable enters the model (and from where). $\endgroup$ – BeanSith Mar 6 at 17:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.