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 
-0.5598518 -1.0624169  1.7672216  0.3559118  0.8064101 -0.9885981  1.3195376 
     scale 
 0.4167581 

$value
[1] -1179411

$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 1372.2885 -6.597689e+02
scale       0      0      0     0     0     0 -659.7689  1.364242e-09

$theta_tilde
         alpha1    alpha2   alpha3     beta0     beta1      beta2    shape
[1,] -0.5598518 -1.062417 1.767222 0.3559118 0.8064101 -0.9885981 1.319538
         scale
[1,] 0.4167581

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(-3,7), #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)

# View the optimized parameter estimates
print(fit_mle_2)

$par
        beta0         beta1         beta2        gamma0        gamma1 
 0.6464515570  0.0323388154 -0.0006104756 -0.4194350423 -0.2247135911 
       gamma2       lambda0       lambda1       lambda2        alpha0 
 0.0025081998 -1.8701123189 -0.1255456542  0.0020705542  0.0942775457 
       alpha1        alpha2         shape         scale 
-0.2375350214  0.0028685194  1.8575548932  0.9531061340 

$value
[1] -1028.952

$return_code
[1] 0

$hessian
                beta0         beta1         beta2        gamma0        gamma1
beta0       -510.4419 -2.358671e+04      -9446214     -139.0758 -6.565968e+03
beta1     -23586.7145 -1.214945e+06    -714451701    -6565.9683 -3.382448e+05
beta2   -9446213.9882 -7.144517e+08 -106814489718 -3046610.9881 -2.384079e+08
gamma0      -139.0758 -6.565968e+03      -3046611     -139.0758 -6.565968e+03
gamma1     -6565.9683 -3.382448e+05    -238407919    -6565.9683 -3.382448e+05
gamma2  -3046610.9881 -2.384079e+08  -37150322646 -3046610.9881 -2.384079e+08
lambda0     -211.6907 -9.412945e+03      -3862417        0.0000  0.000000e+00
lambda1    -9412.9454 -4.815500e+05    -288453468        0.0000  0.000000e+00
lambda2 -3862416.8336 -2.884535e+08  -42331029599        0.0000  0.000000e+00
alpha0      -117.5396 -5.442339e+03      -2141634        0.0000  0.000000e+00
alpha1     -5442.3387 -2.788804e+05    -162242757        0.0000  0.000000e+00
alpha2  -2141633.7236 -1.622428e+08  -24336927125        0.0000  0.000000e+00
shape      -1641.0333 -6.998503e+04     -16971940     -605.1684 -2.828857e+04
scale        948.1741  4.381363e+04      17546861      258.3409  1.219665e+04
              gamma2       lambda0       lambda1      lambda2        alpha0
beta0       -3046611     -211.6907 -9.412945e+03     -3862417     -117.5396
beta1     -238407919    -9412.9454 -4.815500e+05   -288453468    -5442.3387
beta2   -37150322646 -3862416.8336 -2.884535e+08 -42331029599 -2141633.7236
gamma0      -3046611        0.0000  0.000000e+00            0        0.0000
gamma1    -238407919        0.0000  0.000000e+00            0        0.0000
gamma2  -37150322646        0.0000  0.000000e+00            0        0.0000
lambda0            0     -211.6907 -9.412945e+03     -3862417        0.0000
lambda1            0    -9412.9454 -4.815500e+05   -288453468        0.0000
lambda2            0 -3862416.8336 -2.884535e+08 -42331029599        0.0000
alpha0             0        0.0000  0.000000e+00            0     -117.5396
alpha1             0        0.0000  0.000000e+00            0    -5442.3387
alpha2             0        0.0000  0.000000e+00            0 -2141633.7236
shape       -9484990     -547.0479 -1.996054e+04     -3515303     -479.0589
scale        5659247      393.2272  1.748507e+04      7174651      218.3363
               alpha1       alpha2         shape         scale
beta0   -5.442339e+03     -2141634 -1.641033e+03      948.1741
beta1   -2.788804e+05   -162242757 -6.998503e+04    43813.6259
beta2   -1.622428e+08 -24336927125 -1.697194e+07 17546861.4762
gamma0   0.000000e+00            0 -6.051684e+02      258.3409
gamma1   0.000000e+00            0 -2.828857e+04    12196.6490
gamma2   0.000000e+00            0 -9.484990e+06  5659247.2766
lambda0  0.000000e+00            0 -5.470479e+02      393.2272
lambda1  0.000000e+00            0 -1.996054e+04    17485.0663
lambda2  0.000000e+00            0 -3.515303e+06  7174651.4712
alpha0  -5.442339e+03     -2141634 -4.790589e+02      218.3363
alpha1  -2.788804e+05   -162242757 -2.102181e+04    10109.4449
alpha2  -1.622428e+08 -24336927125 -3.855065e+06  3978202.3083
shape   -2.102181e+04     -3855065 -7.262819e+03     3067.7118
scale    1.010944e+04      3978202  3.067712e+03    -1761.2858

$theta_tilde
         beta0      beta1         beta2    gamma0     gamma1    gamma2
[1,] 0.6464516 0.03233882 -0.0006104756 -0.419435 -0.2247136 0.0025082
       lambda0    lambda1     lambda2     alpha0    alpha1      alpha2    shape
[1,] -1.870112 -0.1255457 0.002070554 0.09427755 -0.237535 0.002868519 1.857555
         scale
[1,] 0.9531061

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     scale2 
 1.2343700 -0.3302439 -0.5081484  0.5943482  0.5140426  3.3881508  0.3146353 
    shape3     scale3     shape4     scale4 
 0.1580331  5.3887654  0.7483702  0.7487255 

$value
[1] -623066

$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 147.06114 -6.835004e+01    0.0000  0.000000e+00
scale1     0     0     0 -68.35004 -4.263256e-11    0.0000  0.000000e+00
shape2     0     0     0   0.00000  0.000000e+00 1570.7759 -4.675648e+02
scale2     0     0     0   0.00000  0.000000e+00 -467.5648 -8.062671e-07
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.00000  0.000000e+00   0.00000  0.000000e+00
beta1    0.00000  0.000000e+00   0.00000  0.000000e+00
beta2    0.00000  0.000000e+00   0.00000  0.000000e+00
shape1   0.00000  0.000000e+00   0.00000  0.000000e+00
scale1   0.00000  0.000000e+00   0.00000  0.000000e+00
shape2   0.00000  0.000000e+00   0.00000  0.000000e+00
scale2   0.00000  0.000000e+00   0.00000  0.000000e+00
shape3 -25.73455 -3.208072e+01   0.00000  0.000000e+00
scale3 -32.08072 -7.105427e-12   0.00000  0.000000e+00
shape4   0.00000  0.000000e+00 -29.13261 -3.292829e+01
scale4   0.00000  0.000000e+00 -32.92829  1.421085e-11

$theta_tilde
       beta0      beta1      beta2    shape1    scale1   shape2    scale2
[1,] 1.23437 -0.3302439 -0.5081484 0.5943482 0.5140426 3.388151 0.3146353
        shape3   scale3    shape4    scale4
[1,] 0.1580331 5.388765 0.7483702 0.7487255

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 
 1.092611525  1.342161113 -0.280058214  1.088995584 -1.772265674  0.216876137 
      gamma0       gamma1       gamma2      lambda0      lambda1      lambda2 
 1.979560089 -0.240936469 -0.629026436 -1.064573197  1.851909447  0.008795181 
      shape1       scale1       shape2       scale2       shape3       scale3 
 0.563326913  0.171467570  6.198661294  2.127730793  0.698974223  5.584424166 
      shape4       scale4 
 1.463150043  0.981991639 

$value
[1] -9.956879e-24

$return_code
[1] 70

$hessian
                 beta0          beta1          beta2        alpha0
beta0    -9.956880e-24  -2.091029e-22  -4.462880e-21 -9.540265e-29
beta1    -2.091029e-22  -4.391341e-21  -9.372438e-20 -2.558227e-27
beta2    -4.462880e-21  -9.372438e-20  -1.999846e-18 -7.170322e-26
alpha0   -9.540265e-29  -2.558227e-27  -7.170322e-26 -9.540265e-29
alpha1   -2.558227e-27  -6.864912e-26  -1.925559e-24 -2.558227e-27
alpha2   -7.170322e-26  -1.925559e-24  -5.395195e-23 -7.170322e-26
gamma0  -1.205978e-193 -2.773993e-192 -6.681360e-191  0.000000e+00
gamma1  -2.773749e-192 -6.380185e-191 -1.536713e-189  0.000000e+00
gamma2  -6.379623e-191 -1.467443e-189 -3.534440e-188  0.000000e+00
lambda0  -9.956785e-24  -2.091003e-22  -4.462808e-21  0.000000e+00
lambda1  -2.091003e-22  -4.391272e-21  -9.372245e-20  0.000000e+00
lambda2  -4.462808e-21  -9.372245e-20  -1.999792e-18  0.000000e+00
shape1   -1.179416e-28  -3.153677e-27  -8.813715e-26 -1.179416e-28
scale1   -8.764038e-28  -2.327950e-26  -6.461558e-25 -8.764038e-28
shape2  -1.502584e-192 -3.456248e-191 -8.324618e-190  0.000000e+00
scale2  -3.654936e-193 -8.407092e-192 -2.024908e-190  0.000000e+00
shape3   -1.829149e-23  -3.841358e-22  -8.198576e-21  0.000000e+00
scale3   -6.652973e-24  -1.397179e-22  -2.981991e-21  0.000000e+00
shape4    1.925830e-51   4.429799e-50   1.066949e-48  0.000000e+00
scale4   -8.654004e-54  -1.990596e-52  -4.794492e-51  0.000000e+00
               alpha1        alpha2         gamma0         gamma1
beta0   -2.558227e-27 -7.170322e-26 -1.205978e-193 -2.773749e-192
beta1   -6.864912e-26 -1.925559e-24 -2.773993e-192 -6.380185e-191
beta2   -1.925559e-24 -5.395195e-23 -6.681360e-191 -1.536713e-189
alpha0  -2.558227e-27 -7.170322e-26   0.000000e+00   0.000000e+00
alpha1  -6.864912e-26 -1.925559e-24   0.000000e+00   0.000000e+00
alpha2  -1.925559e-24 -5.395195e-23   0.000000e+00   0.000000e+00
gamma0   0.000000e+00  0.000000e+00 -2.411956e-193 -5.547743e-192
gamma1   0.000000e+00  0.000000e+00 -5.547743e-192 -1.276037e-190
gamma2   0.000000e+00  0.000000e+00 -1.306098e-190 -3.004155e-189
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  -3.153677e-27 -8.813715e-26   0.000000e+00   0.000000e+00
scale1  -2.327950e-26 -6.461558e-25   0.000000e+00   0.000000e+00
shape2   0.000000e+00  0.000000e+00 -3.005217e-192 -6.912303e-191
scale2   0.000000e+00  0.000000e+00 -7.309871e-193 -1.681344e-191
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
beta0   -6.379623e-191 -9.956785e-24 -2.091003e-22 -4.462808e-21 -1.179416e-28
beta1   -1.467443e-189 -2.091003e-22 -4.391272e-21 -9.372245e-20 -3.153677e-27
beta2   -3.534440e-188 -4.462808e-21 -9.372245e-20 -1.999792e-18 -8.813715e-26
alpha0    0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00 -1.179416e-28
alpha1    0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00 -3.153677e-27
alpha2    0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00 -8.813715e-26
gamma0  -1.306098e-190  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
gamma1  -3.004155e-189  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
gamma2  -7.068879e-188  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
lambda0   0.000000e+00 -9.956785e-24 -2.091003e-22 -4.462808e-21  0.000000e+00
lambda1   0.000000e+00 -2.091003e-22 -4.391272e-21 -9.372245e-20  0.000000e+00
lambda2   0.000000e+00 -4.462808e-21 -9.372245e-20 -1.999792e-18  0.000000e+00
shape1    0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00 -1.699379e-28
scale1    0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00 -1.122934e-27
shape2  -1.627355e-189  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
scale2  -3.958369e-190  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
shape3    0.000000e+00 -1.829149e-23 -3.841358e-22 -8.198576e-21  0.000000e+00
scale3    0.000000e+00 -6.652973e-24 -1.397179e-22 -2.981991e-21  0.000000e+00
shape4    0.000000e+00  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  0.000000e+00
               scale1         shape2         scale2        shape3        scale3
beta0   -8.764038e-28 -1.502584e-192 -3.654936e-193 -1.829149e-23 -6.652973e-24
beta1   -2.327950e-26 -3.456248e-191 -8.407092e-192 -3.841358e-22 -1.397179e-22
beta2   -6.461558e-25 -8.324618e-190 -2.024908e-190 -8.198576e-21 -2.981991e-21
alpha0  -8.764038e-28   0.000000e+00   0.000000e+00  0.000000e+00  0.000000e+00
alpha1  -2.327950e-26   0.000000e+00   0.000000e+00  0.000000e+00  0.000000e+00
alpha2  -6.461558e-25   0.000000e+00   0.000000e+00  0.000000e+00  0.000000e+00
gamma0   0.000000e+00 -3.005217e-192 -7.309871e-193  0.000000e+00  0.000000e+00
gamma1   0.000000e+00 -6.912303e-191 -1.681344e-191  0.000000e+00  0.000000e+00
gamma2   0.000000e+00 -1.627355e-189 -3.958369e-190  0.000000e+00  0.000000e+00
lambda0  0.000000e+00   0.000000e+00   0.000000e+00 -1.829149e-23 -6.652973e-24
lambda1  0.000000e+00   0.000000e+00   0.000000e+00 -3.841358e-22 -1.397179e-22
lambda2  0.000000e+00   0.000000e+00   0.000000e+00 -8.198576e-21 -2.981991e-21
shape1  -1.122934e-27   0.000000e+00   0.000000e+00  0.000000e+00  0.000000e+00
scale1  -8.152567e-27   0.000000e+00   0.000000e+00  0.000000e+00  0.000000e+00
shape2   0.000000e+00 -4.050819e-191 -9.240305e-192  0.000000e+00  0.000000e+00
scale2   0.000000e+00 -9.240305e-192 -1.542884e-192  0.000000e+00  0.000000e+00
shape3   0.000000e+00   0.000000e+00   0.000000e+00 -4.298028e-23 -1.337460e-23
scale3   0.000000e+00   0.000000e+00   0.000000e+00 -1.337460e-23  8.782929e-25
shape4   0.000000e+00   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  0.000000e+00
               shape4        scale4
beta0    1.925830e-51 -8.654004e-54
beta1    4.429799e-50 -1.990596e-52
beta2    1.066949e-48 -4.794492e-51
alpha0   0.000000e+00  0.000000e+00
alpha1   0.000000e+00  0.000000e+00
alpha2   0.000000e+00  0.000000e+00
gamma0   0.000000e+00  0.000000e+00
gamma1   0.000000e+00  0.000000e+00
gamma2   0.000000e+00  0.000000e+00
lambda0  0.000000e+00  0.000000e+00
lambda1  0.000000e+00  0.000000e+00
lambda2  0.000000e+00  0.000000e+00
shape1   0.000000e+00  0.000000e+00
scale1   0.000000e+00  0.000000e+00
shape2   0.000000e+00  0.000000e+00
scale2   0.000000e+00  0.000000e+00
shape3   0.000000e+00  0.000000e+00
scale3   0.000000e+00  0.000000e+00
shape4  -1.334180e-50  8.096156e-53
scale4   8.096156e-53  1.694419e-53

$theta_tilde
        beta0    beta1      beta2   alpha0    alpha1    alpha2  gamma0
[1,] 1.092612 1.342161 -0.2800582 1.088996 -1.772266 0.2168761 1.97956
         gamma1     gamma2   lambda0  lambda1     lambda2    shape1    scale1
[1,] -0.2409365 -0.6290264 -1.064573 1.851909 0.008795181 0.5633269 0.1714676
       shape2   scale2    shape3   scale3  shape4    scale4
[1,] 6.198661 2.127731 0.6989742 5.584424 1.46315 0.9819916

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 
 0.29716757 -1.51245917  1.24009246 -1.26041857 -1.75606663  1.08907075 
    lambda0     lambda1     lambda2      alpha0      alpha1      alpha2 
 0.01878918  0.17918058  0.87162355  0.95851644 -1.33304918 -1.92364958 

$value
[1] 1702797

$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
[1,] 0.2971676 -1.512459 1.240092 -1.260419 -1.756067 1.089071 0.01878918
       lambda1   lambda2    alpha0    alpha1   alpha2
[1,] 0.1791806 0.8716235 0.9585164 -1.333049 -1.92365

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.