Skip to content
Closed
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
134 changes: 131 additions & 3 deletions Comp/r-sas_ancova.qmd
Original file line number Diff line number Diff line change
Expand Up @@ -8,13 +8,141 @@ title: "R vs SAS ANCOVA"
knitr::opts_chunk$set(echo = TRUE)
```

# ANCOVA Comparison
# Summary
Comment thread
DrLynTaylor marked this conversation as resolved.

## Goal

The goal of this comparison is to evaluate whether ANCOVA produces equivalent results in R and SAS using the same dataset and model specification. Results from the R `stats::lm() function`, the `sasLM::GLM()` function, and SAS `PROC GLM` are compared for sums of squares and least-squares means.

## Scope

::::::: columns
:::: {.column width="45%"}
::: {.callout-note appearance="minimal" collapse="false"}
## Methodologies

- Analysis of covariance (ANCOVA)

- Type I sums of squares

- Type III sums of squares

- Least-squares means
:::
::::

:::: {.column width="55%"}
::: {.callout-note appearance="minimal" collapse="false"}
## Technical implementations

- SAS: `PROC GLM`
- R: `stats::lm()`
- R: `sasLM::GLM()` with `EMEANS=TRUE`
:::
::::
:::::::

# Prerequisites

## R packages

```{r}
#| output: false
#| message: false
library(stats) # lm()
library(sasLM) # GLM()
library(tibble) # tibble()
library(tidyverse) # tidyverse functions
```

## Data

The same dataset is used for the R and SAS analyses.

R:

```{r}
df_sas <- tibble(
drug = factor(c( "A", "A", "A", "A", "A", "A", "A", "A", "A", "A",
"D", "D", "D", "D", "D", "D", "D", "D", "D", "D",
"F", "F", "F", "F", "F", "F", "F", "F", "F", "F"
)),
pre = c(
11, 8, 5, 14, 19, 6, 10, 6, 11, 3,
6, 6, 7, 8, 18, 8, 19, 8, 5, 15,
16, 13, 11, 9, 21, 16, 12, 12, 7, 12
),
post = c( 6, 0, 2, 8, 11, 4, 13, 1, 8, 0,
0, 2, 3, 1, 18, 4, 14, 9, 1, 9,
13, 10, 18, 5, 23, 12, 5, 16, 1, 20
)
)
```

SAS:

``` sas
data df_sas;
input drug $ pre post;
datalines;
A 11 6
A 8 0
A 5 2
A 14 8
A 19 11
A 6 4
A 10 13
A 6 1
A 11 8
A 3 0
D 6 0
D 6 2
D 7 3
D 8 1
D 18 18
D 8 4
D 19 14
D 8 9
D 5 1
D 15 9
F 16 13
F 13 10
F 11 18
F 9 5
F 21 23
F 16 12
F 12 5
F 12 16
F 7 1
F 12 20
;
run;
```

The ANCOVA model includes drug as the categorical treatment variable and pre as the continuous covariate.

The general ANCOVA model is:

$$ \text{post} = \text{drug} + \text{pre} $$

where:

- **post** is the response variable

- **drug** is the treatment variable

- **pre** is the continuous covariate.

## ANCOVA Comparison

The following table shows the types of One Sample t-test analysis, the capabilities of each language, and whether or not the results from each language match.

| Analysis | Supported in R | Supported in SAS | Results Match | Notes |
|--------------|--------------|--------------|--------------|-----------------|
|---------------|---------------|---------------|---------------|---------------|
| ANCOVA using general linear model and lsmeans | [Yes](../R/ancova.html) | [Yes](../SAS/ancova.html) | [Yes](#comp) | `GLM()` function from **sasLM** with `EMEANS=TRUE` is the easiest to use and matches SAS |
| Type I Sum of Squares | [Yes](#comp) | [Yes](#comp) | [Yes](#comp) | Sequential sums of squares |
| Type II Sum of Squares | [Yes](#comp) | [Yes](#comp) | [Yes](#comp) | `sasLM::GLM()` provides SAS-style output |
| Type III Sum of Squares | [Yes](#comp) | [Yes](#comp) | [Yes](#comp) | `sasLM::GLM()` with `EMEANS=TRUE` matches SAS `LSMEANS` |

## Comparison Results {#comp}

Expand All @@ -40,4 +168,4 @@ R `lm()` documentation: <https://www.rdocumentation.org/packages/stats/versions/

R `GLM()` documentation: <https://cran.r-project.org/web/packages/sasLM/sasLM.pdf>

SAS `PROC GLM` documentation: <https://documentation.sas.com/doc/en/statug/15.2/statug_glm_syntax01.htm>
SAS `PROC GLM` documentation: <https://documentation.sas.com/doc/en/statug/15.2/statug_glm_syntax01.htm>
54 changes: 38 additions & 16 deletions SAS/weighted_logrank_test.qmd
Original file line number Diff line number Diff line change
Expand Up @@ -6,32 +6,37 @@ execute:
jupyter: false
---

```{r}
#| label: setup
#| include: false
library(tibble)
```

# Introduction

In clinical studies with time-to-event outcomes, it is commonly assumed that the hazard functions of two groups are proportional. The standard log-rank test is widely used to test the equivalence of survival functions. It assumes no relationship between the censoring indicator and the forecast, such that dropouts or individuals who are not fully observed are perceived to be the same in risks (sensitivity to proportional risk); the events are also weighted equally regardless of time of occurrence ( early/ later). However, several scenarios can lead to non-proportional hazards (NPH). For example, a delayed treatment effect may be observed in the treatment arm, which can lead to a departure from proportionality of the survival curves. Log-rank test loses power under the non-proportional hazard assumption. Thus, there are many tests available in the literature that can handle such scenarios. Most commonly used tests are as follows:

- Weighted log-rank test.
- Restricted Mean Survival Time (RMST).
- Milestone survival.
- Max-Combo test.
- Modestly Weighted log-rank test.
- Weighted log-rank test.
- Restricted Mean Survival Time (RMST).
- Milestone survival.
- Max-Combo test.
- Modestly Weighted log-rank test.

This particular document focuses mainly on the weighted log-rank test, Max-combo test and Modestly Weighted log-rank.

# Data.

\
Data source: <https://stats.idre.ucla.edu/sas/seminars/sas-survival/>

The data include 500 subjects from the Worcester WHAS WHAS Study. This study examined several factors, such as age, gender and BMI, that may influence survival time after an WHAS. Follow-up time for all participants begins at the time of hospital admission after an WHAS and ends with death or loss to follow-up (censoring). The variables used here are:
The data include 500 subjects from the Worcester WHAS Study. This study examined several factors, such as age, gender and BMI, that may influence survival time after an WHAS. Follow-up time for all participants begins at the time of hospital admission after an WHAS and ends with death or loss to follow-up (censoring). The variables used here are:

- lenfol: length of follow-up, terminated either by death or censoring - time variable
- lenfol: length of follow-up, terminated either by death or censoring - time variable

- fstat: loss to followup = 0, death = 1 - censoring variable
- fstat: loss to followup = 0, death = 1 - censoring variable

- afb: atrial fibrillation, no = 0, 1 = yes - Covariate
- afb: atrial fibrillation, no = 0, 1 = yes - Covariate

- gender: males = 0, females = 1 - stratification factor
- gender: males = 0, females = 1 - stratification factor

**Kaplan-Meier Survival Curve.**

Expand Down Expand Up @@ -320,10 +325,10 @@ knitr::include_graphics("../images/weighted_log_rank/stratifed MaxC test stat.pn

To change the choice of weight, modify the `weight =` list in the macro to any desired ($\rho,\gamma$) pairs. The macro supports combined multiple weights `weights = %str(0,0 1,0 0,1 1,1)` by default. These correspond to the Fleming–Harrington class $G(\rho,\gamma)$, where:

- $G(0,0)$: standard log-rank test (equal weighting over time)
- $G(1,0)$: emphasises early differences
- $G(0,1)$: emphasises late differences
- $G(1,1)$: emphasises middle-time differences
- $G(0,0)$: standard log-rank test (equal weighting over time)
- $G(1,0)$: emphasises early differences
- $G(0,1)$: emphasises late differences
- $G(1,1)$: emphasises middle-time differences

```{sas}
%combo_wlr(
Expand All @@ -344,7 +349,24 @@ To change the choice of weight, modify the `weight =` list in the macro to any d
knitr::include_graphics("../images/weighted_log_rank/weight choice.png")
```

# 3. Modestly Weighted Log-rank test.
# 3. Summary Results for Weighted Log-rank.

Here are the results using SAS statistical methods with the data from the Worcester WHAS study.

| Test | Function | Chi-square | P-Value |
|----------------|------------------------|----------------|----------------|
| **Standard Log-Rank** | `PROC LIFETEST / TEST=LOGRANK` | 7.7911 | 0.0053 |
| **Peto - Peto** | `PROC LIFETEST / TEST=PETO` | 6.7602 | 0.0093 |
| **Modified Peto-Peto** | `PROC LIFETEST / TEST=MODPETO` | 6.7611 | 0.0093 |
| **Tarone-Ware** | `PROC LIFETEST / TEST=TARONEWARE` | 6.6664 | 0.0098 |
| **Gehan- Breslow / Wilcoxon** | `PROC LIFETEST / TEST=WILCOXON` | 5.5370 | 0.0186 |
| **Fleming-Harrington (0.5, 0.5)** | `PROC LIFETEST / TEST=FH(0.5, 0.5)` | 8.6661 | 0.0032 |
| **Fleming-Harrington (1, 1)** | `PROC LIFETEST / TEST=FH(1, 1)` | 8.2442 | 0.0041 |
| **Fleming-Harrington (0, 1)** | `PROC LIFETEST / TEST=FH(0, 1)` | 8.0896 | 0.0045 |
| **Fleming-Harrington (0.5, 2)** | `PROC LIFETEST / TEST=FH(0.5, 2)` | 7.1724 | 0.0074 |
| **Fleming-Harrington (1, 0)** | `PROC LIFETEST / TEST=FH(1, 0)` | 6.7309 | 0.0095 |

# 4. Modestly Weighted Log-rank test.

This test is suitable for a delayed effect scenario since it gives more weight to later event times; regardless, weighting is controlled, such that the worse effect of treatment is not rewarded at an early time point. However, in a scenario where there is an early treatment effect which diminishes over time, modestly weighted Log-Rank losses sizeable power compared to the log-rank test. To the knowledge of CAMIS contributors, there is no direct implementation of this test in `PROC LIFETEST` for SAS.

Expand Down
Loading