Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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?

share|improve this question
1  
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

1 Answer 1

up vote 10 down vote accepted

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.

share|improve this answer
    
Thank you. It was one of those mistakes that seems obvious in hindsight. –  Mr Alpha Nov 9 '12 at 8:31

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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