1 Example Analysis for Continuous Outcomes
1.1 Training the XGBoostSub_con Model
We will begin by training the XGBoostSub_con model for a continuous outcome. This involves using the A-learning or Weight-learning loss function. Note that the treatment variable must be converted to (1, -1). If your treatment variable is (1, 0), you will need to convert it accordingly.
X_feature=tutorial_data[,c("x1", "x2" ,"x3" ,"x4","x5","x6","x7","x8","x9","x10")]
y_label=tutorial_data$y.con
# Convert treatment variable to (1 -1) format (1 for treatment, -1 for control)
trt=ifelse(tutorial_data$treatment==1, 1, -1)
true_label=tutorial_data$subgroup_label
# Estimate the propensity score using logistic regression
data_logit=tutorial_data[,c("treatment","x1", "x2" ,"x3" ,"x4" , "x5", "x6" , "x7","x8","x9","x10")]
logit_model <- glm(treatment~ ., data = data_logit, family = binomial)
pi <- predict(logit_model, type = "response")
# Train the XGBoostSub_con model using the specified parameters
model <- XGBoostSub_con(X_feature, y_label, trt ,pi,Loss_type = "A_learning", params = list(learning_rate = 0.005, max_depth = 1, lambda = 5, tree_method = 'hist'), nrounds = 200, disable_default_eval_metric = 0, verbose = FALSE)After training the model, you can print the loss value at each epoch using the following code:
cat("Evaluation loss:\n")
#> Evaluation loss:
print(model$evaluation_log)
#> iter train_A_loss
#> <num> <num>
#> 1: 1 1.375787
#> 2: 2 1.375424
#> 3: 3 1.375065
#> 4: 4 1.374708
#> 5: 5 1.374355
#> ---
#> 196: 196 1.339571
#> 197: 197 1.339486
#> 198: 198 1.339402
#> 199: 199 1.339319
#> 200: 200 1.339236You can also evaluate the loss on independent datasets using the eval_metric_con function. In this example, we’ll use the training data as test data for illustration purposes. However, you may use a new dataset to see the loss based on the trained model.
eval_metric_test <- eval_metric_con(model, X_feature, y_label, pi, trt, Loss_type = "A_learning")
cat("Testing Evaluation Result:\n")
#> Testing Evaluation Result:
print(eval_metric_test$value)
#> [1] 1.3392361.2 Evaluating Predictive Biomarker Importance with XGBoostSub_con
Next, we evaluate the importance of predictive biomarkers using the predictive_biomarker_imp function.
biomarker_imp=predictive_biomarker_imp(model)
kable_styling(kable(x=biomarker_imp, format="html", caption = "Biomarker importance table"),
bootstrap_options="striped",font_size=16)| Biomarker | Importance |
|---|---|
| x2 | 0.4904407 |
| x1 | 0.4793201 |
| x10 | 0.0302392 |
From the results, we observe that the most important predictive biomarker is x2.
1.3 Obtaining subgroup results with XGBoostSub_con
You can obtain the subgroup results using the get_subgroup_results function.
subgroup_results=get_subgroup_results(model, X_feature, subgroup_label=true_label, cutoff = 0.5)
kable_styling(kable(x=head(subgroup_results$assignment), format="html", caption = "The first 6 subjects subgroup assigentment"),
bootstrap_options="striped",font_size=16)| Patient_ID | Treatment_Assignment |
|---|---|
| 1 | 1 |
| 2 | 1 |
| 3 | 1 |
| 4 | 0 |
| 5 | 1 |
| 6 | 1 |
cat("Prediction accuracy of true subgroup label:\n")
#> Prediction accuracy of true subgroup label:
cat(subgroup_results$acc)
#> 0.634Note that when using the get_subgroup_results function with real data, the true subgroup label is typically unknown. In such cases, set subgroup_label = NULL. This will return only the subgroup assignment without the prediction accuracy of the subgroup.
1.4 Commonly used subgroup/biomarker analysis tools
Given that x2 was identified as the most important predictive biomarker by XGBoostSub_con, we will use the cdf_plot function to evaluate the distribution of this biomarker by assessing the percentage of subjects falling within different cutoff values.
cdf_plot (xvar="x2", data=tutorial_data, y.int=5, xlim=NULL, xvar.display="Biomarker_X2", group=NULL)
We will also evaluate x2 prognostic capability. To determine if this biomarker is prognostic, we can assess the association between the response variable y.con and the biomarker x2 using the scat_cont_plot function.
scat_cont_plot(
yvar='y.con',
xvar='x2',
data=tutorial_data,
ybreaks = NULL,
xbreaks = NULL,
yvar.display = 'y',
xvar.display = "Biomarker_X2"
)
#> $Cor.Pearson
#> [1] -0.12
#>
#> $Cor.Spearman
#> [1] -0.12
#>
#> $fig
#>
#> $slope
#> [1] -0.1013384
#>
#> $intercept
#> [1] 0.6769244Selecting the cutoff value for a predictive biomarker is important in clinical practice, as it can, for example, help enrich responders. We select the optimal cutoff for the identified biomarker x2 from a list of candidate cutoff values using the fixcut_con function.
cutoff_result_con=fixcut_con(yvar='y.con', xvar="x2", dir=">", cutoffs=c(-0.1,-0.3,-0.5,0.1,0.3,0.5), data=tutorial_data, method="t.test", yvar.display='y.con', xvar.display='biomarker x2', vert.x=F)
cutoff_result_con$fig
From the results, the optimal cutoff value is 0.5. Next, we define the biomarker positive group using a cutoff of 0.5. The biomarker positive group includes subjects with x2 values less than or equal to 0.5, while the biomarker negative group includes subjects with x2 values greater than 0.5. We first use the cut_perf() function to evaluate the performance between the biomarker positive and negative groups.
res=cut_perf (yvar="y.con", censorvar=NULL, xvar="x2", cutoff=c(0.5), dir="<=", xvars.adj=NULL, data=tutorial_data, type='c', yvar.display='y.con', xvar.display="biomarker x2")
kable_styling(kable(x=res$comp.res.display[,-2:-4], format="html", caption = "Performace at optimal cutoff"),
bootstrap_options="striped",font_size=16)| Response | Median.Diff | Mean.Diff.CI | t.pvalue | WRS.pvalue |
|---|---|---|---|---|
| y.con | -0.3 | -0.38 (NA , NA) | NA | NA |
Next, we will further assess the predictive model performance by examining the differences between the treatment and control groups within both the biomarker positive and negative groups.
tutorial_data$biogroup=ifelse(tutorial_data$x2<=0.5,'biomarker_positive','biomarker_negative')
res = subgrp_perf_pred(yvar="y.time", censorvar="y.event", grpvar="biogroup", grpname=c("biomarker_positive",'biomarker_negative'),trtvar="treatment_categorical", trtname=c("Placebo" , "Treatment"), xvars.adj=NULL,data=tutorial_data, type="s")
kable_styling(kable(x=res$group.res.display, format="html", caption = "BioSubgroup Stat"),
bootstrap_options="striped",font_size=16)| Response | N | Events | Mean.Surv.CI | Med.Surv.CI | |
|---|---|---|---|---|---|
| treatment_categorical=Placebo,biogroup=biomarker_positive | y.time | 321 | 127 | 3.7 (2.49 , 4.92) | 1.99 (1.65 , 2.7) |
| treatment_categorical=Placebo,biogroup=biomarker_negative | y.time | 192 | 76 | 4.49 (2.37 , 6.61) | 1.83 (1.58 , 2.27) |
| treatment_categorical=Treatment,biogroup=biomarker_positive | y.time | 321 | 133 | 4.84 (3.66 , 6.02) | 2.85 (2.11 , 3.52) |
| treatment_categorical=Treatment,biogroup=biomarker_negative | y.time | 166 | 75 | 2.63 (1.93 , 3.32) | 1.51 (0.84 , 2.14) |
res$fig
From the Kaplan-Meier (KM) curve, a difference between the treatment and control groups is observed in the biomarker positive group. This difference is not observed in the biomarker negative group, validating that biomarker x2 identified by XGBoostSub_con is predictive.
Here, users can select the optimal cutoff based on various metrics such as sensitivity, specificity, PPV (Positive Predictive Value), and NPV (Negative Predictive Value). The choice of metric depends on the specific research context. For instance, if the goal is to enrich responders, focusing on PPV might be more appropriate. Based on the highest PPV value, users can then select the optimal cutoff.
From the KM curve, it appears that 0.5 for 
