Understanding R survdiff() - where do expected values come from?

Question

I found myself repeating an analysis I don't fully understand.

Example data:

chisq_data <- data.frame(rx = c(rep(1, 7), rep(2,7))
                         ,Vital.Status = rep(1, 14)
                         ,EFS = (c(150, 260, 110, 111, 550, 1, 1, 1,140 , 60, 1, 70, 1, 1))
)
 

I then want to know whether survival of the two groups defined in the rx column as 1 and 2 is different. I use the survdiff() function.

survdiff(Surv(chisq_data$EFS, chisq_data$Vital.Status) ~ rx, 
                 data = chisq_data)

And get this output:

Call:
survdiff(formula = Surv(chisq_data$EFS, chisq_data$Vital.Status) ~ 
    rx, data = chisq_data)

     N Observed Expected (O-E)^2/E (O-E)^2/V
rx=1 7        7     9.72     0.762      4.01
rx=2 7        7     4.28     1.733      4.01

 Chisq= 4  on 1 degrees of freedom, p= 0.05 

I think I understand what the p= 0.04 means. It tests the hypothesis that my groups 1 and 2 come from the same population/distribution. What I don't understand is where do these Expected values come from (and in fact the other values, but it's the Expected column that I would mostly like to understand). What do they mean? Because chi-square test is used, and it works on categorical data, I would also like to know how the 'categories' for this test are created.

In the end, I would like to do power analysis using this dataset as my preliminary dataset. How many more samples do I need to get p<0.01 or p<0.001? I am not sure how to go about it, either, but I feel understanding this first step is necessary to know what to do next.

Edit2:
Following @EdM's answer, I have some follow up/clarifying questions.
@EdM now addresses all what I wrote below in his edited response. What I write below is incorrect because I calculated cumulative number of events rather than number at each time point - see the answer below for a correct data. I am leaving this here as it might be helpful to some to see where my thinking was incorrect

If I understood the answer correctly, the number of events in each of three groups (all, rx1, rx2) is calculated at each time point. It would produce a table as follows:

  EFS prob_dead_all prob_dead_r1 prob_dead_r2
1   1     0.4285714    0.5714286    0.2857143
2  60     0.5000000    0.7142857    0.2857143
3  70     0.5714286    0.8571429    0.2857143
4 110     0.6428571    0.8571429    0.4285714
5 111     0.7142857    0.8571429    0.5714286
6 140     0.7857143    1.0000000    0.5714286
7 150     0.8571429    1.0000000    0.7142857
8 260     0.9285714    1.0000000    0.8571429
9 550     1.0000000    1.0000000    1.0000000

Columns:
EFS - time points

prob_dead_all - proportion of events that have occurred by each timepoint for the whole dataset. For example, at timepoint 1, 6 out of 14 people die, and therefore the number is 6/14 = 0.4285714

prob_dead_r1 - proportion of events that occurred by each timepoint in the rx1 group. For example, at timepoint 1, 4 out of 7 people die, and therefore the number is 4/7 = 0.5714286

prob_dead_r2 - proportion of events that occurred by each timepoint in the rx2 group. For example, at timepoint 1, 2 out of 7 people die, and therefore the number is 2/7 = 0.2857143

Then, the Expected values (shown in the output of the survdiff()) for each of two groups,rx1 and rx2, are calculated at the sum of columns prob_dead_r1 and prob_dead_r2, respectively. However, when I sum up these columns, I see similar, but not identical, results to the survdiff(). Why?

> colSums(chisq_data2 %>% select(-EFS))
prob_dead_all  prob_dead_r2  prob_dead_r1 
     6.428571      5.000000      7.857143 

Answer 1

The default log-rank test performed by survdiff() compares the number of observed events in each group against the number that would have been "expected" if the groups had the same survival functions. There is no problem here with the "observed" values; they are the 7 events specified for each value of rx (no censoring).

The "expected" values are a bit trickier, as they take into account the actual numbers at risk in each group at each event time. Say the the event times are $t_i$. Say that at $t_i$ there are $d_i$ total deaths out of $Y_i$ total individuals at risk. The fraction of all at-risk individuals having an event at time $t_i$ is thus $(d_i/Y_i)$.

Now, say that there are $Y_{ij}$ individuals at risk specifically in group j at time $t_i$. If all groups had the same survival function (the null hypothesis), then the same fraction of those still at risk in each group should have an event at each event time. That is, the number in group $j$ expected to have an event at time $t_i$ should be

$$Y_{ij}\left( \frac{d_i}{Y_i}\right).$$

The "expected" value for group j is the sum, over all event times, of those values expected under the null hypothesis. The "categorization" is thus just by the groups j that you specify in the model, not by the event times (EFS in your data).

With your data, the rx = 1 group has typically longer survival. The $Y_{i1}$ values stay higher longer over time than the $Y_{i2}$ values, so the sum of "expected" deaths in rx = 1 is higher than for rx = 2 or than the actual number of observed deaths for rx = 1.

The comparison of observed against expected values is done with a chi-square test. The last two columns of your output show the associated calculations.

Example with these data

The following data frame summarizes your data by event time.

eventTimes <- unique(chisq_data$EFS) %>% sort()
eventTimeFrame <-data.frame(eventTime=eventTimes,
    atRiskR1= c(7,3,2,1,1,1,0,0,0),
    atRiskR2=c(7,5,5,5,4,3,3,2,1),
    eventNumbers=c(6,1,1,1,1,1,1,1,1))
eventTimeFrame[,"totalAtRisk"] <- eventTimeFrame$atRiskR1 + eventTimeFrame$atRiskR2
eventTimeFrame
#   eventTime atRiskR1 atRiskR2 eventNumbers totalAtRisk
# 1         1        7        7            6          14
# 2        60        3        5            1           8
# 3        70        2        5            1           7
# 4       110        1        5            1           6
# 5       111        1        4            1           5
# 6       140        1        3            1           4
# 7       150        0        3            1           3
# 8       260        0        2            1           2
# 9       550        0        1            1           1

Then do the per-event-time calculations of "expected" events in each group, based on the number still at risk in each group and the overall fraction of cases having an event at that time.

eventTimeFrame[,"expectedR1"] <- eventTimeFrame$atRiskR1 *
    eventTimeFrame$eventNumbers /eventTimeFrame$totalAtRisk
eventTimeFrame[,"expectedR2"] <- eventTimeFrame$atRiskR2 *
    eventTimeFrame$eventNumbers /eventTimeFrame$totalAtRisk
sum(eventTimeFrame[,"expectedR1"])
# [1] 4.277381
sum(eventTimeFrame[,"expectedR2"])
# [1] 9.722619

Those agree with the survdiff() output.

Power estimates

Even if you now understand the function, working with survdiff() is not a very efficient way to get power/sample-size estimates. You need to take into account censoring and the variability in sampling from the population that provided your pilot data. I find the cpower() function in the R Hmisc package to be helpful for that in simple studies with defined accrual patterns and follow-up times. The package's spower() function allows for more complex designs. The CRAN survival task view has additional suggestions.