3
$\begingroup$

I'm fitting some data to a Logit model in both Mathematica and R and I'm getting slightly different results.

R code:

data = read.table("http://www.rni.helsinki.fi/~kja/epid12/BCG.dat",header=T)
logitmodel = glm(cbind(D,H)~BCG,data=data,family=binomial(link="logit"))
summary(logitmodel)

Mathematica code:

{d, h, bcg, age} = 
 Transpose@
  Rest@Import["http://www.rni.helsinki.fi/~kja/epid12/BCG.dat", 
    "Data"];
logitmodel = 
 GeneralizedLinearModelFit[Transpose@{bcg, d/h}, BCG, BCG, 
  ExponentialFamily -> "Binomial", LinkFunction -> Automatic, 
  Weights -> d + h]
logitmodel["ParameterTable"]

The results are estimates that are almost the same but differ in a few decimals. For example, the R estimate and standard error for BCG are -0.74152 and 0.12744 respectively. The corresponding Mathematica results are -0.742492 and 0.127755. Am I making a subtle mistake somewhere or is this the result of some numerical approximation?

$\endgroup$
  • 1
    $\begingroup$ I get {-0.7453603285, 0.1271282958} (Win7-64). Probably a numerical accuracy issue, but does it matter? The size of the differences is much smaller than the error in the estimates themselves. $\endgroup$ – Sjoerd C. de Vries Nov 8 '12 at 20:39
  • $\begingroup$ With precision and accuracy set to 200 I get {-0.74536032846461989040..., 0.12712830428988635844565...} $\endgroup$ – Sjoerd C. de Vries Nov 8 '12 at 20:45
  • $\begingroup$ @Sjoerd That's a valid point. However, subtle mistakes can occur in statistical analysis, so whenever a difference--no matter how small--cannot definitively be assigned to numerical imprecision, it ought to be tracked down. Kudos to Mr Alpha for being careful. $\endgroup$ – whuber Nov 8 '12 at 21:43
11
$\begingroup$

To match what you're doing in R, you need to use d/(d+h) for the fractions rather than d/h:

logitmodel = GeneralizedLinearModelFit[Transpose@{bcg, d/(d + h)}, BCG, BCG, 
  ExponentialFamily -> "Binomial", LinkFunction -> Automatic, Weights -> d + h];
logitmodel["ParameterTable"]

Parameter table

The clues that led to this realization were that the deviances, log likelihoods, and AICs could not be made to match, even after trying to make glm.fit more precise in R (and starting the optimization at the Mathematica values for the parameters): this indicated there really was something different between the two calculations. Once it became apparent--upon inspecting the data--that h was huge compared to d, the source of the difference was evident.

$\endgroup$
  • $\begingroup$ Thank you. It was one of those mistakes that seems obvious in hindsight. $\endgroup$ – Mr Alpha Nov 9 '12 at 8:31

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.