<jupyter_start><jupyter_text>Replication - High Dimensional Case2Here we provide a notebook to replicate the simulation results for the high dimensional case2 simulation. As replicating all of the runs is rather time consuming, we provide this script to replicate a specified simulation run, which can then compared to the data provided on github.The notebook replicates the results in:- /out/simulation/sim_hd2The main script can be found at: - /scripts/simulation/highdimensional_case2.py- /scripts/simulation/highdimensional_case2.R  Please choose the settup for replication:<jupyter_code>run_id = 13 # any integer 1-200
suffix = 'rank5' # choose: rank5, rank50
Rsuffix = 'R_Alasso1_theta' # 'R_lasso_theta', 'R_lasso_theta_1se', 'R_Alasso1_theta', 'R_Alasso1_theta_1se', 'R_Alasso2_theta', 'R_Alasso2_theta_1se, c, R_SCAD_theta, R_MCP_theta'
# google colab specific - installing probcox
!pip3 install torch==1.7.0
!pip3 install pyro-ppl==1.5.1
!pip3 install probcox
# load R 
%load_ext rpy2.ipython
# install packages
%%R 
install.packages('survival')
install.packages('glmnet')
install.packages('coxed')
install.packages('ncvreg')
install.packages('doMC')
install.packages('BVSNLP')

# Modules
# =======================================================================================================================
import os
import sys
import shutil
import subprocess
import tqdm

import numpy as np
import pandas as pd

import torch
from torch.distributions import constraints

import pyro
import pyro.distributions as dist

from pyro.infer import SVI, Trace_ELBO

import matplotlib.pyplot as plt

import warnings
warnings.filterwarnings("ignore")

import probcox as pcox

dtype = torch.FloatTensor

np.random.seed(890)
torch.manual_seed(543)
%%R -i run_id -o simdata

# Generate Data
set.seed(13)
library(survival)
library(glmnet)
library(doMC)
library(coxed)
library(ncvreg)
library(BVSNLP)

# generate paramters - 10 non -zero effects
theta <- c(-0.5, 0.7, 1.2, 0.65, -0.9, 1.4, 0.2, -0.4, -1.3, 0.1)
theta <- c(theta, replicate(990, 0))

# generate datasets
for(ii in seq(200)){
  simdata <- sim.survdata(N=750, T=10000, xvars=1000, censor=0.75, num.data.frames=1, beta=theta)$data
  if(ii==run_id){break}
}

X = torch.from_numpy(np.asarray(simdata.iloc[:, :-2]).astype(float)).type(dtype)
surv = torch.from_numpy(np.concatenate((np.zeros((750, 1)), simdata.iloc[:, -2:]), axis=1).astype(float)).type(dtype)

total_obs = surv.shape[0]
total_events = torch.sum(surv[:, -1] == 1).numpy().tolist()

# Inference Setup
# =======================================================================================================================
# Custom linear predictor - Here: simple linear combination
def predictor(data):
    theta =  pyro.sample("theta", dist.StudentT(1, loc=0, scale=0.001).expand([data[1].shape[1], 1])).type(dtype)
    pred = torch.mm(data[1], theta)
    return(pred)

def evaluate(surv, X, rank, batchsize, sampling_proportion, iter_, predictor=predictor):
    sampling_proportion[1] = batchsize
    eta=5 # paramter for optimization
    run = True # repeat initalization if NAN encounterd while training - gauge correct optimization settings
    while run:
        run = False
        pyro.clear_param_store()
        m = pcox.PCox(sampling_proportion=sampling_proportion, predictor=predictor)
        m.initialize(eta=eta, rank=rank, num_particles=5)
        loss=[0]
        for ii in tqdm.tqdm(range((iter_))):
            idx = np.unique(np.concatenate((np.random.choice(np.where(surv[:, -1]==1)[0], 1, replace=False), np.random.choice(range(surv.shape[0]), batchsize-1, replace=False)))) # random sample of data - force at least one event (no evaluation otherwise)
            data=[surv[idx], X[idx]] # subsampled data
            loss.append(m.infer(data=data))
            # divergence check
            if loss[-1] != loss[-1]:
                eta = eta * 0.1
                run=True
                break
    g = m.return_guide()
    out = g.quantiles([0.025, 0.5, 0.975])
    return(out)
# Run
# =======================================================================================================================
if suffix == 'rank5':
    pyro.clear_param_store()
    out = evaluate(rank=5, batchsize=256, iter_=15000, surv=surv, X=X, sampling_proportion=[total_obs, None, total_events, None])
elif suffix == 'rank50':
    pyro.clear_param_store()
    out = evaluate(rank=50, batchsize=256, iter_=15000, surv=surv, X=X, sampling_proportion=[total_obs, None, total_events, None])


# Download data from github
# =======================================================================================================================

url1 = 'https://raw.githubusercontent.com/alexwjung/ProbCox/main/paper/ProbCox/out/simulation/sim_hd2/probcox' + suffix +'_theta.txt'
url2 = 'https://raw.githubusercontent.com/alexwjung/ProbCox/main/paper/ProbCox/out/simulation/sim_hd2/probcox' + suffix +'_theta_lower.txt'
url3 = 'https://raw.githubusercontent.com/alexwjung/ProbCox/main/paper/ProbCox/out/simulation/sim_hd2/probcox' + suffix +'_theta_upper.txt'


theta_est_sim = pd.read_csv(url1, header=None, sep=';')
theta_est_lower_sim = pd.read_csv(url2, header=None, sep=';')
theta_est_upper_sim = pd.read_csv(url3, header=None, sep=';')<jupyter_output><empty_output><jupyter_text>The estimates can vary to a certain degree - this is due to a numerical instability we observe for the prior distribution. Choosing a scale of lower then 0.01 for the studentT distribution makes it very difficult to get exact replications. However, the overall results are still valid and the replication provides very similar estimates.<jupyter_code># Comparision
# =======================================================================================================================

print('Mean Absolute Error - estmiate: ', np.mean(np.abs(np.squeeze(out['theta'][1].detach().numpy()) - np.squeeze(np.asarray(theta_est_sim.iloc[np.where(theta_est_sim.iloc[:, 0] == run_id)[0], 1:-1])))))

print('Mean Absolute Error - lower bound estmiate: ',np.mean(np.abs(np.squeeze(out['theta'][0].detach().numpy()) - np.squeeze(np.asarray(theta_est_lower_sim.iloc[np.where(theta_est_lower_sim.iloc[:, 0] == run_id)[0], 1:-1])))))

print('Mean Absolute Error - upper bound estmiate: ',np.mean(np.abs(np.squeeze(out['theta'][2].detach().numpy()) - np.squeeze(np.asarray(theta_est_upper_sim.iloc[np.where(theta_est_upper_sim.iloc[:, 0] == run_id)[0], 1:-1])))))

# Comparision - Plot
# =======================================================================================================================
theta = np.concatenate((np.asarray([-0.5, 0.7, 1.2, 0.65, -0.9, 1.4, 0.2, -0.4, -1.3, 0.1]), np.zeros((990,))))[:, None]
theta_est_sim = np.squeeze(np.asarray(theta_est_sim.iloc[np.where(theta_est_sim.iloc[:, 0] == run_id)[0], 1:-1]))[:, None]
theta_est_lower_sim = np.squeeze(np.asarray(theta_est_lower_sim.iloc[np.where(theta_est_lower_sim.iloc[:, 0] == run_id)[0], 1:-1]))[:, None]
theta_est_upper_sim = np.squeeze(np.asarray(theta_est_upper_sim.iloc[np.where(theta_est_upper_sim.iloc[:, 0] == run_id)[0], 1:-1]))[:, None]

theta_est = out['theta'][1].detach().numpy()
theta_est_lower = out['theta'][0].detach().numpy()
theta_est_upper = out['theta'][2].detach().numpy()

fig, ax = plt.subplots(1, 1, figsize=(10, 10), dpi=100)  
ax.errorbar(theta[:, 0], theta_est[:, 0], yerr=(theta_est[:, 0] - theta_est_lower[:, 0], theta_est_upper[:, 0]- theta_est[:, 0]),  ls='', c="red", capsize=2, capthick=0.95, elinewidth=0.95, label='colab')

ax.plot(theta[:, 0], theta_est[:, 0], ls='', c=".3", marker='x', ms=2)
ax.errorbar(theta[:, 0]+0.025, theta_est_sim[:, 0], yerr=(theta_est_sim[:, 0] - theta_est_lower_sim[:, 0], theta_est_upper_sim[:, 0]- theta_est_sim[:, 0]),  ls='', c="#5c161e", capsize=2, capthick=0.95, elinewidth=0.95, label='paper')
ax.plot(theta[:, 0]+0.025, theta_est_sim[:, 0], ls='', c=".3", marker='x', ms=2)

ax.set(xlim=(-2, 2), ylim=(-2, 2))
ax.plot(ax.get_xlim(), ax.get_ylim(), ls="--", c="black", linewidth=0.75)
ax.set_yticks([-2, 0, 2])
ax.set_ylim([-2, 2])
ax.set_xticks([-2, 0, 2])
ax.set_xlim([-2, 2])
ax.set_xlabel('theta')
ax.set_ylabel('theta_hat')
plt.legend()
plt.show()
plt.close()<jupyter_output><empty_output><jupyter_text>R evaluation<jupyter_code>%%R -o R_lasso -o R_lasso_1se -o R_Alasso1 -o R_Alasso1_1se -o R_Alasso2 -o R_Alasso2_1se -o R_SCAD -o R_MCP -o R_BVSNLP

sim <- simdata
yss = Surv(simdata$y, simdata$failed)

# ridge
# Lasso 
cv.fit <- c()
cv.fit <- cv.glmnet(as.matrix(sim[, 1:1000]), yss, family ="cox", nfolds=3, parallel=TRUE, type.measure ="C", alpha=1)
R_lasso <- as.matrix(unname(coef(cv.fit, s=cv.fit$lambda.min)))
R_lasso_1se <- as.matrix(unname(coef(cv.fit, s=cv.fit$lambda.1se)))

w1 <- 1/abs(as.numeric(coef(cv.fit, s=cv.fit$lambda.min)))
w1[w1 == Inf] <- 999999999 

w2 <- 1/abs(as.numeric(coef(cv.fit, s=cv.fit$lambda.1se)))
w2[w2 == Inf] <- 999999999 

# Ridge 
#cv.fit <- c()
#cv.fit <- cv.glmnet(as.matrix(sim[, 1:1000]), yss, family ="cox", nfolds=3, parallel=TRUE, type.measure ="C", alpha=0)
#x = paste(sim_name, paste(unname(coef(cv.fit, s=cv.fit$lambda.min)), collapse="; "), sep='; ')
#write(x, file = paste(ROOT_DIR, '/out/simulation/sim_hd2/R_ridge_theta.txt', sep=''), ncolumns = 1, append = TRUE, sep = ";")
#x = paste(sim_name, paste(unname(coef(cv.fit, s=cv.fit$lambda.1se)), collapse="; "), sep='; ')
#write(x, file = paste(ROOT_DIR, '/out/simulation/sim_hd2/R_ridge_theta_1se.txt', sep=''), ncolumns = 1, append = TRUE, sep = ";")

#w3 <- 1/abs(as.numeric(coef(cv.fit, s=cv.fit$lambda.min)))
#w3[w3 == Inf] <- 999999999 

#w4 <- 1/abs(as.numeric(coef(cv.fit, s=cv.fit$lambda.1se)))
#w4[w4 == Inf] <- 999999999 


# Adaptive Lasso
# Alasso 
cv.fit <- c()
cv.fit <-cv.glmnet(as.matrix(sim[, 1:1000]), yss, family ="cox", nfolds=3, parallel=TRUE, type.measure="C", penalty.factor=w1)
R_Alasso1 <- as.matrix(unname(coef(cv.fit, s=cv.fit$lambda.min)))
R_Alasso1_1se <- as.matrix(unname(coef(cv.fit, s=cv.fit$lambda.1se)))

# Alasso 2
cv.fit <- c()
cv.fit <-cv.glmnet(as.matrix(sim[, 1:1000]), yss, family ="cox", nfolds=3, parallel=TRUE, type.measure="C", penalty.factor=w2)
R_Alasso2 <- as.matrix(unname(coef(cv.fit, s=cv.fit$lambda.min)))
R_Alasso2_1se <- as.matrix(unname(coef(cv.fit, s=cv.fit$lambda.1se)))

# Alasso 3
#cv.fit <- c()
#cv.fit <-cv.glmnet(as.matrix(sim[, 1:1000]), yss, family ="cox", nfolds=3, parallel=TRUE, type.measure="C", penalty.factor=w3)
#x = paste(sim_name, paste(unname(coef(cv.fit, s=cv.fit$lambda.min)), collapse="; "), sep='; ')
#write(x, file = paste(ROOT_DIR, '/out/simulation/sim_hd2/R_Alasso3_theta.txt', sep=''), ncolumns = 1, append = TRUE, sep = ";")
#x = paste(sim_name, paste(unname(coef(cv.fit, s=cv.fit$lambda.1se)), collapse="; "), sep='; ')
#write(x, file = paste(ROOT_DIR, '/out/simulation/sim_hd2/R_Alasso3_theta_1se.txt', sep=''), ncolumns = 1, append = TRUE, sep = ";")

# Alasso 4
#cv.fit <- c()
#cv.fit <-cv.glmnet(as.matrix(sim[, 1:1000]), yss, family ="cox", nfolds=3, parallel=TRUE, type.measure="C", penalty.factor=w4)
#x = paste(sim_name, paste(unname(coef(cv.fit, s=cv.fit$lambda.min)), collapse="; "), sep='; ')
#write(x, file = paste(ROOT_DIR, '/out/simulation/sim_hd2/R_Alasso4_theta.txt', sep=''), ncolumns = 1, append = TRUE, sep = ";")
#x = paste(sim_name, paste(unname(coef(cv.fit, s=cv.fit$lambda.1se)), collapse="; "), sep='; ')
#write(x, file = paste(ROOT_DIR, '/out/simulation/sim_hd2/R_Alasso4_theta_1se.txt', sep=''), ncolumns = 1, append = TRUE, sep = ";")

# SCAD
m <- cv.ncvsurv(as.matrix(sim[, 1:1000]), yss, penalty=c("SCAD"), nfolds=3, se=c('bootstrap'), trace=TRUE)
fit <- m$fit$beta[,m$min]
R_SCAD <- as.matrix(unname(m$fit$beta[,m$min]))

# MCP
m <- cv.ncvsurv(as.matrix(sim[, 1:1000]), yss, penalty=c("MCP"), nfolds=3, se=c('bootstrap'), trace=TRUE)
fit <- m$fit$beta[,m$min]
R_MCP <- as.matrix(unname(m$fit$beta[,m$min]))

# BVSNLP
bout <- bvs(sim[, 1:1000], yss, family = "survival", nlptype = "piMOM", niter = 2000, ncpu=8)
coef <- rep(0, 1000)
coef[bout$HPM] <- bout$beta_hat
R_BVSNLP <- as.matrix(unname(coef))




# Download data from github
# =======================================================================================================================
# choose any R model to load for comparison (need to select corresponding run id.)
url = 'https://raw.githubusercontent.com/alexwjung/ProbCox/main/paper/ProbCox/out/simulation/sim_hd2/' + Rsuffix + '.txt'
theta_R = pd.read_csv(url, header=None, sep=';')


# R_lasso R_lasso_1se R_Alasso1 R_Alasso1_1se R_Alasso2  R_Alasso2 R_SCAD R_MCP choose from the corresponding
R_colab = R_Alasso1 #R_Alasso1_theta
# Comparision
# =======================================================================================================================
print('Mean Absolute Error: ', np.mean(np.abs(R_colab[:, 0] - np.squeeze(np.asarray(theta_R.iloc[np.where(theta_R.iloc[:, 0] == run_id)[0], 1:])))))

plt.scatter(R_colab[:10, 0], np.squeeze(np.asarray(theta_R.iloc[np.where(theta_R.iloc[:, 0] == run_id)[0], 1:]))[:10])<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.24732334047109208, 'fraction_numerical': 0.03647904557968798, 'mean_word_length': 5.699282786885246, 'pattern_counts': {'":': 0, '<': 53, '<?xml version=': 0, '>': 13, 'https://': 4, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 4, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '364429', 'license': 'MIT', 'n_tokens_mistral': 5289, 'n_tokens_neox': 4866, 'n_words': 950, 'path': 'paper/ProbCox/replication/simulations/highdimensional_case2.ipynb', 'repo_name': 'alexwjung/ProbCox', 'size': '13076'}