Kaplan Meier Survival curve results differ between R and SAS?

6
votes

I'm re-running Kaplan-Meier Survival Curves from previously published data, using the exact data set used in the publication (Charpentier et al. 2008 - Inbreeding depression in ring-tailed lemurs (Lemur catta): genetic diversity predicts parasitism, immunocompetence, and survivorship). This publication ran the curves in SAS Version 9, using LIFETEST, to analyze the age at death structured by genetic heterozygosity and sex of the animal (n=64). She reports a Chi square value of 6.31 and a p value of 0.012; however, when I run the curves in R, I get a Chi square value of 0.9 and a p value of 0.821. Can anyone explain this??

R Code used: Age is the time to death, mort is the censorship code, sex is the stratum of gender, and ho2 is the factor delineating the two groups to be compared.

> survdiff(Surv(age, mort1)~ho2+sex,data=mariekmsurv1)
Call:
survdiff(formula = Surv(age, mort1) ~ ho2 + sex, data = mariekmsurv1)

              N Observed Expected (O-E)^2/E (O-E)^2/V
ho2=1, sex=F 18        3     3.23    0.0166    0.0215
ho2=1, sex=M 12        3     2.35    0.1776    0.2140
ho2=2, sex=F 17        5     3.92    0.3004    0.4189
ho2=2, sex=M 17        4     5.50    0.4088    0.6621

Chisq= 0.9  on 3 degrees of freedom, p= 0.821 


> str(mariekmsurv1)
'data.frame':   64 obs. of  6 variables:
 $ id   : Factor w/ 65 levels "","aeschylus",..: 14 31 33 30 47 57 51 39 36 3 ...
 $ sex  : Factor w/ 3 levels "","F","M": 3 2 3 2 2 2 2 2 2 2 ...
 $ mort1: int  0 0 0 0 0 0 0 0 0 0 ...
 $ age  : num  0.12 0.192 0.2 0.23 1.024 ...
 $ sex.1: Factor w/ 3 levels "","F","M": 3 2 3 2 2 2 2 2 2 2 ...
 $ ho2  : int  1 1 1 2 1 1 1 1 1 2 ...
- attr(*, "na.action")=Class 'omit'  Named int [1:141] 65 66 67 68 69 70 71 72 73 74 ...
  .. ..- attr(*, "names")= chr [1:141] "65" "66" "67" "68" ...
2
rsassurvival-analysis
can you put the data in a reproducible form? – Ben Bolker
Does the pub include the SAS code to run PROC LIFETEST? – Joe
When you say chi-square are you referring to a log-rank (Mantel-Haenszel) chi-square test? Was it specified in the paper whether this test was used rather than Peto-Peto, etc? The default in both SAS and R (survdiff with rho=0) is log-rank. – Alex A.
This question appears to be off-topic because it is about statistical methods, not coding. Belongs On CrossValidated.com – IRTFM
Note that the default tie handling method for survival curves is Breslow in SAS and Efron in R. Perhaps try changing your tie handling in R to Breslow to see if it matches the SAS output. – Alex A.

2 Answers

0
votes

Some ideas:

Try running it in SAS -- see if you get the same results as the author. Maybe they didn't send you the exact same dataset they used.

Look into the default values of the relevant SAS PROC and compare to the defaults of the R function you are using.

0
votes

Given the HUGE difference between the Chi-squared (6.81 and 0.9) and P values (0.012 and 0.821) beteween SAS procedure and R procedure for survival analyses; I suspect that you have used wrong variables in the either one of the procedures.

The procedural difference / (data handling difference between SAS and R can cause some very small differences ) .

This is not a software error, this is highly likely to be a human error.