# Getting slightly different results in fitting a logit model in R and Mathematica

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

-
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. –  Sjoerd C. de Vries Nov 8 '12 at 20:39
With precision and accuracy set to 200 I get {-0.74536032846461989040..., 0.12712830428988635844565...} –  Sjoerd C. de Vries Nov 8 '12 at 20:45
@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. –  whuber Nov 8 '12 at 21:43

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,

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.