2
$\begingroup$

Information added: this question is not duplicate, neither my question nor the question in the link had been solved at all, I added more info to see if you can see the problem

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.

Edited to answer:

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

My answer:

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.

Edited to add more info:

My dataset: https://drive.google.com/file/d/1z-Hgl4eAUoNgotrDdsCBnmghVywwbEnO/view?usp=sharing

My spreadsheet: https://docs.google.com/spreadsheets/d/1WvgfCmL_h4DVOlAMt4LB0Cdz4nUsNrUy7CW8DmleGEU/edit#gid=0

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.

$\endgroup$
2
  • $\begingroup$ The help does tell you (with references) what it estimates. $\endgroup$
    – mdewey
    Commented Sep 30, 2023 at 13:38
  • $\begingroup$ 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. $\endgroup$
    – Estrella
    Commented Sep 30, 2023 at 16:08

1 Answer 1

4
$\begingroup$

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.

Here's the head of your data frame, ordered by event times.

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.

Access to the inner code

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.

$\endgroup$
1
  • 1
    $\begingroup$ 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! :-) $\endgroup$
    – Estrella
    Commented Oct 1, 2023 at 22:56

Your Answer

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

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