{"id":1196,"date":"2015-08-27T10:37:04","date_gmt":"2015-08-27T09:37:04","guid":{"rendered":"http:\/\/pcool.dyndns.org:8080\/statsbook\/?page_id=1196"},"modified":"2025-06-26T09:14:05","modified_gmt":"2025-06-26T08:14:05","slug":"answers-survival-analysis","status":"publish","type":"page","link":"https:\/\/pcool.dyndns.org\/index.php\/answers-survival-analysis\/","title":{"rendered":"Answers Survival Analysis"},"content":{"rendered":"\n

The table below shows 20 patients who have been diagnosed with cancer. The first column shows the follow up (in years) of the patients who were alive at review. In the second column, the follow up till time of death is indicated. The third column shows the time to last review in the patients who were lost to follow up:<\/em><\/p>\n\n\n\n

The data is also available in Q9.rda <\/a>and for convenience the table is shown below:<\/p>\n\n\n\n

\"\"<\/figure>\n\n\n\n

1. Calculate the 5-year survival in the best-case scenario,<\/em> using life table analysis.<\/strong><\/p>\n\n\n\n

The best-case scenario<\/em><\/strong> is as follows:<\/p>\n\n\n\n

\"\"<\/figure>\n\n\n\n

So, the 5-year survival in the best-case scenario<\/em> as estimated with life table analysis is 41%.<\/p>\n\n\n\n

2. Calculate the 5-year survival in the worst-case scenario,<\/em> using life table analysis.<\/strong><\/p>\n\n\n\n

The worst-case scenario <\/em><\/strong>life table is as follows:<\/p>\n\n\n\n

\"\"<\/figure>\n\n\n\n

So, the 5-year survival in the worst-case scenario<\/em> as estimated with life table analysis is 20%.<\/p>\n\n\n\n

3. Show the best-case scenario <\/em>and worst-case scenario <\/em>survival curves in one graph, using life table analysis.<\/strong><\/p>\n\n\n\n

Using the tables constructed in Q1 and Q2, survival curves can be drawn:<\/p>\n\n\n\n

\"\"<\/figure>\n\n\n\n

The best-case scenario<\/em> <\/strong>is shown in black and the worst-case scenario<\/em> <\/strong>is shown in grey.<\/p>\n\n\n\n

4. Perform the Kaplan-Meier survival analysis in the best-case scenario<\/em>.<\/strong><\/p>\n\n\n\n

In the best-case scenario<\/em><\/strong>, the patients lost to follow up count as a success. The Kaplan-Meier table is as follows:<\/p>\n\n\n\n

\"\"<\/figure>\n\n\n\n

In R:<\/p>\n\n\n\n

Once the data frame is loaded, the data should be visible with the variables ‘Number’, ‘FU’ and ‘Outcome’. <\/p>\n\n\n\n

load(\"\/path\/to_file\/Q9 (1).rda\")\nQ9<\/mark>\n   Number  FU Outcome\n1       1 0.2   alive\n2       2 0.4   alive\n3       3 1.5   alive\n4       4 2.2   alive\n5       5 2.5   alive\n6       6 3.1   alive\n7       7 3.5   alive\n8       8 4.1   alive\n9       9 0.6    dead\n10     10 1.2    dead\n11     11 1.4    dead\n12     12 1.9    dead\n13     13 2.1    dead\n14     14 2.5    dead\n15     15 3.8    dead\n16     16 0.8    lost\n17     17 1.4    lost\n18     18 1.8    lost\n19     19 2.1    lost\n20     20 3.6    lost<\/mark><\/em><\/code><\/pre>\n\n\n\n
It is of course possible to create a new logical variable that is needed for survival analysis (called ‘Best’) and go through the data manually; setting the value to FALSE (0) if the ‘Outcome’ variable is ‘alive’ or ‘lost’ and TRUE (1) if the ‘Outcome’ variable is ‘dead’. However, this would become tedious for larger data sets. It would be easier to recode the variable ‘Outcome’ in a new variable ‘Best’.<\/p>\n\n\n\n
Q9$Best <- ifelse(Q9$Outcome == 'dead', 1, 0)\nQ9<\/mark>\n   Number  FU Outcome Best\n1       1 0.2   alive    0\n2       2 0.4   alive    0\n3       3 1.5   alive    0\n4       4 2.2   alive    0\n5       5 2.5   alive    0\n6       6 3.1   alive    0\n7       7 3.5   alive    0\n8       8 4.1   alive    0\n9       9 0.6    dead    1\n10     10 1.2    dead    1\n11     11 1.4    dead    1\n12     12 1.9    dead    1\n13     13 2.1    dead    1\n14     14 2.5    dead    1\n15     15 3.8    dead    1\n16     16 0.8    lost    0\n17     17 1.4    lost    0\n18     18 1.8    lost    0\n19     19 2.1    lost    0\n20     20 3.6    lost    0<\/mark><\/em>\nstr(Q9)<\/mark>\n'data.frame':\t20 obs. of  4 variables:\n $ Number : num  1 2 3 4 5 6 7 8 9 10 ...\n $ FU     : num  0.2 0.4 1.5 2.2 2.5 3.1 3.5 4.1 0.6 1.2 ...\n $ Outcome: Factor w\/ 3 levels \"alive\",\"dead\",..: 1 1 1 1 1 1 1 1 2 2 ...\n $ Best   : num  0 0 0 0 0 0 0 0 1 1 ...<\/mark><\/em><\/code><\/pre>\n\n\n\n
If the variable Outcome (Q9$Outcome) is equal to (==) ‘dead’, then the Best variable (Q9$Best) = 1, else it is 0.<\/p>\n\n\n\n
As can be seen in the structure of the data frame (str), the newly created variable ‘Best’ is numerical and takes values 0 and 1.<\/p>\n\n\n\n
Now perform the survival analysis with the survival package1<\/a><\/sup>. Make sure the package survival is installed<\/a>.<\/p>\n\n\n\n
my_survival <- Surv(Q9$FU, Q9$Best)\nmy_survival\n<\/mark> [1] 0.2+ 0.4+ 1.5+ 2.2+ 2.5+ 3.1+ 3.5+ 4.1+ 0.6  1.2  1.4  1.9  2.1  2.5  3.8  0.8+ 1.4+ 1.8+\n[19] 2.1+ 3.6+<\/mark>\nsummary(my_survival)<\/mark>\n      time           status    \n Min.   :0.200   Min.   :0.00  \n 1st Qu.:1.350   1st Qu.:0.00  \n Median :2.000   Median :0.00  \n Mean   :2.035   Mean   :0.35  \n 3rd Qu.:2.650   3rd Qu.:1.00  \n Max.   :4.100   Max.   :1.00  <\/mark><\/em><\/code><\/pre>\n\n\n\n
And calculate the survival curve:<\/p>\n\n\n\n
my_curve <- survfit(my_survival ~ 1, conf.int = 0.95, data=Q9)\nsummary(my_curve)<\/mark>\nCall: survfit(formula = my_survival ~ 1, data = Q9, conf.int = 0.95)\n\n time n.risk n.event survival std.err lower 95% CI upper 95% CI\n  0.6     18       1    0.944  0.0540       0.8443        1.000\n  1.2     16       1    0.885  0.0763       0.7477        1.000\n  1.4     15       1    0.826  0.0913       0.6655        1.000\n  1.9     11       1    0.751  0.1096       0.5644        1.000\n  2.1     10       1    0.676  0.1217       0.4751        0.962\n  2.5      7       1    0.580  0.1374       0.3642        0.922\n  3.8      2       1    0.290  0.2161       0.0672        1.000<\/mark><\/em><\/code><\/pre>\n\n\n\n
Plot the curve:<\/p>\n\n\n\n
par(mar=c(2,2,2,2)) # to increase plot margins in RStudio<\/em>\nplot(my_curve, col='black', main='Kaplan Meier Curve',  \n  xlab='Follow Up Duration', ylab='Cumulative Survival Probability')<\/em><\/mark><\/code><\/pre>\n\n\n\n
\"\"<\/figure>\n\n\n\n
Please note the survival estimate is similar to the table above, but the last two values are slightly different due to rounding error.<\/p>\n\n\n\n
5. What is the survival in Q4 at 4.1 years?<\/strong><\/p>\n\n\n\n
The survival at 4.1 years in the best-case scenario<\/em> as estimated with Kaplan-Meier analysis as calculated in the table in Q4 above is 28%.<\/p>\n\n\n\n
In R:<\/p>\n\n\n\n
Using the output from R in Q4 above: 29%.<\/span><\/p>\n\n\n\n
Or alternatively in the console on the recoded variable ‘Best’ following on from above:<\/p>\n\n\n\n
summary(my_curve,time=4.1)<\/mark>\nCall: survfit(formula = my_survival ~ 1, conf.int = 0.95)\n\n time n.risk n.event survival std.err lower 95% CI upper 95% CI\n  4.1      1       7     0.29 <\/strong>  0.216       0.0672            1<\/mark><\/em><\/code><\/pre>\n\n\n\n
Therefore, the survival at 4.1 years is 29%.<\/p>\n\n\n\n
6. What is the median survival in Q4?<\/strong><\/p>\n\n\n\n
For the answer to this question we need to plot the survival curve:<\/p>\n\n\n\n
\"\"<\/figure>\n\n\n\n
<\/a>It can be seen from the graph that the median survival \u2248 3.75 years.<\/p>\n\n\n\n
In R:<\/p>\n\n\n\n
Alternatively, using the console method (as above):<\/span><\/span><\/p>\n\n\n\n
my_curve<\/mark>\nCall: survfit(formula = my_survival ~ 1, data = Q9, conf.int = 0.95)\n\n      n events median 0.95LCL 0.95UCL\n[1,] 20      7    3.8<\/strong>     2.1      NA<\/mark><\/em><\/code><\/pre>\n\n\n\n
Therefore, the median survival is 3.8 years.<\/p>\n\n\n\n
7. Perform the Kaplan-Meier survival analysis in the worst-case scenario<\/em>.<\/strong><\/p>\n\n\n\n
In the worst-case scenario<\/em><\/strong>, the patients lost to follow up count as a failure. <\/p>\n\n\n\n
\"\"<\/figure>\n\n\n\n
Similarly as in Q4, the ‘Outcome’ variable can be recoded into a new variable ‘Worst’; when the ‘Outcome’ is ‘alive’ ‘Worst’ is coded as 0 and when the ‘Outcome’ is ‘dead’ or ‘lost’ ‘Worst’ is coded as 1:<\/span><\/span><\/p>\n\n\n\n
Q9$Worst <- ifelse(Q9$Outcome == 'alive', 0, 1)<\/em>\nQ9<\/em><\/mark>\n   Number  FU Outcome Best Worst\n1       1 0.2   alive    0     0\n2       2 0.4   alive    0     0\n3       3 1.5   alive    0     0\n4       4 2.2   alive    0     0\n5       5 2.5   alive    0     0\n6       6 3.1   alive    0     0\n7       7 3.5   alive    0     0\n8       8 4.1   alive    0     0\n9       9 0.6    dead    1     1\n10     10 1.2    dead    1     1\n11     11 1.4    dead    1     1\n12     12 1.9    dead    1     1\n13     13 2.1    dead    1     1\n14     14 2.5    dead    1     1\n15     15 3.8    dead    1     1\n16     16 0.8    lost    0     1\n17     17 1.4    lost    0     1\n18     18 1.8    lost    0     1\n19     19 2.1    lost    0     1\n20     20 3.6    lost    0     1<\/mark><\/em>\nstr(Q9)<\/mark>\n'data.frame':\t20 obs. of  5 variables:\n $ Number : num  1 2 3 4 5 6 7 8 9 10 ...\n $ FU     : num  0.2 0.4 1.5 2.2 2.5 3.1 3.5 4.1 0.6 1.2 ...\n $ Outcome: Factor w\/ 3 levels \"alive\",\"dead\",..: 1 1 1 1 1 1 1 1 2 2 ...\n $ Best   : num  0 0 0 0 0 0 0 0 1 1 ...\n $ Worst  : num  0 0 0 0 0 0 0 0 1 1 ...<\/mark><\/em><\/code><\/pre>\n\n\n\n
Next perform the survival analysis as described above.<\/p>\n\n\n\n
my_survival <- Surv(Q9$FU, Q9$Worst)\nmy_survival\n<\/em><\/mark> [1] 0.2+ 0.4+ 1.5+ 2.2+ 2.5+ 3.1+ 3.5+ 4.1+ 0.6  1.2  1.4  1.9 \n[13] 2.1  2.5  3.8  0.8  1.4  1.8  2.1  3.6 \n<\/mark><\/em>summary(my_survival)\n<\/em><\/mark>      time           status   \n Min.   :0.200   Min.   :0.0  \n 1st Qu.:1.350   1st Qu.:0.0  \n Median :2.000   Median :1.0  \n Mean   :2.035   Mean   :0.6  \n 3rd Qu.:2.650   3rd Qu.:1.0  \n Max.   :4.100   Max.   :1.0  <\/mark><\/em>\nmy_curve <- survfit(my_survival ~ 1, data= Q9, conf.int=0.95)\nsummary(my_curve)\n<\/em><\/mark>Call: survfit(formula = my_survival ~ 1, data = Q9, conf.int = 0.95)\n\n time n.risk n.event survival std.err lower 95% CI upper 95% CI\n  0.6     18       1    0.944  0.0540       0.8443        1.000\n  0.8     17       1    0.889  0.0741       0.7549        1.000\n  1.2     16       1    0.833  0.0878       0.6778        1.000\n  1.4     15       2    0.722  0.1056       0.5423        0.962\n  1.8     12       1    0.662  0.1126       0.4743        0.924\n  1.9     11       1    0.602  0.1174       0.4107        0.882\n  2.1     10       2    0.481  0.1209       0.2944        0.788\n  2.5      7       1    0.413  0.1216       0.2316        0.735\n  3.6      3       1    0.275  0.1385       0.1026        0.738\n  3.8      2       1    0.138  0.1194       0.0251        0.754<\/mark><\/em>\npar(mar=c(2,2,2,2)) # to increase plot margins in RStudio\nplot(my_curve, col='black', main='Kaplan Meier Curve', <\/em>\n  xlab='Follow Up Duration', ylab='Cumulative Survival Probability')<\/em><\/mark>\n<\/em><\/code><\/pre>\n\n\n\n
\"\"<\/figure>\n\n\n\n
8. What is the survival in Q7 at 4.1 years?<\/strong><\/p>\n\n\n\n
Using the table in Q7: 13.8%.<\/span><\/p>\n\n\n\n
Or alternatively in the console on the recoded variable ‘Worst’:<\/p>\n\n\n\n
summary(my_curve,time=4.1)<\/mark>\nCall: survfit(formula = my_survival ~ 1, data = Q9, conf.int = 0.95)\n\n time n.risk n.event survival std.err lower 95% CI upper 95% CI\n  4.1      1      12    0.138<\/strong>   0.119       0.0251        0.754<\/mark><\/em><\/code><\/pre>\n\n\n\n
Therefore, the survival at 4.1 years is 13.8%.<\/p>\n\n\n\n
9. What is the median survival in Q7?<\/strong><\/p>\n\n\n\n
For the answer to this question we need to plot the survival curve:<\/p>\n\n\n\n
\"\"<\/figure>\n\n\n\n
It can be estimated from the graph that the median survival \u2248 2.25 years.<\/p>\n\n\n\n
In R:<\/p>\n\n\n\n
my_curve<\/mark>\nCall: survfit(formula = my_survival ~ 1, data = Q9, conf.int = 0.95)\n\n      n events median 0.95LCL 0.95UCL\n[1,] 20     12    2.1<\/strong>     1.8      NA<\/mark><\/code><\/pre>\n\n\n\n
Therefore, the median survival is 2.1 years.<\/p>\n\n\n\n
The outcome data of 20 patients who had a total ankle replacement are listed in  <\/em>survivalankle.rda.<\/a><\/p>\n\n\n\n
10. Calculate the 10-year Kaplan-Meier survival of the prosthesis using revision for aseptic loosening<\/em> as \u2018hard end point\u2019.<\/strong><\/p>\n\n\n\n
Only patients who had a \u2018Revision\u2019 for aseptic loosening<\/em><\/strong> are counted as failures. The Kaplan-Meier table is as follows:<\/p>\n\n\n\n
\"\"<\/figure>\n\n\n\n
<\/a>The 10-year Kaplan-Meier survival for aseptic loosening is therefore 80%:<\/p>\n\n\n\n
\"\"<\/figure>\n\n\n\n
In R:<\/p>\n\n\n\n
First recode the ‘Outcome’ variable to a variable called ‘Censor’ as described above (Q4).<\/p>\n\n\n\n
load(\"\/path\/to_file\/survivalankle.rda\")\nsurvivalankle<\/mark>\n   Number   FU                     Outcome\n1       1  1.1                        Well\n2       2  1.2         Revised (infection)\n3       3  1.9                        Well\n4       4  2.1                        Well\n5       5  2.9                        Well\n6       6  3.1                        Well\n7       7  3.9                        Well\n8       8  5.0                        Well\n9       9  5.3 Revised (aseptic loosening)\n10     10  5.5                        Well\n11     11  6.1                  Lost to FU\n12     12  6.3                        Well\n13     13  6.9 Revised (aseptic loosening)\n14     14  7.1                        Well\n15     15  7.3                  Lost to FU\n16     16  7.4                  Lost to FU\n17     17  7.5                  Lost to FU\n18     18  8.9                        Well\n19     19  9.5                  Lost to FU\n20     20 10.2                        Well\n<\/mark>survivalankle$Censor <- ifelse(survivalankle$Outcome == 'Revised (aseptic loosening)', 1, 0)<\/mark>\n<\/mark>survivalankle<\/mark>\n   Number   FU                     Outcome Censor\n1       1  1.1                        Well      0\n2       2  1.2         Revised (infection)      0\n3       3  1.9                        Well      0\n4       4  2.1                        Well      0\n5       5  2.9                        Well      0\n6       6  3.1                        Well      0\n7       7  3.9                        Well      0\n8       8  5.0                        Well      0\n9       9  5.3 Revised (aseptic loosening)      1\n10     10  5.5                        Well      0\n11     11  6.1                  Lost to FU      0\n12     12  6.3                        Well      0\n13     13  6.9 Revised (aseptic loosening)      1\n14     14  7.1                        Well      0\n15     15  7.3                  Lost to FU      0\n16     16  7.4                  Lost to FU      0\n17     17  7.5                  Lost to FU      0\n18     18  8.9                        Well      0\n19     19  9.5                  Lost to FU      0\n20     20 10.2                        Well      0<\/mark><\/em>\nsurvival_ankle <- Surv(survivalankle$FU, survivalankle$Censor)<\/mark>\n<\/mark>summary(survival_ankle)<\/mark>\n      time            status   \n Min.   : 1.100   Min.   :0.0  \n 1st Qu.: 3.050   1st Qu.:0.0  \n Median : 5.800   Median :0.0  \n Mean   : 5.460   Mean   :0.1  \n 3rd Qu.: 7.325   3rd Qu.:0.0  \n Max.   :10.200   Max.   :1.0  \n<\/mark>ankle_curve <- survfit(survival_ankle ~ 1, data=survivalankle, conf.int=0.95)<\/mark>\n<\/mark>summary(ankle_curve)<\/mark>\nCall: survfit(formula = survival_ankle ~ 1, data = survivalankle, conf.int = 0.95)\n\n time n.risk n.event survival std.err lower 95% CI upper 95% CI\n  5.3     12       1    0.917  0.0798        0.773            1\n  6.9      8       1    0.802  0.1279        0.587            1<\/mark><\/em><\/code><\/pre>\n\n\n\n
Plot the curve:<\/p>\n\n\n\n
par(mar=c(2,2,2,2)) # to increase plot margins in RStudio\nplot(ankle_curve, col='black', main='Kaplan Meier Curve', <\/em>\n  xlab='Follow Up Duration', ylab='Cumulative Survival Probability')<\/em><\/mark><\/code><\/pre>\n\n\n\n
\"\"<\/figure>\n\n\n\n
Therefore, the 10 year survival is 80%.<\/p>\n\n\n\n
In the console:<\/p>\n\n\n\n
summary(ankle_curve,time=10)<\/mark>\nCall: survfit(formula = survival_ankle ~ 1, data = survivalankle, conf.int = 0.95)\n\n time n.risk n.event survival std.err lower 95% CI upper 95% CI\n   10      1       2    0.802<\/strong>   0.128        0.587            1<\/mark><\/em><\/code><\/pre>\n\n\n\n
11. Calculate the 10-year Kaplan-Meier survival of the prosthesis using revision as \u2018hard end point\u2019.<\/strong><\/p>\n\n\n\n
All patients who had a \u2018Revision\u2019<\/strong> are counted as failures. The Kaplan-Meier table is as follows:<\/p>\n\n\n\n
\"\"<\/figure>\n\n\n\n
The 10-year Kaplan-Meier survival for revision is therefore 76%:<\/p>\n\n\n\n
\"\"<\/figure>\n\n\n\n
In R:<\/p>\n\n\n\n
Perform recoding similar as to above into a new variable ‘Censor2’, but now counting all revisions as failures and perform the analysis.<\/p>\n\n\n\n
survivalankle$Censor2 <- ifelse(survivalankle$Outcome == 'Revised (aseptic loosening)', 1, 0)\nsurvivalankle$Censor2 <- ifelse(survivalankle$Outcome == 'Revised (infection)', 1, survivalankle$Censor2)\nsurvivalankle<\/mark>\n   Number   FU                     Outcome Censor Censor2\n1       1  1.1                        Well      0       0\n2       2  1.2         Revised (infection)      0       1\n3       3  1.9                        Well      0       0\n4       4  2.1                        Well      0       0\n5       5  2.9                        Well      0       0\n6       6  3.1                        Well      0       0\n7       7  3.9                        Well      0       0\n8       8  5.0                        Well      0       0\n9       9  5.3 Revised (aseptic loosening)      1       1\n10     10  5.5                        Well      0       0\n11     11  6.1                  Lost to FU      0       0\n12     12  6.3                        Well      0       0\n13     13  6.9 Revised (aseptic loosening)      1       1\n14     14  7.1                        Well      0       0\n15     15  7.3                  Lost to FU      0       0\n16     16  7.4                  Lost to FU      0       0\n17     17  7.5                  Lost to FU      0       0\n18     18  8.9                        Well      0       0\n19     19  9.5                  Lost to FU      0       0\n20     20 10.2                        Well      0       0<\/mark><\/em><\/code><\/pre>\n\n\n\n
Survival analysis:<\/p>\n\n\n\n
survival_ankle <- Surv(survivalankle$FU, survivalankle$Censor2)\nsurvival_ankle<\/mark>\n [1]  1.1+  1.2   1.9+  2.1+  2.9+  3.1+  3.9+  5.0+  5.3   5.5+\n[11]  6.1+  6.3+  6.9   7.1+  7.3+  7.4+  7.5+  8.9+  9.5+ 10.2+\n<\/mark>summary(survival_ankle)<\/mark>\n      time            status    \n Min.   : 1.100   Min.   :0.00  \n 1st Qu.: 3.050   1st Qu.:0.00  \n Median : 5.800   Median :0.00  \n Mean   : 5.460   Mean   :0.15  \n 3rd Qu.: 7.325   3rd Qu.:0.00  \n Max.   :10.200   Max.   :1.00  \n<\/mark>ankle_curve <- survfit(survival_ankle ~ 1, data=survivalankle, conf.int=0.95)\nsummary(ankle_curve)<\/mark>\nCall: survfit(formula = survival_ankle ~ 1, data = survivalankle, conf.int = 0.95)\n\n time n.risk n.event survival std.err lower 95% CI upper 95% CI\n  1.2     19       1    0.947  0.0512        0.852            1\n  5.3     12       1    0.868  0.0890        0.710            1\n  6.9      8       1    0.760  0.1280        0.546            1<\/mark><\/em><\/code><\/pre>\n\n\n\n
Plot the curve:<\/p>\n\n\n\n
par(mar=c(2,2,2,2)) # to increase plot margins in RStudio\nplot(ankle_curve, col='black', main='Kaplan Meier Curve', <\/em>\n  xlab='Follow Up Duration', ylab='Cumulative Survival Probability')<\/em><\/mark><\/code><\/pre>\n\n\n\n
\"\"<\/figure>\n\n\n\n
<\/a>Therefore, the 10 year survival is 76%.<\/p>\n\n\n\n
12. Calculate the 10-year worst-case scenario<\/em> Kaplan-Meier survival of the prosthesis using revision as \u2018hard end point\u2019.<\/strong><\/p>\n\n\n\n
In the worst-case scenario<\/em><\/strong>, all patients who are lost to follow up are also counted as a failure:<\/p>\n\n\n\n
\"\"<\/figure>\n\n\n\n
<\/a>The worst-case scenario <\/em>10-year Kaplan-Meier survival for revision is therefore 17%:<\/p>\n\n\n\n
\"\"<\/figure>\n\n\n\n
In R:<\/p>\n\n\n\n
Again, perform recoding similar to above into a new variable called ‘Censor3’, now counting all revisions and all patients lost to follow up as a failure.<\/p>\n\n\n\n
survivalankle$Censor3 <- ifelse(survivalankle$Outcome == 'Revised (aseptic loosening)', 1, 0)\nsurvivalankle$Censor3 <- ifelse(survivalankle$Outcome == 'Revised (infection)', 1, survivalankle$Censor3)\nsurvivalankle$Censor3 <- ifelse(survivalankle$Outcome == 'Lost to FU', 1, survivalankle$Censor3)\nsurvivalankle<\/mark>\n   Number   FU                     Outcome Censor Censor2 Censor3\n1       1  1.1                        Well      0       0       0\n2       2  1.2         Revised (infection)      0       1       1\n3       3  1.9                        Well      0       0       0\n4       4  2.1                        Well      0       0       0\n5       5  2.9                        Well      0       0       0\n6       6  3.1                        Well      0       0       0\n7       7  3.9                        Well      0       0       0\n8       8  5.0                        Well      0       0       0\n9       9  5.3 Revised (aseptic loosening)      1       1       1\n10     10  5.5                        Well      0       0       0\n11     11  6.1                  Lost to FU      0       0       1\n12     12  6.3                        Well      0       0       0\n13     13  6.9 Revised (aseptic loosening)      1       1       1\n14     14  7.1                        Well      0       0       0\n15     15  7.3                  Lost to FU      0       0       1\n16     16  7.4                  Lost to FU      0       0       1\n17     17  7.5                  Lost to FU      0       0       1\n18     18  8.9                        Well      0       0       0\n19     19  9.5                  Lost to FU      0       0       1\n20     20 10.2                        Well      0       0       0<\/mark><\/em><\/code><\/pre>\n\n\n\n
Perform survival analysis:<\/p>\n\n\n\n
survival_ankle = Surv(survivalankle$FU, survivalankle$Censor3)\nsurvival_ankle<\/mark>\n [1]  1.1+  1.2   1.9+  2.1+  2.9+  3.1+  3.9+  5.0+  5.3   5.5+\n[11]  6.1   6.3+  6.9   7.1+  7.3   7.4   7.5   8.9+  9.5  10.2+\n<\/mark>summary(survival_ankle)<\/mark>\n      time            status   \n Min.   : 1.100   Min.   :0.0  \n 1st Qu.: 3.050   1st Qu.:0.0  \n Median : 5.800   Median :0.0  \n Mean   : 5.460   Mean   :0.4  \n 3rd Qu.: 7.325   3rd Qu.:1.0  \n Max.   :10.200   Max.   :1.0  \n<\/mark>ankle_curve <- survfit(survival_ankle ~ 1, data=survivalankle, conf.int=0.95)\nsummary(ankle_curve)<\/mark>\nCall: survfit(formula = survival_ankle ~ 1, data = survivalankle, conf.int = 0.95)\n\n time n.risk n.event survival std.err lower 95% CI upper 95% CI\n  1.2     19       1    0.947  0.0512       0.8521        1.000\n  5.3     12       1    0.868  0.0890       0.7104        1.000\n  6.1     10       1    0.782  0.1149       0.5859        1.000\n  6.9      8       1    0.684  0.1359       0.4633        1.000\n  7.3      6       1    0.570  0.1538       0.3358        0.967\n  7.4      5       1    0.456  0.1598       0.2294        0.906\n  7.5      4       1    0.342  0.1552       0.1404        0.833\n  9.5      2       1    0.171  0.1437       0.0329        0.888<\/mark><\/em><\/code><\/pre>\n\n\n\n
Plot the survival curve:<\/p>\n\n\n\n
par(mar=c(2,2,2,2)) # to increase plot margins in RStudio\nplot(ankle_curve, col='black', main='Kaplan Meier Curve', <\/em>\n  xlab='Follow Up Duration', ylab='Cumulative Survival Probability')<\/em><\/mark><\/code><\/pre>\n\n\n\n
\"\"<\/figure>\n\n\n\n
Again, giving the result 17%.<\/p>\n\n\n\n
These questions illustrate the importance of choosing a \u2018hard end point\u2019. They also show that if all patients lost to follow up are counted as a failure (worst-case scenario<\/em>), the survival curve drops steeply.<\/p>\n\n\n\n
There was not much difference in the 10-year survival for \u2018revision\u2019 and \u2018revision for aseptic loosening\u2019. This is because the one patient who had an infection had this relatively early (at 1.2 years). At this time there were 18 patients who had a longer follow up. Consequently, the effect on the failure rate (1 out of 19) is not as big as it would have been if only 1 patient had a longer follow up (1 out of 2).<\/p>\n\n\n\n
So, a failure at the \u2018tail end\u2019 (longer follow up) of the survival curve has a far more pronounced effect on the cumulative survival than a failure at the beginning (short follow up).<\/strong><\/p>\n\n\n\n
13. The plot can be created using console with the survival2<\/a><\/sup> package .<\/strong><\/p>\n\n\n\n
Download the data<\/a> and perform the analysis:<\/p>\n\n\n\n
load(\"\/path\/to_file\/survivalhip1.rda\")\nsurvivalhip1<\/mark>\n   Number   FU Censor Group\n1       1  5.1  FALSE Hip 1\n2       2  6.4   TRUE Hip 1\n3       3  7.8   TRUE Hip 1\n4       4  8.2   TRUE Hip 1\n5       5  9.3   TRUE Hip 1\n6       6 10.4  FALSE Hip 1\n7       7 10.5   TRUE Hip 1\n8       8 10.7  FALSE Hip 1\n9       9 10.9   TRUE Hip 1\n10     10 11.1  FALSE Hip 1\n11     11  4.6  FALSE Hip 2\n12     12  5.2  FALSE Hip 2\n13     13  6.9  FALSE Hip 2\n14     14  7.6  FALSE Hip 2\n15     15  8.1  FALSE Hip 2\n16     16  8.4  FALSE Hip 2\n17     17  8.9  FALSE Hip 2\n18     18  9.7   TRUE Hip 2\n19     19 10.6   TRUE Hip 2\n20     20 11.3  FALSE Hip 2\n<\/mark>hip_survival <- Surv(survivalhip1$FU, survivalhip1$Censor)\nhip_survival<\/mark>\n [1]  5.1+  6.4   7.8   8.2   9.3  10.4+ 10.5  10.7+ 10.9  11.1+\n[11]  4.6+  5.2+  6.9+  7.6+  8.1+  8.4+  8.9+  9.7  10.6  11.3+\n<\/mark>summary(hip_survival)<\/mark>\n      time            status   \n Min.   : 4.600   Min.   :0.0  \n 1st Qu.: 7.425   1st Qu.:0.0  \n Median : 8.650   Median :0.0  \n Mean   : 8.585   Mean   :0.4  \n 3rd Qu.:10.525   3rd Qu.:1.0  \n Max.   :11.300   Max.   :1.0  \n<\/mark>hip_curve <- survfit(hip_survival ~ Group, data=survivalhip1, conf.int=0.95)\nsummary(hip_curve)<\/mark>\nCall: survfit(formula = hip_survival ~ Group, data = survivalhip1, \n    conf.int = 0.95)\n\n                Group=Hip 1 \n time n.risk n.event survival std.err lower 95% CI upper 95% CI\n  6.4      9       1    0.889   0.105       0.7056        1.000\n  7.8      8       1    0.778   0.139       0.5485        1.000\n  8.2      7       1    0.667   0.157       0.4200        1.000\n  9.3      6       1    0.556   0.166       0.3097        0.997\n 10.5      4       1    0.417   0.173       0.1847        0.940\n 10.9      2       1    0.208   0.171       0.0418        1.000\n\n                Group=Hip 2 \n time n.risk n.event survival std.err lower 95% CI upper 95% CI\n  9.7      3       1    0.667   0.272       0.2995            1\n 10.6      2       1    0.333   0.272       0.0673            1<\/mark><\/em><\/code><\/pre>\n\n\n\n
Plot the curve:<\/p>\n\n\n\n
par(mar=c(2,2,2,2)) # to increase plot margins in RStudio\nplot(hip_curve, col=c('red', 'blue'), main='Kaplan Meier Curve', <\/em>\n  xlab='Follow Up Duration', ylab='Cumulative Survival Probability')<\/em><\/mark><\/code><\/pre>\n\n\n\n
\"\"<\/figure>\n\n\n\n
Hip 1 is red and Hip 2 is blue<\/p>\n\n\n\n
R labels the curves alphabetically and it may be necessary to check which curve is which or add a legend (i.e. with the survminer3<\/a><\/sup> package). From the output above, it can be seen that there were 6 failures in ‘Hip 1’ and two failures in ‘Hip 2’. Therefore, There should be 6 steps in the red curve and only two in the blue which is correct. The drops in the curve are steeper for ‘Hip 2’ because 7 patients were censored at the time of the first failure.<\/p>\n\n\n\n
14. This can be done using the console and the survival package4<\/a><\/sup>.<\/strong><\/p>\n\n\n\n
Using the hip_survival object from question 13:<\/p>\n\n\n\n
hip_survival <- Surv(survivalhip1$FU, survivalhip1$Censor)\nhip_survival<\/mark>\n [1]  5.1+  6.4   7.8   8.2   9.3  10.4+ 10.5  10.7+ 10.9  11.1+\n[11]  4.6+  5.2+  6.9+  7.6+  8.1+  8.4+  8.9+  9.7  10.6  11.3+\n<\/mark>summary(hip_survival)<\/mark>\n      time            status   \n Min.   : 4.600   Min.   :0.0  \n 1st Qu.: 7.425   1st Qu.:0.0  \n Median : 8.650   Median :0.0  \n Mean   : 8.585   Mean   :0.4  \n 3rd Qu.:10.525   3rd Qu.:1.0  \n Max.   :11.300   Max.   :1.0  <\/mark><\/em>\n<\/code><\/pre>\n\n\n\n
Perform the logrank test:<\/p>\n\n\n\n
hip_logrank <- survdiff(hip_survival ~ Group, data=survivalhip1)<\/em><\/span>\nhip_logrank<\/em><\/span>\nCall:<\/em><\/span>\nsurvdiff(formula = hip_survival ~ Group, data = survivalhip1)<\/em><\/span>\n             N Observed Expected (O-E)^2\/E (O-E)^2\/V<\/em><\/span>\nGroup=Hip 1 10        6     4.91     0.242     0.634<\/em><\/span>\nGroup=Hip 2 10        2     3.09     0.385     0.634<\/em><\/span>\n\nChisq= 0.6  on 1 degrees of freedom, p= 0.4<\/strong><\/span> <\/em><\/span><\/code><\/pre>\n\n\n\n
Perform Cox Proportional Hazards:<\/p>\n\n\n\n
hip_coxph <- coxph(hip_survival ~ Group, data=survivalhip1)<\/em><\/span>\nhip_coxph<\/em><\/span>\nCall:<\/em><\/span>\ncoxph(formula = hip_survival ~ Group, data = survivalhip1)<\/em><\/span>\n             coef exp(coef) se(coef)     z    p<\/em><\/span>\nGroupHip 2 -0.642     0.526    0.820 -0.78 0.43<\/em><\/span>\nLikelihood ratio test=0.67  on 1 df, p=0.412<\/strong><\/em><\/span>\nn= 20, number of events= 8<\/span> <\/em><\/span><\/code><\/pre>\n\n\n\n
And parametric Weibull:<\/p>\n\n\n\n
hip_parametric <- survreg(hip_survival ~ Group , dist = \"weibull\"<\/em><\/span>, <\/mark><\/em>\n  data=survivalhip1)<\/mark><\/em>\nhip_parametric<\/em><\/span>\nCall:<\/em><\/span>\nsurvreg(formula = hip_survival ~ Group, data = survivalhip1,<\/em><\/span>\n    dist = \"weibull\")<\/em><\/span>\nCoefficients:<\/em><\/span>\n(Intercept)  GroupHip 2 <\/em><\/span>\n 2.36525220  0.08988363 <\/em><\/span>\nScale= 0.1416005 <\/em><\/span>\nLoglik(model)= -20.5   Loglik(intercept only)= -20.8<\/em><\/span>\n    Chisq= 0.66 on 1 degrees of freedom, p= 0.4<\/strong><\/em><\/span>16<\/mark><\/strong>\nn= 20<\/span> <\/em><\/span><\/code><\/pre>\n\n\n\n
Similarly, perform the exponential parametric test which should show the following output (parts omitted for brevity):<\/p>\n\n\n\n
hip_parametric <- survreg(hip_survival ~ Group , dist = \"exponential<\/em><\/span>\"<\/em><\/span>, <\/mark><\/em>\n  data=survivalhip1)<\/mark><\/em>\nhip_parametric<\/em><\/span>\nCall:<\/em><\/span>\nsurvreg(formula = hip_survival ~ Group, data = survivalhip1,<\/em><\/span>\n    dist = \"exponential\")<\/em><\/span>\nCoefficients:<\/em><\/span>\n(Intercept)  GroupHip 2 <\/em><\/span>\n   2.712485    0.992514 <\/em><\/span>\nScale fixed at 1 <\/em><\/span>\nLoglik(model)= -31.7   Loglik(intercept only)= -32.5<\/em><\/span>\n    Chisq= 1.69 on 1 degrees of freedom, p= 0.193<\/strong> <\/em><\/span>\nn= 20<\/em><\/span><\/code><\/pre>\n\n\n\n
There is no statistical significant difference using the log rank test (p=0.42), cox proportional hazards model (p=0.412) and parametric regression analysis using the exponential (p=0.193) and Weibull (p=0.416) distributions between the survival curves of both hips. Therefore, the null hypothesis can’t be rejected and it should be concluded that there is no difference in the outcome of both hips.<\/p>\n\n\n\n
<\/p>\n","protected":false},"excerpt":{"rendered":"
The table below shows 20 patients who have been diagnosed with cancer. The first column shows the follow up (in years) of the patients who were alive at review. In the second column, the follow up till time of death is indicated. The third column shows the time to last review in the patients who […]<\/p>\n","protected":false},"author":2,"featured_media":0,"parent":0,"menu_order":0,"comment_status":"closed","ping_status":"closed","template":"","meta":{"inline_featured_image":false,"footnotes":""},"class_list":["post-1196","page","type-page","status-publish","hentry"],"_links":{"self":[{"href":"https:\/\/pcool.dyndns.org\/index.php\/wp-json\/wp\/v2\/pages\/1196","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/pcool.dyndns.org\/index.php\/wp-json\/wp\/v2\/pages"}],"about":[{"href":"https:\/\/pcool.dyndns.org\/index.php\/wp-json\/wp\/v2\/types\/page"}],"author":[{"embeddable":true,"href":"https:\/\/pcool.dyndns.org\/index.php\/wp-json\/wp\/v2\/users\/2"}],"replies":[{"embeddable":true,"href":"https:\/\/pcool.dyndns.org\/index.php\/wp-json\/wp\/v2\/comments?post=1196"}],"version-history":[{"count":29,"href":"https:\/\/pcool.dyndns.org\/index.php\/wp-json\/wp\/v2\/pages\/1196\/revisions"}],"predecessor-version":[{"id":4447,"href":"https:\/\/pcool.dyndns.org\/index.php\/wp-json\/wp\/v2\/pages\/1196\/revisions\/4447"}],"wp:attachment":[{"href":"https:\/\/pcool.dyndns.org\/index.php\/wp-json\/wp\/v2\/media?parent=1196"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}