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

Question

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?

Answer 1

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

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.