# What kind of logrank test is performed by survdiff function in R?

I was trying to understand how logrank works. For that purpose I used a toy dataset with only a few rows (30 rows). I calculated by hand (with Google Spreadsheet) the logrank statistic (chi squared) and I obtained a value approximately 9.01. Later I loaded the dataset in R, and used survdiff command to do the same task. It returned a value of 7.9 for the same statistic. So I am trying to understand why this could be. I wonder if R uses some sort of weighted logrank, or maybe it performs some other correction, like an estimate value fitted with some sort of curve. Can anyone help me with this issue? Thanks in advance.

Some people suggested this post answers my question: Understanding R survdiff() - where do expected values come from?

No, I actually have the same problem of the original poster, I have slightly different numbers when I calculate the expected values by hand (with Google Spreadsheet or Excell) compaired with the expected values and variance that survdiff outputs. So I don't believe the algorithm is exactly the one explained in this post. Actually, I got closer to the survdiff value when I use the total at risk as a weight for both variance and expected value, but still slightly different. So I guess R is using something different, if someone has access to the inner code of "survdiff" maybe it would help.

My R code:

datos <- read.csv("/Users/usuario/Downloads/survival.csv")

datos_diff <- survdiff(Surv(time, delta) ~ type, data = datos)

datos_diff


R output (the --- are added to fix the spaces):

Type N Observed Expected (O-E)^2/E (O-E)^2/V
1 15 12 17.43 1.69 7.89
2 15 12 6.57 4.48 7.89

$$\chi^2 = 7.9$$ $$\text{df = 1}$$, $$p = 0.005$$

The algorithm I am using by hand:

Expected values for each event time:

$$E_1 = N_1*(M_1+M_2)/(N_1+N_2) \\ E_2 = N_2*(M_1+M_2)/(N_1+N_2)$$

Variance for each event time:

$$V= N_1*N_2*(M_1+M_2)*(N_1+N_2-M_1-M_2)/((N_1+N_2-1)*(N_1+N_2)^2)$$

Where

$$N_1$$ = persons at risk in Group 1 at a given time

$$N_2$$ = persons at risk in Group 2 at a given time

$$M_1$$ = events in Group 1 at a given time

$$M_2$$ = events in Group 1 at a given time

Then I sum the all the differences of observed and expected events (O-E) for all times of one group, and later I square that sum. Finally, I take this result I divide by the sum of variance for all times. This gives me the statistic value, $$9.01$$.

If I weight each of the sum members of numerator with the number of total persons at risk at this time $$(N_1 + N_2)$$ and I weight each of the sum members in the denominator (variances) by the square $$(N_1 + N_2)^2$$, I get a result which is closer R output, $$7.44$$, but still not the same. So I think maybe R is using some other weighting method.

• The help does tell you (with references) what it estimates. Sep 30, 2023 at 13:38
• No, I actually have the same problem of the original poster, I have slightly different numbers when I calculate the expected values by hand (with Google Spreadsheet or Excell) compaired with the expected values and variance that survdiff outputs. So I don't believe the algorithm is exactly the one explained in this post. Actually, I got closer to the survdiff value when I use the total at risk as a weight for both variance and expected value, but still slightly different. So I guess R is using something different, if someone has access to the inner code of "survdiff" maybe it would help. Sep 30, 2023 at 16:08

The R survdiff() function can perform weighted log-rank tests (see its rho argument), but that's not the default. In this case, there seems to be a slight mis-coding of the data, leading to an error in the numbers at risk.

head(datos[order(datos\$time),],7)
#    type time delta
# 26    2    1     1
# 30    2    1     1
# 5     2    2     1
# 22    2    2     1
# 29    1    3     1
# 13    2    4     0
# 23    2    5     1


One of these first 7 cases has a right-censored event time with the last observation at time = 4. That individual thus should be included in the at-risk members of the type=2 group for the event at time = 3. All 4 events prior to time = 3 were in the type=2 group, with no censoring prior to time = 3. With 15 originally in each group, at time = 3 there should be 11 at risk in the type=2 group.

In your linked spreadsheet, the numbers at risk in the type=2 group are in column E. The entry in that column at time = 3 is only 10. That case with the right-censored observation at time = 4 seems to have been omitted, listed as censored at time = 2.

I suspect that the problem came from trying to restrict the rows of the spreadsheet to only include times with events. It looks like your formulas would have worked OK if you had allowed there to be rows at right-censored times as well. At times with no events, the expected events are 0 and the variance is also 0.

I applaud your interest in trying to reproduce such calculations by hand. It's unfortunately very easy to make small mistakes like this. In my experience, relying on Excel makes mistakes even more probable. I'd recommend doing "by-hand" calculations in R instead.

You can often examine the code for an R function by typing the name of the function at the command prompt. If the package isn't loaded or the function isn't exported from the package, you might have to prepend the package name appropriately to see the R code. Sometimes there is code compiled from a different language that is called by the R function, in which case you need to inspect the corresponding source package to see the source for that compiled code.

That's illustrated nicely in this case. Typing survdiff at the command prompt gives the main function, but if you examine the code you see that it calls survdiff.fit(). Typing survdiff.fit at the command prompt doesn't provide a useful result; you can get that R code by typing survival:::survdiff.fit.

Then you will see that survdiff.fit() calls C code in the following line.

xx <- .C(Csurvdiff2, as.integer(n), as.integer(ngroup), as.integer(nstrat),
as.double(rho), as.double(y[ord, 1]), as.integer(y[ord,
2]), as.integer(x[ord]), as.integer(strat2), observed = double(ngroup *
nstrat), expected = double(ngroup * nstrat), var.e = double(ngroup *
ngroup), double(ngroup), double(n))


If you inspect the source package for survival, you can find the C source in the survdiff2.c file. In this case, the C code is so simple that even I could almost understand it.

• Thank you very much for such a detailed answer. You are completely right and this has helped me a lot! This is an excellent explanation addressing all the issues I had. Thank you again! :-) Oct 1, 2023 at 22:56