3
$\begingroup$

My colleague and I are working with complex survey data (i.e., with weights). He is doing the analysis in SPSS, and I am trying to double-check his work, but I use R.

We have a binary explanatory variable and a binary outcome, and we'd like to estimate the association (odds ratio, OR), using logistic regression. Eventually we will do multivariable analysis, but we are first trouble-shooting univariate results (1 binary explanatory variable).

We've noticed that the point estimates (coefficients and ORs, i.e., exp(beta)) from the two softwares are identical, but standard errors (SE) and confidence intervals differ.

For example, using the same variable definitions, he gets the following logistic regression OR and CI, using SPSS: 

1.885 (1.611 – 2.206)

B=0.634, SE=0.080

[full SPSS outputs included at bottom of this post]

And I get, using R (svyglm): 

1.885457 (1.407309 - 2.52606)

B=0.63417, SE=0.14923

[full R outputs included at bottom of this post]

In other words, the SEs and CIs in SPSS are larger than those estimate in R.

Any ideas why this is the case? We're wondering if we've made an error or if the softwares actually calculate the SEs in different ways.

Thank you for any help!

SPSS syntax:
WEIGHT BY New_Weight.
DATASET ACTIVATE DataSet1.

SAVE OUTFILE=
    '/Users/jd/datafilec3.sav'
  /COMPRESSED.
LOGISTIC REGRESSION VARIABLES SRMH
  /METHOD=ENTER sexual_orientation
  /CONTRAST (sexual_orientation)=Indicator(1)
  /PRINT=CI(95)
  /CRITERIA=PIN(0.05) POUT(0.10) ITERATE(20) CUT(0.5).

SPSS outputs: enter image description here

R syntax:

> summary(svyglm(srmh.r ~ as.factor(so.r), design=joshunsvy, family=binomial))

R outputs:

Call:
svyglm(formula = srmh.r ~ as.factor(so.r), design = joshunsvy, 
    family = binomial)

Survey design:
svydesign(id = ~1, weights = ~WTS_M, data = joshun)

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)      -2.75422    0.02257 -122.06  < 2e-16
as.factor(so.r)1  0.63417    0.14923    4.25 2.14e-05
as.factor(so.r)2  1.35152    0.08300   16.28  < 2e-16

(Dispersion parameter for binomial family taken to be 0.9997891)

Number of Fisher Scoring iterations: 5
| cite | improve this question | |
$\endgroup$
  • 2
    $\begingroup$ It might be easier to interpret the results with the complete output and functions used. Now we just see that the size of the confidence intervals differs by a factor 1.86 (and probably it is the same for the standard errors although that can not be seen here) but that can be anything. $\endgroup$ – Sextus Empiricus May 27 at 15:02
  • $\begingroup$ Thank you--I have added the complete outputs to the post, along with the function (for R). $\endgroup$ – Travis May 27 at 22:48
  • $\begingroup$ It would also help to show the SPSS syntax (even if you did it via menus, I think the code appears somewhere). That would help with working out if the two programs are trying to do the same thing. Also, your predictor variable can't actually be binary, or R would be giving you one coefficient rather than two -- does it have some zero or negative values? $\endgroup$ – Thomas Lumley May 27 at 23:06
  • $\begingroup$ Thank you Thomas. I have added the SPSS syntax, above. $\endgroup$ – Travis May 27 at 23:20
  • $\begingroup$ OK. That's your problem. You're using SPSS LOGISTIC procedure, which doesn't handle sampling weights correctly. The CSLOGISTIC procedure in the Complex Samples module does handle sampling weights. I'm a bit surprised the SPSS standard errors are larger rather than smaller, but that will depend on how the weights are scaled. $\endgroup$ – Thomas Lumley May 27 at 23:26
5
$\begingroup$

SPSS LOGISTIC does not handle sampling weights correctly for computing standard errors.

If you have weights $w_i$ for each observation, SPSS will work out the loglikelihood contribution $\ell_i(\beta)$ for each observation, and maximise the weighted sum $\hat\ell(\beta) = \sum_i w_i\ell_i(\beta)$. So will R. The point estimates will agree exactly.

SPSS, however, will calculate the variance matrix for $\hat\beta$ by treating $\hat\ell(\beta)$ as a real loglikelihood. The estimated variance matrix will be the inverse of the second derivative of $\hat\ell(\beta)$. This would be correct if the $w_i$ were precision weights, (so a weight of 10 means 10 times lower variance) or frequency weights (so a weight of 10 means there were 10 identical observations that are stored this way to save space). In both these settings, a higher weight means the observation intrinsically carries more information and so reduces the standard errors by more.

It is not correct for sampling weights, where a weight of 10 means the observation represents 10 observations in the population. Here, an observation with a high weight doesn't reduce the standard errors more -- perhaps even the reverse, since it will have high influence.

SPSS does have a "Complex Samples" add-on that does correct standard error estimation with sampling weights, giving a CSLOGISTIC procedure. You get the same phenomenon with SAS, where PROC LOGISTIC does not handle survey weights correctly, but PROC SURVEYLOGISTIC does. And, I suppose, with R, where glm doesn't but svyglm does.

| cite | improve this answer | |
$\endgroup$

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.