INLA

Contents

• Quarto
• Running Code
• Load required libraries

Quarto

Quarto enables you to weave together content and executable code into a finished document. To learn more about Quarto see https://quarto.org.

Running Code

When you click the Render button a document will be generated that includes both content and the output of embedded code. You can embed code like this:

library(mvtnorm)
library(MASS)
library(pracma)

You can add options to executable code like this

Attaching package: 'Matrix'

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

    expm, lu, tril, triu

Load required libraries

library(INLA)

Loading required package: sp

This is INLA_24.06.27 built 2024-06-27 03:00:51 UTC.
 - See www.r-inla.org/contact-us for how to get help.
 - List available models/likelihoods/etc with inla.list.models()
 - Use inla.doc(<NAME>) to access documentation

library(INLA)
library(ggplot2)
library(reshape2)


Y=y

# Assuming Y is your 3D array of binary outcomes: patients x diseases x time
# N: number of patients, D: number of diseases, T: number of time points, K: number of factors

# 1. Data Preparation
N <- dim(Y)[1]
D <- dim(Y)[2]
T <- dim(Y)[3]
K <- 3  # Number of factors, adjust as needed

data_long <- melt(Y)
colnames(data_long) <- c("patient", "disease", "time", "y")
data_long$patient_factor <- paste0(data_long$patient, "_", data_long$time)
data_long$disease_factor <- paste0(data_long$disease, "_", data_long$time)


# 2. Construct precision matrices
construct_precision_matrix <- function(times, sigma, l) {
  n <- length(times)
  D <- as.matrix(dist(times))
  K <- sigma^2 * exp(-D^2 / (2*l^2))
  Q <- solve(K + diag(1e-6, n))
  return(Q)
}

unique_times <- 1:T
Time_K <- construct_precision_matrix(unique_times, 1, T/4)

# 3. Define INLA formula
formula <- y ~ f(disease, model="iid") +
  f(patient_factor, model="iid", group=patient, control.group=list(model="ar1")) +
  f(disease_factor, model="iid", group=disease, control.group=list(model="ar1", order=K))

# 4. Fit INLA model
result <- inla(formula,
               family = "binomial",
               data = data_long,
               control.predictor = list(compute = TRUE),
               control.compute = list(dic = TRUE, waic = TRUE))

# 5. Extract results
mu_d <- result$summary.random$disease$mean

# Extract disease factor loadings (φ_kd(t))
phi_posterior <- result$summary.random$disease_factor
phi_matrix <- array(phi_posterior$mean, dim=c(D, K, T))

# Extract patient-specific factors (λ_ik(t))
lambda_posterior <- result$summary.random$patient_factor
lambda_matrix <- array(lambda_posterior$mean, dim=c(N, K, T))

# 6. Compute time-varying disease correlations
disease_cor <- array(NA, dim=c(D, D, T))
for(t in 1:T) {
  disease_cor[,,t] <- cor(t(phi_matrix[,,t]))
}

# 7. Visualization functions
plot_disease_cor <- function(t) {
  cor_data <- melt(disease_cor[,,t])
  ggplot(cor_data, aes(Var1, Var2, fill=value)) +
    geom_tile() +
    scale_fill_gradient2(low="blue", high="red", mid="white", midpoint=0, limits=c(-1,1)) +
    theme_minimal() +
    ggtitle(paste("Disease Correlations at Time", t))
}

plot_factor_loadings <- function(k) {
  factor_data <- melt(phi_matrix[,k,])
  colnames(factor_data) <- c("Disease", "Time", "Loading")
  ggplot(factor_data, aes(x=Time, y=Loading, color=factor(Disease))) +
    geom_line() +
    theme_minimal() +
    ggtitle(paste("Disease Loadings for Factor", k))
}

# 8. Generate visualizations
# Plot disease correlations at a specific time point
print(plot_disease_cor(1))

# Plot factor loadings for each factor
for(k in 1:K) {
  print(plot_factor_loadings(k))
}

# Create an animation of disease correlations over time

  for(t in 1:T) {
    print(plot_disease_cor(t))
  }

# 9. Model evaluation
print(result$dic$dic)

[1] 2279.436

print(result$waic$waic)

[1] 2272.068

# 10. Print summary of key parameters
print(summary(result))

Time used:
    Pre = 0.871, Running = 74.9, Post = 6.46, Total = 82.2 
Fixed effects:
              mean    sd 0.025quant 0.5quant 0.975quant   mode kld
(Intercept) -3.989 0.123     -4.241   -3.985     -3.756 -3.985   0

Random effects:
  Name    Model
    disease IID model
   patient_factor IID model
   disease_factor IID model

Model hyperparameters:
                                  mean       sd 0.025quant 0.5quant 0.975quant
Precision for disease        22072.946 2.51e+04   1414.365 1.42e+04   8.86e+04
Precision for patient_factor 22088.845 2.46e+04   1450.207 1.44e+04   8.75e+04
GroupRho for patient_factor     -0.001 6.99e-01     -0.988 0.00e+00   9.87e-01
Precision for disease_factor     1.375 3.19e-01      0.852 1.34e+00   2.10e+00
GroupRho for disease_factor      0.003 6.99e-01     -0.988 7.00e-03   9.88e-01
                                 mode
Precision for disease        3820.573
Precision for patient_factor 3938.580
GroupRho for patient_factor    -0.998
Precision for disease_factor    1.274
GroupRho for disease_factor    -0.998

Deviance Information Criterion (DIC) ...............: 2279.44
Deviance Information Criterion (DIC, saturated) ....: 2264.97
Effective number of parameters .....................: 59.53

Watanabe-Akaike information criterion (WAIC) ...: 2272.07
Effective number of parameters .................: 51.74

Marginal log-Likelihood:  -1168.32 
 is computed 
Posterior summaries for the linear predictor and the fitted values are computed
(Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')