HW3a

      # Dependencies
      library(rstan)

Loading required package: StanHeaders

rstan version 2.35.0.9000 (Stan version 2.35.0)

For execution on a local, multicore CPU with excess RAM we recommend calling
options(mc.cores = parallel::detectCores()).
To avoid recompilation of unchanged Stan programs, we recommend calling
rstan_options(auto_write = TRUE)
For within-chain threading using `reduce_sum()` or `map_rect()` Stan functions,
change `threads_per_chain` option:
rstan_options(threads_per_chain = 1)

      library(arrow)

Attaching package: 'arrow'

The following object is masked from 'package:utils':

    timestamp

      library(survival)
      library(dplyr)

Attaching package: 'dplyr'

The following objects are masked from 'package:stats':

    filter, lag

The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union

      library(ggplot2)
      library(eha)

Loading Data

# Load Data
      df <- read_parquet("age-stage-survival.parquet")

      # Prepare data for Stan
      df$status <- 1  # All events occurred, no censoring
      df$stage <- as.integer(df$grade)
      df <- df %>% 
        mutate(
          stage1 = ifelse(stage ==1,1,0),
          stage2 = ifelse(stage ==2,1,0),
          stage3 = ifelse(stage ==3,1,0)  
          )
      
      
      stan_data <- list(
        N = nrow(df),            # Number of observations
        K = length(unique(df$grade)),  # Number of unique stages
        time = df$time,          # Survival times
        status = df$status,      # Event indicators (1 = event, 0 = censored)
        age = df$age,            # Age covariate
        stage = df$stage,         # Stage
        stage1 = df$stage1,
        stage2 = df$stage2,
        stage3 = df$stage3
      )

1. A single Weibull baseline hazard with a common age relationship across stage.

1a. Write a mathematical expression of the model.

The log hazard function can be expressed as:

\[ log(h(t|x)) = \beta_0 + \beta_1 \text{age} + \beta_2 \text{age}^2 + \alpha_1 I_1 + \alpha_2 I_2 + \alpha_3 I_3 + log(h_0(t, shape, scale)) \]

Let’s write the model. In stan for MLE, we need to define these 3 sections

writeLines("

data {
  int<lower=1> N;  // Number of observations
  vector[N] age;  // Age variable
  array[N] int<lower=0, upper=1> stage1;
  array[N] int<lower=0, upper=1> stage2;
  array[N] int<lower=0, upper=1> stage3;
  vector[N] time;  // Time until death variable
}

parameters {
  real alpha1;  // Coefficients for stages 2, 3, and 4 (stage 1 is the baseline)
  real alpha2;
  real alpha3;
  real beta0;
  real beta1;  // Coefficient for age
  real beta2;  // Coefficient for squared age
  real<lower=0> shape;  // Shape parameter of the Weibull distribution
  real<lower=0> scale;  // Scale parameter of the Weibull distribution
}

model {
  for (i in 1:N) {  // Loop through each observation
    // Calculate the linear predictor based on age and stage
    // Add the stage coefficient; stage 4 is the baseline, so we only adjust for stages 1-3
    real linear_pred = beta0 + beta1 * age[i] + beta2 * square(age[i])
                        + alpha1 * stage1[i]
                        + alpha2 * stage2[i]
                        + alpha3 * stage3[i]
                        ;
    
    target += weibull_lpdf(time[i] | shape, scale) +linear_pred + (exp(linear_pred)-1) * weibull_lccdf(time[i] | shape, scale);  
  }
}
", con = "model1.stan")

Fitting the model

# Compile the Stan model
stan_model_1 <- stan_model(file = "model1.stan")

Trying to compile a simple C file

Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
using C compiler: ‘Apple clang version 14.0.3 (clang-1403.0.22.14.1)’
using SDK: ‘MacOSX13.3.sdk’
clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/opt/R/arm64/include    -fPIC  -falign-functions=64 -Wall -g -O2  -c foo.c -o foo.o
In file included from <built-in>:1:
In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
#include <cmath>
         ^~~~~~~
1 error generated.
make: *** [foo.o] Error 1

# Perform MLE using the optimizing function
fit_mle_1 <- optimizing(stan_model_1, data = stan_data, hessian = TRUE)

Warning in .local(object, ...): non-zero return code in optimizing

# View the optimized parameter estimates
print(fit_mle_1)

$par
    alpha1     alpha2     alpha3      beta0      beta1      beta2      shape 
-1.7152824 -0.6525107  1.5985501  0.7469437  0.9868907 -0.5541252  0.5455413 
     scale 
 1.7387412 

$value
[1] -649445

$return_code
[1] 70

$hessian
       alpha1 alpha2 alpha3 beta0 beta1 beta2     shape         scale
alpha1      0      0      0     0     0     0    0.0000  0.000000e+00
alpha2      0      0      0     0     0     0    0.0000  0.000000e+00
alpha3      0      0      0     0     0     0    0.0000  0.000000e+00
beta0       0      0      0     0     0     0    0.0000  0.000000e+00
beta1       0      0      0     0     0     0    0.0000  0.000000e+00
beta2       0      0      0     0     0     0    0.0000  0.000000e+00
shape       0      0      0     0     0     0  177.7217 -2.727707e+02
scale       0      0      0     0     0     0 -272.7707  8.526513e-11

$theta_tilde
        alpha1     alpha2  alpha3     beta0     beta1      beta2     shape
[1,] -1.715282 -0.6525107 1.59855 0.7469437 0.9868907 -0.5541252 0.5455413
        scale
[1,] 1.738741

1b. Generate a partial effect plot of log relative hazard by age, stratified by stage. (Uncertainty interval not required for now.)

# Extract MLE estimates for parameters
beta0 <- fit_mle_1$par["beta0"]
beta1 <- fit_mle_1$par["beta1"]
beta2 <- fit_mle_1$par["beta2"]
alpha1 <- fit_mle_1$par["alpha1"]
alpha2 <- fit_mle_1$par["alpha2"]
alpha3 <- fit_mle_1$par["alpha3"]
shape <- fit_mle_1$par["shape"]
scale <- fit_mle_1$par["scale"]


# Define the age range for plotting
age_range <- seq(20, 80, by = 1)  # Modify this as needed

# Weibull baseline hazard for each age point
weibull_baseline <- (shape - 1) * log(age_range / scale)

# Calculate log relative hazards for each stage
log_relative_hazard_stage1 <- beta0 + beta1 * age_range + beta2 * age_range^2 + alpha1 + weibull_baseline
log_relative_hazard_stage2 <- beta0 + beta1 * age_range + beta2 * age_range^2 + alpha2 + weibull_baseline
log_relative_hazard_stage3 <- beta0 + beta1 * age_range + beta2 * age_range^2 + alpha3 + weibull_baseline
log_relative_hazard_stage4 <- beta0 + beta1 * age_range + beta2 * age_range^2 + weibull_baseline  # baseline stage


# Create the plot
plot(age_range, log_relative_hazard_stage1, type = "l", col = "blue", lty = 1,
     ylim = range(log_relative_hazard_stage1, log_relative_hazard_stage2, log_relative_hazard_stage3, log_relative_hazard_stage4),
     xlab = "Age", ylab = "Log Relative Hazard", main = "Log Relative Hazard vs Age by Stage")
lines(age_range, log_relative_hazard_stage2, col = "red", lty = 2)
lines(age_range, log_relative_hazard_stage3, col = "green", lty = 3)
lines(age_range, log_relative_hazard_stage4, col = "purple", lty = 4)

# Add a legend
legend("topleft", legend = c("Stage 1", "Stage 2", "Stage 3", "Stage 4 (B)"),
       col = c("blue", "red", "green", "purple"), lty = 1:4)

1c. Generate a partial effect plot of median survival by age, stratified by stage. (Uncertainty interval not required for now.) There is no reason to think the log relative hazard is linear in age.

# Extract MLE estimates for parameters
beta0 <- fit_mle_1$par["beta0"]
beta1 <- fit_mle_1$par["beta1"]
beta2 <- fit_mle_1$par["beta2"]
alpha1 <- fit_mle_1$par["alpha1"]
alpha2 <- fit_mle_1$par["alpha2"]
alpha3 <- fit_mle_1$par["alpha3"]
shape <- fit_mle_1$par["shape"]
scale <- fit_mle_1$par["scale"]

# Define age range and stages
ages <- seq(20, 90, by = 5)  # Adjust age range as needed
stages <- c(1, 2, 3, 4)

# Function to calculate median survival using exponential_lccdf
median_survival <- function(age, stage_indicator) {
  # Select the appropriate alpha for each stage
  stage_alpha <- c(0, alpha1, alpha2, alpha3)[stage_indicator + 1]  # +1 to align with R indexing
  
  # Calculate the log hazard rate for the given age and stage & convert log-hazard to median survival using the exponential lccdf
  log_hazard <- beta0 + beta1 * age + beta2 * age^2 #+ stage_alpha-exponential_lccdf(0.5, log_hazard, shape, scale)
  
}

# Generate data for median survival by age and stage
survival_data <- expand.grid(age = ages, stage = stages)
survival_data$median_survival <- mapply(median_survival, survival_data$age, survival_data$stage)

# Plot
library(ggplot2)

ggplot(survival_data, aes(x = age, y = median_survival, color = as.factor(stage))) +
  geom_line() +
  labs(
    x = "Age",
    y = "Median Survival",
    color = "Stage",
    title = "Partial Effect of Age on Median Survival by Stage"
  ) +
  theme_minimal()

2. A single Weibull baseline hazard with a common age relationship across stage.

2a.Write a mathematical expression of the model. The hazard function is defined as:

\[ log(h(t|x)) = \beta_0 + \beta_1 (\text{age}) + \beta_2 (\text{age}^2) + \alpha_0 I_1 + \alpha_1 I_1 (\text{age}) + \alpha_2 I_1 (\text{age}^2) + \gamma_0 I_2 + \gamma_1 I_2 (\text{age}) + \gamma_2 I_2 (\text{age}^2) + \lambda_0 I_3 + \lambda_1 I_3 (\text{age}) + \lambda_2 I_3 (\text{age}^2) + log(h_0(t,shape,scale)) \]

writeLines("
 
data {
   int<lower=1> N;  // Number of observations
   vector[N] age;  // Age variable
   array[N] int<lower=0, upper=1> stage2;
   array[N] int<lower=0, upper=1> stage3;
   array[N] int<lower=0, upper=1> stage1;
   vector[N] time;  // Time until death variable
 }
 
 parameters {
   real beta0; // constant coefficient stage 4
   real beta1;  // Coefficient for age stage 4
   real beta2;  // Coefficient for squared age stage 4
   real gamma0; //coefficient for stage 2 constant
   real gamma1; //cofefficient for stage 2 age
   real gamma2; // coefficient for stage 2 age squared
   real lambda0; //coefficient for stage 3 constant
   real lambda1; //cofefficient for stage 3 age
   real lambda2; // coefficient for stage 3 age squared
   real alpha0; //coefficient for stage 1 constant
   real alpha1; //cofefficient for stage 1 age
   real alpha2; // coefficient for stage 1 age squared
   real<lower=0> shape;  // Shape parameter of the Weibull distribution
   real<lower=0> scale;  // Scale parameter of the Weibull distribution
 }
 

 
 model {
   for (i in 1:N) {  // Loop through each observation
     // Calculate the linear predictor based on age and stage
     // Add the stage coefficient; stage 4 is the baseline, so we only adjust for stages 1-3
     real linear_pred = beta0 + beta1 * age[i] + beta2 * square(age[i])
                         + alpha0  * stage1[i]  + alpha1 * stage1[i] * age[i] + alpha2 * stage1[i] * square(age[i])
                         + gamma0 * stage2[i]  + gamma1 *  stage2[i] * age[i] + gamma2  * stage2[i] * square(age[i])
                         +  lambda0  *stage3[i]  + lambda1 *  stage3[i] * age[i] + lambda2  * stage3[i] * square(age[i])
                         ;
     
     target += weibull_lpdf(time[i] | shape, scale) +linear_pred + (exp(linear_pred)-1) * weibull_lccdf(time[i] | shape, scale);  
   }
 }
 ", con = "model2.stan")

# Compile the Stan model
stan_model_2 <- stan_model(file = "model2.stan")

Trying to compile a simple C file

Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
using C compiler: ‘Apple clang version 14.0.3 (clang-1403.0.22.14.1)’
using SDK: ‘MacOSX13.3.sdk’
clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/opt/R/arm64/include    -fPIC  -falign-functions=64 -Wall -g -O2  -c foo.c -o foo.o
In file included from <built-in>:1:
In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
#include <cmath>
         ^~~~~~~
1 error generated.
make: *** [foo.o] Error 1

# Perform MLE using the optimizing function
fit_mle_2 <- optimizing(stan_model_2, data = stan_data, hessian = TRUE)

Warning in .local(object, ...): non-zero return code in optimizing

# View the optimized parameter estimates
print(fit_mle_2)

$par
      beta0       beta1       beta2      gamma0      gamma1      gamma2 
-0.33224930  0.28368981 -0.69530382  1.98379009 -1.05746328 -0.01125933 
    lambda0     lambda1     lambda2      alpha0      alpha1      alpha2 
-1.59085175  1.00515259 -0.10365003 -0.37114492  0.27803046 -0.48446769 
      shape       scale 
 4.29897227  2.80848936 

$value
[1] -1021448

$return_code
[1] 70

$hessian
        beta0 beta1 beta2 gamma0 gamma1 gamma2 lambda0 lambda1 lambda2 alpha0
beta0       0     0     0      0      0      0       0       0       0      0
beta1       0     0     0      0      0      0       0       0       0      0
beta2       0     0     0      0      0      0       0       0       0      0
gamma0      0     0     0      0      0      0       0       0       0      0
gamma1      0     0     0      0      0      0       0       0       0      0
gamma2      0     0     0      0      0      0       0       0       0      0
lambda0     0     0     0      0      0      0       0       0       0      0
lambda1     0     0     0      0      0      0       0       0       0      0
lambda2     0     0     0      0      0      0       0       0       0      0
alpha0      0     0     0      0      0      0       0       0       0      0
alpha1      0     0     0      0      0      0       0       0       0      0
alpha2      0     0     0      0      0      0       0       0       0      0
shape       0     0     0      0      0      0       0       0       0      0
scale       0     0     0      0      0      0       0       0       0      0
        alpha1 alpha2      shape         scale
beta0        0      0     0.0000  0.000000e+00
beta1        0      0     0.0000  0.000000e+00
beta2        0      0     0.0000  0.000000e+00
gamma0       0      0     0.0000  0.000000e+00
gamma1       0      0     0.0000  0.000000e+00
gamma2       0      0     0.0000  0.000000e+00
lambda0      0      0     0.0000  0.000000e+00
lambda1      0      0     0.0000  0.000000e+00
lambda2      0      0     0.0000  0.000000e+00
alpha0       0      0     0.0000  0.000000e+00
alpha1       0      0     0.0000  0.000000e+00
alpha2       0      0     0.0000  0.000000e+00
shape        0      0   369.8346 -2.149486e+03
scale        0      0 -2149.4863  1.202807e-07

$theta_tilde
          beta0     beta1      beta2  gamma0    gamma1      gamma2   lambda0
[1,] -0.3322493 0.2836898 -0.6953038 1.98379 -1.057463 -0.01125933 -1.590852
      lambda1  lambda2     alpha0    alpha1     alpha2    shape    scale
[1,] 1.005153 -0.10365 -0.3711449 0.2780305 -0.4844677 4.298972 2.808489

2b.Generate a partial effect plot of log relative hazard by age, stratified by stage. (Uncertainty interval not required for now.)

# Extract MLE estimates for parameters
beta0 <- fit_mle_2$par["beta0"]
beta1 <- fit_mle_2$par["beta1"]
beta2 <- fit_mle_2$par["beta2"]

gamma0 <- fit_mle_2$par["gamma0"]
gamma1 <- fit_mle_2$par["gamma1"]
gamma2 <- fit_mle_2$par["gamma2"]
lambda0 <- fit_mle_2$par["lambda0"]
lambda1 <- fit_mle_2$par["lambda1"]
lambda2 <- fit_mle_2$par["lambda2"]
alpha0 <- fit_mle_2$par["alpha0"]
alpha1 <- fit_mle_2$par["alpha1"]
alpha2 <- fit_mle_2$par["alpha2"]

shape <- fit_mle_2$par["shape"]
scale <- fit_mle_2$par["scale"]


# Define the age range for plotting
age_range <- seq(20, 80, by = 1)  # Modify this as needed

# Weibull baseline hazard for each age point
weibull_baseline <- (shape - 1) * log(age_range / scale)

# Calculate log relative hazards for each stage
log_relative_hazard_stage1 <- beta0 + beta1 * age_range + beta2 * age_range^2 + alpha0 + alpha1 * age_range + alpha2 * age_range^2 + weibull_baseline
log_relative_hazard_stage2 <- beta0 + beta1 * age_range + beta2 * age_range^2 + gamma0 + gamma1 * age_range + gamma2 * age_range^2 + weibull_baseline
log_relative_hazard_stage3 <- beta0 + beta1 * age_range + beta2 * age_range^2 + lambda0 + lambda1 * age_range + lambda2 * age_range^2 + weibull_baseline
log_relative_hazard_stage4 <- beta0 + beta1 * age_range + beta2 * age_range^2 + weibull_baseline  # baseline stage

#print(log_relative_hazard_stage1)
#print(log_relative_hazard_stage2)
#print(log_relative_hazard_stage3)
#print(log_relative_hazard_stage4)


# Create the plot
plot(age_range, log_relative_hazard_stage1, type = "l", col = "blue", lty = 1,
     ylim = range(log_relative_hazard_stage1, log_relative_hazard_stage2, log_relative_hazard_stage3, log_relative_hazard_stage4),
     xlab = "Age", ylab = "Log Relative Hazard", main = "Log Relative Hazard vs Age by Stage")
lines(age_range, log_relative_hazard_stage2, col = "red", lty = 2)
lines(age_range, log_relative_hazard_stage3, col = "green", lty = 3)
lines(age_range, log_relative_hazard_stage4, col = "purple", lty = 4)

# Add a legend
legend("bottomleft", legend = c("Stage 1", "Stage 2", "Stage 3", "Stage 4 (Baseline)"),
       col = c("blue", "red", "green", "purple"), lty = 1:4)

2c. Generate a partial effect plot of median survival by age, stratified by stage. (Uncertainty interval not required for now.) There is no reason to think the log relative hazard is linear in age.

3. Parametric baseline hazard stratified by stage, a common age relationship.

3a.Write a mathematical expression of the model. The hazard function is defined as:

The log hazard function can be expressed as:

\[ \log(h(t|x)) = \beta_0 + \beta_1 \text{age} + \beta_2 \text{age}^2 + \log(h_0^{(1)}(t; \text{shape}_1, \text{scale}_1)) I_1 + \log(h_0^{(2)}(t; \text{shape}_2, \text{scale}_2)) I_2 + \log(h_0^{(3)}(t; \text{shape}_3, \text{scale}_3)) I_3 + \log(h_0^{(4)}(t; \text{shape}_4, \text{scale}_4)) \]

writeLines("

data {
  int<lower=1> N;  // Number of observations
  vector[N] age;  // Age variable
  array[N] int<lower=0, upper=1> stage2;
  array[N] int<lower=0, upper=1> stage3;
  array[N] int<lower=0, upper=1> stage1;
  array[N] int<lower=1, upper=4> stage; 
  vector[N] time;  // Time until death variable
}

parameters {
  real beta0;                 // Intercept
  real beta1;                 // Coefficient for age
  real beta2;                 // Coefficient for age^2
  real<lower=0> shape1;       // Shape parameter for stage 1
  real<lower=0> scale1;       // Scale parameter for stage 1
  real<lower=0> shape2;       // Shape parameter for stage 2
  real<lower=0> scale2;       // Scale parameter for stage 2
  real<lower=0> shape3;       // Shape parameter for stage 3
  real<lower=0> scale3;       // Scale parameter for stage 3
  real<lower=0> shape4;       // Shape parameter for stage 4
  real<lower=0> scale4;       // Scale parameter for stage 4
}

model {
  for (i in 1:N) {  // 
    real linear_pred = beta0+ beta1* age[i] + beta2 * square(age[i]);
    if (stage[i] == 1){
    target += weibull_lpdf(time[i] | shape1, scale1) +linear_pred + (exp(linear_pred)-1) * weibull_lccdf(time[i] | shape1, scale1); 
    } else if (stage[i] == 2){
    target += weibull_lpdf(time[i] | shape2, scale2) +linear_pred + (exp(linear_pred)-1) * weibull_lccdf(time[i] | shape2, scale2); 
    } else if (stage[i] == 3){
    target += weibull_lpdf(time[i] | shape3, scale3) +linear_pred + (exp(linear_pred)-1) * weibull_lccdf(time[i] | shape3, scale3); 
    } else{
    target += weibull_lpdf(time[i] | shape4, scale4) +linear_pred + (exp(linear_pred)-1) * weibull_lccdf(time[i] | shape4, scale4); 
    }
  }
}
", con = "model3.stan")

# Compile the Stan model
stan_model_3 <- stan_model(file = "model3.stan")

Trying to compile a simple C file

Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
using C compiler: ‘Apple clang version 14.0.3 (clang-1403.0.22.14.1)’
using SDK: ‘MacOSX13.3.sdk’
clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/opt/R/arm64/include    -fPIC  -falign-functions=64 -Wall -g -O2  -c foo.c -o foo.o
In file included from <built-in>:1:
In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
#include <cmath>
         ^~~~~~~
1 error generated.
make: *** [foo.o] Error 1

# Perform MLE using the optimizing function
fit_mle_3 <- optimizing(stan_model_3, data = stan_data, hessian = TRUE)

Warning in .local(object, ...): non-zero return code in optimizing

# View the optimized parameter estimates
print(fit_mle_3)

$par
      beta0       beta1       beta2      shape1      scale1      shape2 
-0.04381929  1.42715097 -0.88984900  0.14087244  0.80079301  3.54921880 
     scale2      shape3      scale3      shape4      scale4 
 2.82453255  1.86674776  1.49953918  2.02842669  1.04906025 

$value
[1] -1045984

$return_code
[1] 70

$hessian
       beta0 beta1 beta2    shape1        scale1    shape2        scale2
beta0      0     0     0   0.00000  0.000000e+00    0.0000  0.000000e+00
beta1      0     0     0   0.00000  0.000000e+00    0.0000  0.000000e+00
beta2      0     0     0   0.00000  0.000000e+00    0.0000  0.000000e+00
shape1     0     0     0  27.67489 -1.620033e+01    0.0000  0.000000e+00
scale1     0     0     0 -16.20033 -3.552714e-12    0.0000  0.000000e+00
shape2     0     0     0   0.00000  0.000000e+00  570.5092 -4.897922e+02
scale2     0     0     0   0.00000  0.000000e+00 -489.7922 -1.261924e-08
shape3     0     0     0   0.00000  0.000000e+00    0.0000  0.000000e+00
scale3     0     0     0   0.00000  0.000000e+00    0.0000  0.000000e+00
shape4     0     0     0   0.00000  0.000000e+00    0.0000  0.000000e+00
scale4     0     0     0   0.00000  0.000000e+00    0.0000  0.000000e+00
          shape3        scale3     shape4        scale4
beta0     0.0000  0.000000e+00    0.00000  0.000000e+00
beta1     0.0000  0.000000e+00    0.00000  0.000000e+00
beta2     0.0000  0.000000e+00    0.00000  0.000000e+00
shape1    0.0000  0.000000e+00    0.00000  0.000000e+00
scale1    0.0000  0.000000e+00    0.00000  0.000000e+00
shape2    0.0000  0.000000e+00    0.00000  0.000000e+00
scale2    0.0000  0.000000e+00    0.00000  0.000000e+00
shape3  180.7504 -3.789498e+02    0.00000  0.000000e+00
scale3 -378.9498  2.273737e-10    0.00000  0.000000e+00
shape4    0.0000  0.000000e+00 -109.06502 -8.925078e+01
scale4    0.0000  0.000000e+00  -89.25078 -7.105427e-12

$theta_tilde
           beta0    beta1     beta2    shape1   scale1   shape2   scale2
[1,] -0.04381929 1.427151 -0.889849 0.1408724 0.800793 3.549219 2.824533
       shape3   scale3   shape4  scale4
[1,] 1.866748 1.499539 2.028427 1.04906

3b.Generate a partial effect plot of log relative hazard by age, stratified by stage. (Uncertainty interval not required for now.)

# Extract MLE estimates for parameters
beta0 <- fit_mle_3$par["beta0"]
beta1 <- fit_mle_3$par["beta1"]
beta2 <- fit_mle_3$par["beta2"]

shape1 <- fit_mle_3$par["shape1"]
scale1 <- fit_mle_3$par["scale1"]
shape2 <- fit_mle_3$par["shape2"]
scale2 <- fit_mle_3$par["scale2"]
shape3 <- fit_mle_3$par["shape3"]
scale3 <- fit_mle_3$par["scale3"]
shape4 <- fit_mle_3$par["shape4"]
scale4 <- fit_mle_3$par["scale4"]


# Define the age range for plotting
age_range <- seq(20, 80, by = 1)  # Modify this as needed

# Weibull baseline hazard for each age point
weibull_baseline1 <- (shape1 - 1) * log(age_range / scale1)
weibull_baseline2 <- (shape2 - 1) * log(age_range / scale2)
weibull_baseline3 <- (shape3 - 1) * log(age_range / scale3)
weibull_baseline4 <- (shape4 - 1) * log(age_range / scale4)

# Calculate log relative hazards for each stage
log_relative_hazard_stage1 <- beta0 + beta1 * age_range + beta2 * age_range^2 + weibull_baseline1
log_relative_hazard_stage2 <- beta0 + beta1 * age_range + beta2 * age_range^2 + weibull_baseline2
log_relative_hazard_stage3 <- beta0 + beta1 * age_range + beta2 * age_range^2 + weibull_baseline3
log_relative_hazard_stage4 <- beta0 + beta1 * age_range + beta2 * age_range^2 + weibull_baseline4  # baseline stage

#print(log_relative_hazard_stage1)
#print(log_relative_hazard_stage2)
#print(log_relative_hazard_stage3)
#print(log_relative_hazard_stage4)


# Create the plot
plot(age_range, log_relative_hazard_stage1, type = "l", col = "blue", lty = 1,
     ylim = range(log_relative_hazard_stage1, log_relative_hazard_stage2, log_relative_hazard_stage3, log_relative_hazard_stage4),
     xlab = "Age", ylab = "Log Relative Hazard", main = "Log Relative Hazard vs Age by Stage")
lines(age_range, log_relative_hazard_stage2, col = "red", lty = 2)
lines(age_range, log_relative_hazard_stage3, col = "green", lty = 3)
lines(age_range, log_relative_hazard_stage4, col = "purple", lty = 4)

# Add a legend
legend("topleft", legend = c("Stage 1", "Stage 2", "Stage 3", "Stage 4 (Baseline)"),
       col = c("blue", "red", "green", "purple"), lty = 1:4)

3c. Generate a partial effect plot of median survival by age, stratified by stage. (Uncertainty interval not required for now.) There is no reason to think the log relative hazard is linear in age.

4.Parametric baseline hazard stratified by stage, a stage-specific age relationship.

4a.Write a mathematical expression of the model.

The log hazard function is: \[ \log(h(t|x)) = \beta_0 + \beta_1 \text{age} + \beta_2 \text{age}^2 + \begin{cases} \alpha_0 + \alpha_1 \text{age} + \alpha_2 \text{age}^2 + \log(h_0^{(1)}(t; \text{shape}_1, \text{scale}_1)), & \text{if stage 1} \\ \gamma_0 + \gamma_1 \text{age} + \gamma_2 \text{age}^2 + \log(h_0^{(2)}(t; \text{shape}_2, \text{scale}_2)), & \text{if stage 2} \\ \lambda_0 + \lambda_1 \text{age} + \lambda_2 \text{age}^2 + \log(h_0^{(3)}(t; \text{shape}_3, \text{scale}_3)), & \text{if stage 3} \\ \log(h_0^{(4)}(t; \text{shape}_4, \text{scale}_4)), & \text{if stage 4} \end{cases} \]

writeLines("

data {
  int<lower=1> N;  // Number of observations
  vector[N] age;  // Age variable
  array[N] int<lower=0, upper=1> stage2;
  array[N] int<lower=0, upper=1> stage3;
  array[N] int<lower=0, upper=1> stage1;
  array[N] int<lower=1, upper=4> stage; 
  vector[N] time;  // Time until death variable
}

parameters {
  // Common age-related coefficients
  real beta0;  // intercept stage 4
  real beta1;  // age coefficient (linear term) stage 4
  real beta2;  // age coefficient (quadratic term) stage 4

  // Stage-specific age relationship coefficients
  real alpha0;  // stage 1 constant
  real alpha1;  // stage 1 age coefficient (linear)
  real alpha2;  // stage 1 age coefficient (quadratic)

  real gamma0;  // stage 2 constant
  real gamma1;  // stage 2 age coefficient (linear)
  real gamma2;  // stage 2 age coefficient (quadratic)

  real lambda0; // stage 3 constant
  real lambda1; // stage 3 age coefficient (linear)
  real lambda2; // stage 3 age coefficient (quadratic)

  // Stage-specific baseline hazard parameters (Weibull)
  real<lower=0> shape1;  // shape for stage 1
  real<lower=0> scale1;  // scale for stage 1

  real<lower=0> shape2;  // shape for stage 2
  real<lower=0> scale2;  // scale for stage 2

  real<lower=0> shape3;  // shape for stage 3
  real<lower=0> scale3;  // scale for stage 3

  real<lower=0> shape4;  // shape for stage 4
  real<lower=0> scale4;  // scale for stage 4
}

model {
  for (i in 1:N) {
    real log_baseline_hazard;
    
    // Common age-related effects
    real age_effect = beta0 + beta1 * age[i] + beta2 * square(age[i]);

    // Stage-specific effects
    if (stage[i] == 1) {
      age_effect += alpha0 + alpha1 * age[i] + alpha2 * square(age[i]);
      log_baseline_hazard = log(shape1) + (shape1 - 1) * log(time[i]) - time[i] / scale1;
    } else if (stage[i] == 2) {
      age_effect += gamma0 + gamma1 * age[i] + gamma2 * square(age[i]);
      log_baseline_hazard = log(shape2) + (shape2 - 1) * log(time[i]) - time[i] / scale2;
    } else if (stage[i] == 3) {
      age_effect += lambda0 + lambda1 * age[i] + lambda2 * square(age[i]);
      log_baseline_hazard = log(shape3) + (shape3 - 1) * log(time[i]) - time[i] / scale3;
    } else {
      log_baseline_hazard = log(shape4) + (shape4 - 1) * log(time[i]) - time[i] / scale4;
    }

    
    target += - exp(log_baseline_hazard + age_effect); 
  }
}
", con = "model4.stan")

# Compile the Stan model
stan_model_4 <- stan_model(file = "model4.stan")

Trying to compile a simple C file

Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
using C compiler: ‘Apple clang version 14.0.3 (clang-1403.0.22.14.1)’
using SDK: ‘MacOSX13.3.sdk’
clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/opt/R/arm64/include    -fPIC  -falign-functions=64 -Wall -g -O2  -c foo.c -o foo.o
In file included from <built-in>:1:
In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
#include <cmath>
         ^~~~~~~
1 error generated.
make: *** [foo.o] Error 1

# Perform MLE using the optimizing function
fit_mle_4 <- optimizing(stan_model_4, data = stan_data, hessian = TRUE)

Warning in .local(object, ...): non-zero return code in optimizing

# View the optimized parameter estimates
print(fit_mle_4)

$par
      beta0       beta1       beta2      alpha0      alpha1      alpha2 
-0.43550266 -0.60070079 -1.41903047 -0.18899425  0.68898917  1.06105235 
     gamma0      gamma1      gamma2     lambda0     lambda1     lambda2 
 1.22282789  0.51766160  0.44307481 -0.02309947  0.78144732 -1.72174805 
     shape1      scale1      shape2      scale2      shape3      scale3 
 0.57103710  0.27660174  0.28790450  2.19587606  1.43097855  1.75980966 
     shape4      scale4 
 6.07637241  0.28994549 

$value
[1] -6.198385e-96

$return_code
[1] 70

$hessian
                 beta0          beta1          beta2        alpha0
beta0    -6.198386e-96  -1.487684e-94  -3.670631e-93 -6.198386e-96
beta1    -1.487684e-94  -3.570612e-93  -8.809926e-92 -1.487684e-94
beta2    -3.670631e-93  -8.809926e-92  -2.172092e-90 -3.670631e-93
alpha0   -6.198386e-96  -1.487684e-94  -3.670631e-93 -6.198386e-96
alpha1   -1.487684e-94  -3.570612e-93  -8.809926e-92 -1.487684e-94
alpha2   -3.670631e-93  -8.809926e-92  -2.172092e-90 -3.670631e-93
gamma0  -1.372886e-227 -3.157916e-226 -7.606064e-225  0.000000e+00
gamma1  -3.157638e-226 -7.263206e-225 -1.749395e-223  0.000000e+00
gamma2  -7.262567e-225 -1.670537e-223 -4.023608e-222  0.000000e+00
lambda0   0.000000e+00   0.000000e+00   0.000000e+00  0.000000e+00
lambda1   0.000000e+00   0.000000e+00   0.000000e+00  0.000000e+00
lambda2   0.000000e+00   0.000000e+00   0.000000e+00  0.000000e+00
shape1   -1.083558e-95  -2.600665e-94  -6.416740e-93 -1.083558e-95
scale1   -8.306312e-95  -1.993610e-93  -4.918926e-92 -8.306312e-95
shape2  -2.019092e-227 -4.644320e-226 -1.118618e-224  0.000000e+00
scale2  -3.253131e-227 -7.482860e-226 -1.802300e-224  0.000000e+00
shape3    0.000000e+00   0.000000e+00   0.000000e+00  0.000000e+00
scale3    0.000000e+00   0.000000e+00   0.000000e+00  0.000000e+00
shape4    0.000000e+00   0.000000e+00   0.000000e+00  0.000000e+00
scale4    0.000000e+00   0.000000e+00   0.000000e+00  0.000000e+00
               alpha1        alpha2         gamma0         gamma1
beta0   -1.487684e-94 -3.670631e-93 -1.372886e-227 -3.157638e-226
beta1   -3.570612e-93 -8.809926e-92 -3.157916e-226 -7.263206e-225
beta2   -8.809926e-92 -2.172092e-90 -7.606064e-225 -1.749395e-223
alpha0  -1.487684e-94 -3.670631e-93   0.000000e+00   0.000000e+00
alpha1  -3.570612e-93 -8.809926e-92   0.000000e+00   0.000000e+00
alpha2  -8.809926e-92 -2.172092e-90   0.000000e+00   0.000000e+00
gamma0   0.000000e+00  0.000000e+00 -2.745772e-227 -6.315554e-226
gamma1   0.000000e+00  0.000000e+00 -6.315554e-226 -1.452641e-224
gamma2   0.000000e+00  0.000000e+00 -1.486863e-224 -3.419932e-223
lambda0  0.000000e+00  0.000000e+00   0.000000e+00   0.000000e+00
lambda1  0.000000e+00  0.000000e+00   0.000000e+00   0.000000e+00
lambda2  0.000000e+00  0.000000e+00   0.000000e+00   0.000000e+00
shape1  -2.600665e-94 -6.416740e-93   0.000000e+00   0.000000e+00
scale1  -1.993610e-93 -4.918926e-92   0.000000e+00   0.000000e+00
shape2   0.000000e+00  0.000000e+00 -4.038185e-227 -9.288234e-226
scale2   0.000000e+00  0.000000e+00 -6.506261e-227 -1.496506e-225
shape3   0.000000e+00  0.000000e+00   0.000000e+00   0.000000e+00
scale3   0.000000e+00  0.000000e+00   0.000000e+00   0.000000e+00
shape4   0.000000e+00  0.000000e+00   0.000000e+00   0.000000e+00
scale4   0.000000e+00  0.000000e+00   0.000000e+00   0.000000e+00
                gamma2 lambda0 lambda1 lambda2        shape1        scale1
beta0   -7.262567e-225       0       0       0 -1.083558e-95 -8.306312e-95
beta1   -1.670537e-223       0       0       0 -2.600665e-94 -1.993610e-93
beta2   -4.023608e-222       0       0       0 -6.416740e-93 -4.918926e-92
alpha0    0.000000e+00       0       0       0 -1.083558e-95 -8.306312e-95
alpha1    0.000000e+00       0       0       0 -2.600665e-94 -1.993610e-93
alpha2    0.000000e+00       0       0       0 -6.416740e-93 -4.918926e-92
gamma0  -1.486863e-224       0       0       0  0.000000e+00  0.000000e+00
gamma1  -3.419932e-223       0       0       0  0.000000e+00  0.000000e+00
gamma2  -8.047216e-222       0       0       0  0.000000e+00  0.000000e+00
lambda0   0.000000e+00       0       0       0  0.000000e+00  0.000000e+00
lambda1   0.000000e+00       0       0       0  0.000000e+00  0.000000e+00
lambda2   0.000000e+00       0       0       0  0.000000e+00  0.000000e+00
shape1    0.000000e+00       0       0       0 -2.357922e-95 -1.452051e-94
scale1    0.000000e+00       0       0       0 -1.452051e-94 -1.030041e-93
shape2  -2.186718e-224       0       0       0  0.000000e+00  0.000000e+00
scale2  -3.523206e-224       0       0       0  0.000000e+00  0.000000e+00
shape3    0.000000e+00       0       0       0  0.000000e+00  0.000000e+00
scale3    0.000000e+00       0       0       0  0.000000e+00  0.000000e+00
shape4    0.000000e+00       0       0       0  0.000000e+00  0.000000e+00
scale4    0.000000e+00       0       0       0  0.000000e+00  0.000000e+00
                shape2         scale2 shape3 scale3 shape4 scale4
beta0   -2.019092e-227 -3.253131e-227      0      0      0      0
beta1   -4.644320e-226 -7.482860e-226      0      0      0      0
beta2   -1.118618e-224 -1.802300e-224      0      0      0      0
alpha0    0.000000e+00   0.000000e+00      0      0      0      0
alpha1    0.000000e+00   0.000000e+00      0      0      0      0
alpha2    0.000000e+00   0.000000e+00      0      0      0      0
gamma0  -4.038185e-227 -6.506261e-227      0      0      0      0
gamma1  -9.288234e-226 -1.496506e-225      0      0      0      0
gamma2  -2.186718e-224 -3.523206e-224      0      0      0      0
lambda0   0.000000e+00   0.000000e+00      0      0      0      0
lambda1   0.000000e+00   0.000000e+00      0      0      0      0
lambda2   0.000000e+00   0.000000e+00      0      0      0      0
shape1    0.000000e+00   0.000000e+00      0      0      0      0
scale1    0.000000e+00   0.000000e+00      0      0      0      0
shape2  -7.237409e-227 -9.624607e-227      0      0      0      0
scale2  -9.624607e-227 -9.425647e-227      0      0      0      0
shape3    0.000000e+00   0.000000e+00      0      0      0      0
scale3    0.000000e+00   0.000000e+00      0      0      0      0
shape4    0.000000e+00   0.000000e+00      0      0      0      0
scale4    0.000000e+00   0.000000e+00      0      0      0      0

$theta_tilde
          beta0      beta1    beta2     alpha0    alpha1   alpha2   gamma0
[1,] -0.4355027 -0.6007008 -1.41903 -0.1889942 0.6889892 1.061052 1.222828
        gamma1    gamma2     lambda0   lambda1   lambda2    shape1    scale1
[1,] 0.5176616 0.4430748 -0.02309947 0.7814473 -1.721748 0.5710371 0.2766017
        shape2   scale2   shape3  scale3   shape4    scale4
[1,] 0.2879045 2.195876 1.430979 1.75981 6.076372 0.2899455

4b.Generate a partial effect plot of log relative hazard by age, stratified by stage. (Uncertainty interval not required for now.)

4c. Generate a partial effect plot of median survival by age, stratified by stage. (Uncertainty interval not required for now.) There is no reason to think the log relative hazard is linear in age.

5.Mean of log days with stage-specific age relationship.

Mean log(time) is directly related to the expected survival time on a logarithmic scale, rather than the instantaneous rate of event occurrence. The coefficients in this model represent direct adjustments to the expected log survival time based on age and stage, without needing a baseline hazard term to describe time-dependent risk.

5a.Write a mathematical expression of the model.

\[ E(\log(\text{time})) = \beta_0 + \beta_1 (\text{age}) + \beta_2 (\text{age}^2) + \alpha_0 I_1 + \alpha_1 I_1 (\text{age}) + \alpha_2 I_1 (\text{age}^2) + \gamma \]

writeLines("
 
data {
   int<lower=1> N;  // Number of observations
   vector[N] age;  // Age variable
   array[N] int<lower=0, upper=1> stage2;
   array[N] int<lower=0, upper=1> stage3;
   array[N] int<lower=0, upper=1> stage1;
   vector[N] time;  // Time until death variable
 }
 
 parameters {
   real beta0; // constant coefficient stage 4
   real beta1;  // Coefficient for age stage 4
   real beta2;  // Coefficient for squared age stage 4
   real gamma0; //coefficient for stage 2 constant
   real gamma1; //cofefficient for stage 2 age
   real gamma2; // coefficient for stage 2 age squared
   real lambda0; //coefficient for stage 3 constant
   real lambda1; //cofefficient for stage 3 age
   real lambda2; // coefficient for stage 3 age squared
   real alpha0; //coefficient for stage 1 constant
   real alpha1; //cofefficient for stage 1 age
   real alpha2; // coefficient for stage 1 age squared
 }
 

 
 model {
   for (i in 1:N) {  // Loop through each observation
     // Calculate the linear predictor based on age and stage
     real linear_pred = beta0 + beta1 * age[i] + beta2 * square(age[i])
                         + alpha0  * stage1[i]  + alpha1 * stage1[i] * age[i] + alpha2 * stage1[i] * square(age[i])
                         + gamma0 * stage2[i]  + gamma1 *  stage2[i] * age[i] + gamma2  * stage2[i] * square(age[i])
                         +  lambda0  *stage3[i]  + lambda1 *  stage3[i] * age[i] + lambda2  * stage3[i] * square(age[i])
                         ;
     
     target += linear_pred ;  
   }
 }
 ", con = "model5.stan")

# Compile the Stan model
stan_model_5 <- stan_model(file = "model5.stan")

Trying to compile a simple C file

Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
using C compiler: ‘Apple clang version 14.0.3 (clang-1403.0.22.14.1)’
using SDK: ‘MacOSX13.3.sdk’
clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/opt/R/arm64/include    -fPIC  -falign-functions=64 -Wall -g -O2  -c foo.c -o foo.o
In file included from <built-in>:1:
In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
#include <cmath>
         ^~~~~~~
1 error generated.
make: *** [foo.o] Error 1

# Perform MLE using the optimizing function
fit_mle_5 <- optimizing(stan_model_5, data = stan_data, hessian = TRUE)

Warning in .local(object, ...): non-zero return code in optimizing

# View the optimized parameter estimates
print(fit_mle_5)

$par
      beta0       beta1       beta2      gamma0      gamma1      gamma2 
 1.25138855  1.26567211  1.17775996  1.18613198  1.67223558 -1.93708778 
    lambda0     lambda1     lambda2      alpha0      alpha1      alpha2 
-0.35606700  0.18850521  1.52807883 -1.66984174  1.34602946 -0.09825555 

$value
[1] 1533943

$return_code
[1] 70

$hessian
        beta0 beta1 beta2 gamma0 gamma1 gamma2 lambda0 lambda1 lambda2 alpha0
beta0       0     0     0      0      0      0       0       0       0      0
beta1       0     0     0      0      0      0       0       0       0      0
beta2       0     0     0      0      0      0       0       0       0      0
gamma0      0     0     0      0      0      0       0       0       0      0
gamma1      0     0     0      0      0      0       0       0       0      0
gamma2      0     0     0      0      0      0       0       0       0      0
lambda0     0     0     0      0      0      0       0       0       0      0
lambda1     0     0     0      0      0      0       0       0       0      0
lambda2     0     0     0      0      0      0       0       0       0      0
alpha0      0     0     0      0      0      0       0       0       0      0
alpha1      0     0     0      0      0      0       0       0       0      0
alpha2      0     0     0      0      0      0       0       0       0      0
        alpha1 alpha2
beta0        0      0
beta1        0      0
beta2        0      0
gamma0       0      0
gamma1       0      0
gamma2       0      0
lambda0      0      0
lambda1      0      0
lambda2      0      0
alpha0       0      0
alpha1       0      0
alpha2       0      0

$theta_tilde
        beta0    beta1   beta2   gamma0   gamma1    gamma2   lambda0   lambda1
[1,] 1.251389 1.265672 1.17776 1.186132 1.672236 -1.937088 -0.356067 0.1885052
      lambda2    alpha0   alpha1      alpha2
[1,] 1.528079 -1.669842 1.346029 -0.09825555

5b.Generate a partial effect plot of log relative hazard by age, stratified by stage. (Uncertainty interval not required for now.)

# Extract MLE estimates for parameters
beta0 <- fit_mle_5$par["beta0"]
beta1 <- fit_mle_5$par["beta1"]
beta2 <- fit_mle_5$par["beta2"]

gamma0 <- fit_mle_5$par["gamma0"]
gamma1 <- fit_mle_5$par["gamma1"]
gamma2 <- fit_mle_5$par["gamma2"]
lambda0 <- fit_mle_5$par["lambda0"]
lambda1 <- fit_mle_5$par["lambda1"]
lambda2 <- fit_mle_5$par["lambda2"]
alpha0 <- fit_mle_5$par["alpha0"]
alpha1 <- fit_mle_5$par["alpha1"]
alpha2 <- fit_mle_5$par["alpha2"]




# Define the age range for plotting
age_range <- seq(20, 80, by = 1)  # Modify this as needed



# Calculate log relative hazards for each stage
log_relative_hazard_stage1 <- beta0 + beta1 * age_range + beta2 * age_range^2 + alpha0 + alpha1 * age_range + alpha2 * age_range^2 
log_relative_hazard_stage2 <- beta0 + beta1 * age_range + beta2 * age_range^2 + gamma0 + gamma1 * age_range + gamma2 * age_range^2 
log_relative_hazard_stage3 <- beta0 + beta1 * age_range + beta2 * age_range^2 + lambda0 + lambda1 * age_range + lambda2 * age_range^2 
log_relative_hazard_stage4 <- beta0 + beta1 * age_range + beta2 * age_range^2   # baseline stage

#print(log_relative_hazard_stage1)
#print(log_relative_hazard_stage2)
#print(log_relative_hazard_stage3)
#print(log_relative_hazard_stage4)


# Create the plot
plot(age_range, log_relative_hazard_stage1, type = "l", col = "blue", lty = 1,
     ylim = range(log_relative_hazard_stage1, log_relative_hazard_stage2, log_relative_hazard_stage3, log_relative_hazard_stage4),
     xlab = "Age", ylab = "Log Relative Hazard", main = "Log Relative Hazard vs Age by Stage (Model 5)")
lines(age_range, log_relative_hazard_stage2, col = "red", lty = 2)
lines(age_range, log_relative_hazard_stage3, col = "green", lty = 3)
lines(age_range, log_relative_hazard_stage4, col = "purple", lty = 4)

# Add a legend
legend("topleft", legend = c("Stage 1", "Stage 2", "Stage 3", "Stage 4 (Baseline)"),
       col = c("blue", "red", "green", "purple"), lty = 1:4)

5c. Generate a partial effect plot of median survival by age, stratified by stage. (Uncertainty interval not required for now.) There is no reason to think the log relative hazard is linear in age.