1
$\begingroup$

I've performed a cox regression in rstudio (version 1.0.136) using the coxph function in the package "OIsurv". I've also performed the same analysis on SPSS using the same dataset, but i keep on getting different results between the two programmes, and i cannot figure out why this is.

My data are not typical cox regression data, they are about learning in bumblebees. The experiment investigated whether parasitism has an effect on bee learning. Each bee underwent a series of 15 trials of the learning experiment.

The variables that i am analysing are:

Did bee learn - binary (0 = no, 1 = yes)

Time to first learned response - the number of trials it took the bee to first show a learned response. I observed bees for 15 trials, therefore for bees that did not learn (did bee learn = 0) , i inputted 15 into this column, as this was the trial where they were last censored.

Treatment - parasite or control

Forager or nest - whether the bee was a forager bee or a nest bee

Average thorax width - continuous variable, used as a proxy for bee size

Age - age of the bee when it underwent the learning trials

This is how my model looks in R: mod7<-coxph(Surv(Time.to.first.learned.response, Did.bee.learn) ~ Forager.or.nest*Treatment + Treatment*Average.thorax.width + Age, data=survmotiv3per, method="breslow")

This is the output from R: enter image description here

In SPSS

The time variable: 'time to first learned response' Status: 'did bee learn' (where 1 indicates the event having occurred) Covariates: these are exactly the same as the model in R

This is the call to COXREG: COXREG Time.to.first.learned.response /STATUS=Did.bee.learn(1) /CONTRAST(Forager.Y.N)=Indicator /CONTRAST(Treatment)=Indicator /METHOD=ENTER Forager.Y.N Treatment Average.TW Age.when.tested Forager.Y.N*Treatment Average.TW*Treatment /PRINT=CI(95) /CRITERIA=PIN(.05) POUT(.10) ITERATE(20).

I did not change the default contrast settings on SPSS.

This is the SPSS output: enter image description here

The differences seem to be occuring in the 'average thorax width' variable. This is coded as 'numeric' in R, and 'scale' in SPSS.

I have read that R and SPSS can produce different cox regression outputs due to their use of different methods of estimation. The default in SPSS is "Breslow", whereas in R it is "Efron". I manually changed the default in R and still the problem of having different results occurs.

Does anyone know why i am seeing different results here?

Improve this question
$\endgroup$
5
  • 2
    $\begingroup$ Hint: 0.40144 + 0.45407 = 0.85551 and 0.856 - 0.454 = 0.402. $\endgroup$
    – mdewey
    Feb 22, 2017 at 15:57
  • 1
    $\begingroup$ For reference, please add to your question the call you made to COXREG in SPSS, as you have done for coxph in R. In particular, whether you made any change to the SPSS default CONTRAST setting. $\endgroup$
    – EdM
    Feb 22, 2017 at 19:07
  • $\begingroup$ @mdewey thanks for your reply, but could you clarify what you mean by this? I do not have lots of experience in statistics. $\endgroup$
    – GaryStats
    Feb 23, 2017 at 8:42
  • 2
    $\begingroup$ I suspect @mdewey is suggesting that your two models differ in whether they regard parasite or control as the base treatment, and perhaps similarly with yes/no for forager v. nest. Another hint at this is that the sign of the treatment-average thorax width interaction coefficient reverses from 0.45407 in the R run to -.454 in the SPSS run $\endgroup$
    – Henry
    Feb 23, 2017 at 9:05
  • $\begingroup$ @Henry indeed. Different software can parameterise the same model in different ways as here. I do not know how SPSS does it but since it is product for which David is paying he can presumably ask their technical support. $\endgroup$
    – mdewey
    Feb 23, 2017 at 9:09

1 Answer 1

Reset to default
4
$\begingroup$

It's not only the "average thorax width" coefficient that differs between the models. The signs of the coefficients for the categorical predictor variables differ, too. These differences arise from the different choices about the reference levels for categorical variables between R and SPSS.

For categorical variables, R by default chooses the lowest level as the reference while SPSS chooses the highest. So for a 2-category predictor like your "Treatment" the signs of the coefficients will be opposite.

When you include an interaction term in a model, you need to be very careful in thinking about what the coefficients mean.

The coefficient for a continuous variable that has an interaction with a categorical variable will be reported as its value when the categorical variable is at its reference level. To get its coefficient for a non-reference level of the categorical variable, you add to that the interaction coefficient.

R and SPSS, as you have used them, differ in their choice of categorical reference levels so they report different values for the "average thorax width" coefficient, and they report opposite signs for the interaction coefficient. As @mdewey implies in comments, R and SPSS thus produce the same coefficients for "average thorax width" at both levels of the "Treatment" variable when you take the difference in categorical reference levels into account.

Added in response to comment:

In your particular application, you are under-powered to test this many coefficients. The usual rule of thumb in Cox or logistic regressions is to have about 15 events per predictor variable being considered. (For this, interaction terms count as predictor variables.) Your 53 events thus would limit you to about 3 predictors, while your model includes 6. Note that your overall model does not reach standard statistical significance (p-value is > 0.05 for the omnibus tests), so you should not be paying much attention to the individual regression coefficients anyway. This model is not significantly different, by standard frequentist criteria, from no model at all.

You might consider a model without interaction terms, particularly as this analysis shows that they are far from significant. That would bring you down to 4 predictors. The p-value calculations wouldn't strictly be correct any more, as this simpler model will have been designed only after you saw the results of this larger model. That distinction, however, is often ignored in practice.

And more data seldom hurt.

Improve this answer
$\endgroup$
5
  • $\begingroup$ I would also check yes/no for forager v. nest $\endgroup$
    – Henry
    Feb 23, 2017 at 10:03
  • $\begingroup$ Thank you very much for your answer. So if i change the reference level of my categorical variables in R, i could recreate the results i am seeing in SPSS? $\endgroup$
    – GaryStats
    Feb 23, 2017 at 11:37
  • $\begingroup$ @EdM Also, is there any particular reason to choose one level or another as the reference level? In this example using different reference levels changes the p-value of the 'average thorax width' variable, which could potentially change the interpretation of the results. How do i chose which reference level is most suitable? Thank you. $\endgroup$
    – GaryStats
    Feb 23, 2017 at 11:41
  • $\begingroup$ Changing the reference level doesn't "change the interpretation of the results"; it changes the particular hypothesis that you are testing. In your R code the displayed p-value indicates whether "average thorax width" is related to outcome when treatment is control; in SPSS, whether it is related to outcome when treatment is parasite. Think carefully about what hypotheses you wish to test. I will edit my answer with a few comments about your particular application. $\endgroup$
    – EdM
    Feb 23, 2017 at 14:30
  • $\begingroup$ @EdM Thanks again for your useful comments. I agree with you about this example having too many predictors. This is just a smaller example dataset, my real dataset is a little larger, but it still probably has too many predictors. With the real dataset the model is also significantly diferent from no model at all. $\endgroup$
    – GaryStats
    Feb 28, 2017 at 10:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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