Kaplan Meier curve shows longer survival time while multiple regression analysis shows more risk to death

[eric some]

Dear all,
I run a survival analysis to compare Lamivudine Vs Kaletra groups and my Kaplan Meier curve shows that the survival curve of Kaletra is above the Lamivudine curve. HOwever, in my multiple regression analysis using flexible parametric model, I found that Kaletra group was more at risk (3 time) to not survive than Lamivudine group. Is that possible? Aren't they contradictory results?:

[Clyde Schechter]

It's not clear what you mean by "multiple regression analysis using flexible parametric model." You should post the code and the output so we can tell exactly what you're talking about. If the analysis is not one of the commands designed for use with censored survival data, then your results from it are almost guaranteed to be wrong because the censored observations are not represented in the data and estimates of survival based only on observed deaths provides severely biased estimates.

If it is one of the -st- commands, then the interpretation depends on which one. While some produce hazard ratios, some produce accelerated failure time coefficients--which are interpreted differently. So we need to see precisely what you did and what you got.

[Steve Samuels]

Yes, and start with the stset command and some description of the data (e.g. with stdes).

[Carlo Lazzaro]

Eric:
while KM curves are non-parametric (i.e. let the data speak for themselves), semi-parametric (Cox regression) and parametric (-streg- suite) regression models consider different predictors (which are not described in your post) acting, adjusted to each other, on the outcome measure.
As already highlighted by Clyde and Steve, posting what you typed and what Stata gave you back (as per FAQ #12) would increase your chances of getting helpful replies

[eric some]

Thanks so much for your responses. I am sorry that I have not been enough clear. The KM curve was to describe the mean duration of exclusive breastfeeding starting from birth until 6 months post-partum (exit date). The study was multi-country in the context of a RCT with two treatment arm (Kaletra Vs Lamivudine). I present here a curve by treatment arm for one country. The command:

sts graph, by(gp_ttt) tmax(50) tmin(0) xtitle(analysis time in weeks) ///
title(South Africa) risktable(0 6 14 18 22 26, order(1 "Lamivudine" 2 "Kaletra" rowtitle("Lamivudine" "Kaletra" ))) ///
xlabel(0 6 10 14 18 22 26) ylabel(0 0.5 1) legend(order(1 "Lamivudine" 2 "Kaletra") ring(0) position(7) rows(4)) ///
saving(survEBFPLUSRandom2, replace) name(EBFPLUSRandom, replace)

. The output: (see attached file)

To analyse the risk factor for non exclusive breastfeeders, I use the stpm2 parametric command because the PH assumption did not hold for an stcox:

xi: stpm2 i.rand i.agegroup i.school i.marital i.job i.parity  i.deliv i.bfinitime,   ///
 df(4)scale(hazard)tvc(rand school) dftvc(rand:3 school:2) eform

and got the following display

  xi: stpm2 i.rand i.agegroup i.school i.marital i.job i.parity  i.deliv i.bfinitime,   ///
>  df(4)scale(hazard)tvc(rand school) dftvc(rand:3 school:2) eform
i.rand            _Irand_0-1          (naturally coded; _Irand_0 omitted)
i.agegroup        _Iagegroup_1-5      (naturally coded; _Iagegroup_1 omitted)
i.school          _Ischool_0-1        (naturally coded; _Ischool_0 omitted)
i.marital         _Imarital_1-2       (naturally coded; _Imarital_1 omitted)
i.job             _Ijob_0-1           (naturally coded; _Ijob_0 omitted)
i.parity          _Iparity_1-2        (naturally coded; _Iparity_1 omitted)
i.deliv           _Ideliv_1-2         (naturally coded; _Ideliv_1 omitted)
i.bfinitime       _Ibfinitime_1-2     (naturally coded; _Ibfinitime_1 omitted)

Iteration 0:   log likelihood = -721.12689
Iteration 1:   log likelihood = -697.10336
Iteration 2:   log likelihood = -684.48807
Iteration 3:   log likelihood = -671.50954
Iteration 4:   log likelihood = -669.56861
Iteration 5:   log likelihood = -669.54753
Iteration 6:   log likelihood = -669.54752

Log likelihood = -669.54752                       Number of obs   =       1272

--------------------------------------------------------------------------------
               |     exp(b)   Std. Err.      z    P>|z|     [95% Conf. Interval]
---------------+----------------------------------------------------------------
xb             |
      _

Irand_1 |   3.017531   .9309034     3.58   0.000      1.64838    5.523903
  _Iagegroup_3 |   1.440863   .2166636     2.43   0.015     1.073067    1.934721
  _Iagegroup_5 |   .7019129   .1202683    -2.07   0.039     .5016899    .9820445
    _Ischool_1 |   .6115013   .1742286    -1.73   0.084     .3498433    1.068861
   _Imarital_2 |   1.595827   .2114848     3.53   0.000     1.230783    2.069141
       _Ijob_1 |   1.280016   .1586774     1.99   0.046     1.003913    1.632055
    _Iparity_2 |   1.580864   .2347573     3.08   0.002     1.181657    2.114936
     _Ideliv_2 |   .9268712   .1251394    -0.56   0.574     .7113713    1.207654
 _Ibfinitime_2 |   .8782526   .1097678    -1.04   0.299     .6874365    1.122035
         _rcs1 |   7.384381   4.327239     3.41   0.001     2.341591     23.2872
         _rcs2 |   .9220875   .2645931    -0.28   0.777     .5254359    1.618171
         _rcs3 |   1.090468   .0602981     1.57   0.117     .9784646    1.215292
         _rcs4 |   .9605387   .0160671    -2.41   0.016     .9295583    .9925515
    _rcs_rand1 |   .1960135    .146456    -2.18   0.029     .0453201    .8477756
    _rcs_rand2 |   .8797732   .2814744    -0.40   0.689     .4699348    1.647039
    _rcs_rand3 |   .8168775   .0484413    -3.41   0.001     .7272441    .9175584
  _rcs_school1 |   2.157699   1.206741     1.38   0.169     .7210041    6.457199
  _rcs_school2 |   1.053126   .3111209     0.18   0.861     .5902189    1.879089
         _cons |   .0528786   .0190705    -8.15   0.000     .0260794    .1072167
--------------------------------------------------------------------------------

.
end of do-file

rand is the variable for treatment arm (0=lamivudine; 1=Kaletra). As you can seek, the hazard ratio=3.017531 (1.64838; 5.523903) while in the figure, the kaletra curve is above the lamivudine one when I would expect it to be below.
I do not know how to interpret this.

[Carlo Lazzaro]

Eric:
did you perform a -logrank- test after KM?

[Steve Samuels]

Your mistake is in thinking that 3.017531 summarizes the treatment difference. It does not, because the model also includes interactions of treatment with time, generated by the tvc(rand) option. The interaction variables are _rcs_rand1, _rcs_rand2, and _rcs_rand3. All three have HR < 1, two significantly so. Thus the effect of treatment diminishes with time..

The graph is interesting. The two curves are essentially identical after 24 weeks, when ~ 30% of women have stopped exclusive breastfeeding. Almost all stop in the next two weeks. If you plot the adjusted hazard functions , you will see that they explode upward at the end. I don't think that hazard models are necessarily the best way to describe such data. Consider also models of the restricted mean ("search restricted mean").

[eric some]

Thank you all. Dear Carlo: I performed the logrank test after the KM but not specifically for this country. It was for all sites and it was p=0.05.
Dear Steve: Thanks for your very helpful comment and suggestions. I shall look at the restricted mean model