Chapter 5 Generate sample GP

library(ggplot2)
library(kernlab)

# Placeholder for GP Predictions and Original Densities
gppreds <- list()
original_densities <- list()

# Sample data (assuming you have 12 diseases for simplicity)
set.seed(42)  # For reproducibility
sampled_disease <- sample(unique(m$diag_icd10), 12)
sampled_pheno <- sample(unique(m$phenotype), 12)

# Define a common grid for x values
common_grid <- seq(30, max(m$age_diag), length.out = 100)

# Loop over the diseases
for(i in 1:length(sampled_disease)){
  # Select data for the specific disease
  disease_data <- m[m$diag_icd10 %in% sampled_disease[i],]
  
  # Calculate density
  density_data <- density(disease_data$age_diag)
  
  # Add a small noise term to the y-values
  #density_data$y <- density_data$y + rnorm(length(density_data$y), mean = 0, sd = 0.001)
  
  # Fit Gaussian Process
  gp_fit <- gausspr(x = density_data$x, y = density_data$y, kernel = "rbfdot", kpar = "automatic")
  
  # Predict using the GP on the common grid
  y_pred <- predict(gp_fit, common_grid)
  
  # Truncate predictions to avoid inflation in unobserved regions
  y_pred[common_grid < min(disease_data$age_diag)] <- 0
  y_pred[common_grid > max(disease_data$age_diag)] <- 0
  
  # Store original density and predictions
  original_densities[[i]] <- data.frame(x = density_data$x, y = density_data$y, phenotype = sampled_pheno[i])
  gppreds[[i]] <- data.frame(x = common_grid, y = y_pred, phenotype = sampled_pheno[i])
}

## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel

# Combine all original densities and predictions into data frames
original_density_df <- do.call(rbind, original_densities)
gp_predictions <- do.call(rbind, gppreds)

# Plot using ggplot2
ggplot() +
  geom_line(data = original_density_df, aes(x = x, y = y, color = phenotype), linetype = "dashed") +
  geom_line(data = gp_predictions, aes(x = x, y = y, color = phenotype)) +
  labs(title = "Density Estimation with Gaussian Process",
       x = "Age at Diagnosis",
       y = "Density") +
  theme_minimal() +
  theme(legend.position = "bottom")

library(ggplot2)
library(splines)
library(kernlab)

# Placeholder for Spline Predictions and Original Densities
spline_preds <- list()
original_densities <- list()
gppreds <- list()

# Sample data (assuming you have 12 diseases for simplicity)
set.seed(42)  # For reproducibility
sampled_disease <- sample(unique(m$diag_icd10), 12)
sampled_pheno <- sample(unique(m$phenotype), 12)

# Define a common grid for x values
common_grid <- seq(30, max(m$age_diag), length.out = 100)

# Loop over the diseases
for(i in 1:length(sampled_disease)){
  # Select data for the specific disease
  disease_data <- m[m$diag_icd10 %in% sampled_disease[i],]
  
  # Calculate density
  density_data <- density(disease_data$age_diag)
  
  # Fit Spline
  spline_fit <- lm(y ~ ns(x, df = 5), data = data.frame(x = density_data$x, y = density_data$y))
  
  # Predict using the spline on the common grid
  y_pred_spline <- predict(spline_fit, newdata = data.frame(x = common_grid))
  
  # Fit Gaussian Process
  gp_fit <- gausspr(x = density_data$x, y = density_data$y, kernel = "rbfdot", kpar = "automatic")
  
  # Predict using the GP on the common grid
  y_pred_gp <- predict(gp_fit, common_grid)
  
  # Truncate predictions to avoid inflation in unobserved regions
  y_pred_gp[common_grid < min(disease_data$age_diag)] <- 0
  y_pred_gp[common_grid > max(disease_data$age_diag)] <- 0
  y_pred_spline[common_grid < min(disease_data$age_diag)] <- 0
  y_pred_spline[common_grid > max(disease_data$age_diag)] <- 0
  
  # Store original density and predictions
  original_densities[[i]] <- data.frame(x = density_data$x, y = density_data$y, phenotype = sampled_pheno[i], type = "Original")
  gppreds[[i]] <- data.frame(x = common_grid, y = y_pred_gp, phenotype = sampled_pheno[i], type = "GP")
  spline_preds[[i]] <- data.frame(x = common_grid, y = y_pred_spline, phenotype = sampled_pheno[i], type = "Spline")
}

## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel

# Combine all original densities, GP predictions, and spline predictions into data frames
original_density_df <- do.call(rbind, original_densities)
gp_predictions <- do.call(rbind, gppreds)
spline_predictions <- do.call(rbind, spline_preds)

# Combine all data frames for plotting
all_data <- rbind(original_density_df, gp_predictions, spline_predictions)

# Plot using ggplot2
ggplot(all_data, aes(x = x, y = y, color = phenotype, linetype = type)) +
  geom_line() +
  labs(title = "Aladyn Estimation with Gaussian Process and Splines",
       x = "Age at Diagnosis",
       y = "Density",
       linetype = "Model Type") +
  theme_minimal() +
  theme(legend.position = "bottom") +
  scale_linetype_manual(values = c("Original" = "dashed", "GP" = "solid", "Spline" = "dotted"))

5.1 Sample from the Gaussian Process

So now since we have an ability to simulate GP …

set.seed(42)  # For reproducibility

top_scorers=m%>%group_by(diag_icd10)%>%summarise(l=length(diag_icd10))%>%arrange(desc(l))%>%head(K*V)


sampled_disease <- sample(unique(top_scorers$diag_icd10), K*V)
sdm=matrix(sampled_disease,nrow=K,ncol=V)

sampled_pheno <- sample(unique(m$phenotype), K*V)

# Define a common grid for x values
common_grid <- seq(30, max(m$age_diag), length.out = 100)
beta=array(NA,dim=c(K,V,T))


# Loop over the diseases
for(k in 1:K)
  for(d in 1:V){
    # Select data for the specific disease
    
    code=sdm[k,d]
    disease_data <- m[m$diag_icd10 %in% code,]
  
  # Calculate density
    density_data <- density(disease_data$age_diag)
  
   # Fit Gaussian Process
  gp_fit <- gausspr(x = density_data$x, y = density_data$y, kernel = "rbfdot", kpar = "automatic")
  
  # Predict using the GP on the common grid
  y_pred_gp <- predict(gp_fit, common_grid)
  
  # Truncate predictions to avoid inflation in unobserved regions
  #y_pred_gp[common_grid < min(disease_data$age_diag)] <- density_data$y[1]
  #y_pred_gp[common_grid > max(disease_data$age_diag)] <- tail(density_data$y,1)

  beta[k,d,] <- pmax(y_pred_gp,0)/2
  

  }

## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel

num_topics_per_disease = rep(1,V)
# Initialize the list to store the primary topics for each disease
disease_primary_topics <- vector("list", length = V)

# Initialize the disease-topic matrix with zeros indicating inactivity
disease_topic_matrix <- matrix(0, nrow = V, ncol = K)

# Assign primary topics to each disease
for (disease in 1:V) {
  primary_topics <- sample(1:K, num_topics_per_disease[disease])
  disease_primary_topics[[disease]] <- primary_topics
  
  # Set the bias value for the disease in its primary topics, can vary
  #bias_values <- runif(num_topics_per_disease[disease], 0, 1)  # Variable bias between 0 and 1
  bias_values = rep(1,num_topics_per_disease[disease])  # Fixed bias of 1
  low_value=rep(0.1,num_topics_per_disease[disease])
  # Fill in the disease-topic matrix with bias values for the primary topics
  disease_topic_matrix[disease, ]=low_value
  disease_topic_matrix[disease, primary_topics] <- bias_values
  
}


image(t(disease_topic_matrix),xlab="Topic",ylab="Disease")

for(v in 1:V){
  for(k in 1:K){
    for(t in 1:T){
      beta[k,v,t]=beta[k,v,t]*disease_topic_matrix[v,k]}}}

matplot(common_grid,t(beta[,1,]),type ="l",main="Topic Specificty")

Ok now we need to do the warping:

# Simulate genetic scores for each individual (assuming 5 genetic scores)
genetic_scores <- matrix(runif(D * 5, min = 0, max = 1), nrow = D, ncol = 5)

# Define the weights for each genetic score (to be learned)
w <- runif(5, min = 0.5, max = 2)

# Compute the warping coefficients for each individual and topic
rho <- genetic_scores %*% w
rho <- matrix(rho, nrow = D, ncol = K, byrow = TRUE)

time_warping <- function(t, rho, T) {
return(max(1,floor((t / T) ^ (1 / rho) * T)))
}


# Initialize inverse warped time array
warped_t_array <- array(0, dim = c(D, K,V, T))

warped_beta <- array(NA, dim = c(D, K, V, T))

# Apply time warping to each individual's beta values
for (i in 1:D) {
  for (k in 1:K) {
    for (v in 1:V) {
      for (t in 1:T) {
        warped_time <- time_warping(t, rho[i, k], T)
        warped_t_array [i, k, v, t] =warped_time
        warped_beta[i, k, v, t] <- beta[k, v, warped_time]
      }
    }
  }
}

5.2 Show that individuals at high genetic risk have early disease:

sd=sample(V,1)
sk=which.max(disease_topic_matrix[sd,])

hi=which.max(rho[,sk])
low=which.min(rho[,sk])
mid=order(rho[,sk])[D/2]

plot(warped_beta[hi,sk,sd,],col="darkred",main=paste("Disease Probabilities for ",sampled_pheno[4]," for high, mid and low genetic risk"))
points(warped_beta[mid,sk,sd,],col="orange")
points(warped_beta[low,sk,sd,],col="darkgreen")

beta_prime=warped_beta