### Install R Packages for BKMR Analysis

Source: https://jenfb.github.io/bkmr/SimData1

Installs necessary R packages for BKMR analysis, including `bkmr`, `readxl`, and `ggplot2`. These packages are required for loading data, fitting the BKMR model, and generating visualizations.

```R
install.packages("bkmr")
install.packages("readxl")
install.packages("ggplot2")
```

--------------------------------

### Extracting Posterior Inclusion Probabilities (PIPs) from BKMR Fits

Source: https://jenfb.github.io/bkmr/overview

This snippet shows how to extract and compare Posterior Inclusion Probabilities (PIPs) from BKMR models fitted with different variable selection strategies. `ExtractPIPs` function is used to get PIPs for component-wise selection and group-wise/conditional PIPs for hierarchical selection. This helps in model interpretation and identifying important predictors.

```R
ExtractPIPs(fitkm_corr)
```

```R
ExtractPIPs(fitkm_hier)
```

--------------------------------

### R: Compare Computation Time of BKMR Fits

Source: https://jenfb.github.io/bkmr/SimData1

This R code compares the computation time per iteration for different BKMR model fits: the full model, and models using 100 and 50 knots. It calculates the time difference between model start and end, divided by the number of iterations, to provide a per-iteration time. The result is presented in a data frame.

```r
time_per_iter <- data.frame(
    full = with(fitkm, 
                difftime(time2, time1, units = "secs")/iter),
    knots100 = with(fitkm_knots100, 
                    difftime(time2, time1, units = "secs")/iter),
    knots50 = with(fitkm_knots50, 
                   difftime(time2, time1, units = "secs")/iter)
)
time_per_iter
```

--------------------------------

### Load `bkmr` R package

Source: https://jenfb.github.io/bkmr/overview

Loads the `bkmr` R package into the current R session. This package is required for all subsequent operations related to Bayesian kernel machine regression.

```R
library(bkmr)
```

--------------------------------

### Download and Load Dataset for BKMR

Source: https://jenfb.github.io/bkmr/SimData1

Checks if 'DataSet1.xls' exists in the working directory. If not, it downloads the dataset from the specified URL. It then loads the dataset using `readxl::read_excel` into a data frame named `dat`.

```R
if (!"DataSet1.xls" %in% list.files(getwd())) {
    download.file("http://www.niehs.nih.gov/about/events/pastmtg/2015/statistical/dataset1xls.xls", "DataSet1.xls")
}
dat <- readxl::read_excel("DataSet1.xls", sheet = 1, col_names =  TRUE)
```

--------------------------------

### Inspecting BKMR Algorithm Tuning Parameters

Source: https://jenfb.github.io/bkmr/overview

This snippet demonstrates how to view the tuning parameters used in the Markov Chain Monte Carlo (MCMC) fitting process of a BKMR model. The `control.params` argument of the `kmbayes` function stores these parameters, which can be accessed and inspected, often by converting them to a data frame. These parameters control the proposal distributions for Metropolis-Hastings steps and influence MCMC mixing and convergence.

```R
data.frame(fitkm$control.params)
```

--------------------------------

### Fit BKMR Model with Variable Selection

Source: https://jenfb.github.io/bkmr/overview

Fits a Bayesian kernel machine regression (BKMR) model using the `kmbayes` function. Includes MCMC sampling for parameter estimation and performs variable selection on predictors.

```R
set.seed(111)
fitkm <- kmbayes(y = y, Z = Z, X = X, iter = 10000, verbose = FALSE, varsel = TRUE)
```

--------------------------------

### Visualize MCMC Trace Plot for Beta Parameters

Source: https://jenfb.github.io/bkmr/overview

Generates a trace plot for the 'beta' parameters of the fitted BKMR model. This plot helps in visually assessing the convergence of the MCMC sampler for the covariate coefficients.

```R
TracePlot(fit = fitkm, par = "beta")
```

--------------------------------

### Visualize MCMC Trace Plot for Smoothness Parameters

Source: https://jenfb.github.io/bkmr/overview

Generates a trace plot for the 'r' parameter(s) (smoothness tuning parameters) for a specified component of the fitted BKMR model. This plot helps assess MCMC convergence for these parameters.

```R
TracePlot(fit = fitkm, par = "r", comp = 1)
```

--------------------------------

### Investigate Model Convergence with Trace Plots

Source: https://jenfb.github.io/bkmr/SimData1

Generates trace plots to visually inspect the convergence of MCMC samplers for different model parameters. Trace plots display how parameter values change over the course of the sampler's execution. Parameters like 'beta', 'sigsq.eps', and 'r' can be plotted.

```r
TracePlot(fit = fitkm, par = "beta")
```

```r
TracePlot(fit = fitkm, par = "sigsq.eps")
```

```r
TracePlot(fit = fitkm, par = "r", comp = 7)
```

--------------------------------

### Plot Overall Risk Summaries using ggplot

Source: https://jenfb.github.io/bkmr/overview

Visualizes overall risk summaries by plotting quantiles against estimated values, with confidence intervals. Requires the 'ggplot2' package and a data frame 'risks.overall' with 'quantile', 'est', and 'sd' columns.

```r
ggplot(risks.overall, aes(quantile, est, ymin = est - 1.96*sd, ymax = est + 1.96*sd)) +
    geom_pointrange()
```

--------------------------------

### Simulate Data for BKMR Analysis

Source: https://jenfb.github.io/bkmr/overview

Generates simulated data for Bayesian kernel machine regression. This function creates outcome variable `y`, exposure matrix `Z`, and covariate matrix `X`, used for demonstrating model fitting and analysis.

```R
set.seed(111)
dat <- SimData(n = 50, M = 4)
y <- dat$y
Z <- dat$Z
X <- dat$X
```

--------------------------------

### Visualize True Exposure-Response Function

Source: https://jenfb.github.io/bkmr/overview

Generates and displays a 3D perspective plot of the true exposure-response function used in data simulation. This helps in understanding the underlying relationship being modeled.

```R
z1 <- seq(min(dat$Z[, 1]), max(dat$Z[, 1]), length = 20)
z2 <- seq(min(dat$Z[, 2]), max(dat$Z[, 2]), length = 20)
hgrid.true <- outer(z1, z2, function(x,y) apply(cbind(x,y), 1, dat$HFun))

res <- persp(z1, z2, hgrid.true, theta = 30, phi = 20, expand = 0.5, 
             col = "lightblue", xlab = "", ylab = "", zlab = "")
```

--------------------------------

### Load R Packages for BKMR and Plotting

Source: https://jenfb.github.io/bkmr/SimData1

Loads the `bkmr` and `ggplot2` R packages into the current R session. `bkmr` is used for fitting the Bayesian Kernel Machine Regression model, and `ggplot2` is used for data visualization.

```R
library(bkmr)
library(ggplot2)
```

--------------------------------

### Prepare Data Matrices for BKMR Model

Source: https://jenfb.github.io/bkmr/SimData1

Prepares the data for the BKMR model by creating matrices for covariates (`covar`) and exposures (`expos`), and extracting the outcome variable `Y`. Assumes `dat` is a data frame containing these variables.

```R
covar <- data.matrix(dat[, "Z"])
expos <- data.matrix(dat[, c("X1", "X2", "X3", "X4", "X5", "X6", "X7")])
Y <- dat$Y
```

--------------------------------

### Visualize MCMC Trace Plot for Residual Variance

Source: https://jenfb.github.io/bkmr/overview

Generates a trace plot for the 'sigsq.eps' parameter (residual variance) of the fitted BKMR model. This plot aids in assessing the MCMC convergence for the error variance.

```R
TracePlot(fit = fitkm, par = "sigsq.eps")
```

--------------------------------

### Display Overall Cumulative Risk Summaries in R

Source: https://jenfb.github.io/bkmr/SimData1

Prints the calculated overall cumulative risk summaries to the console. The output table shows the quantile, estimated effect, and standard deviation for each specified quantile.

```r
risks.overall
```

--------------------------------

### Compare Methods for Estimating h(z) in BKMR

Source: https://jenfb.github.io/bkmr/overview

Compares three different methods for estimating the h(z) function (the expected or posterior mean of the response given exposures) using a fitted BKMR model. It calculates the posterior mean and standard deviation for each method, including an approximate method and two exact methods. Requires median exposure values (med_vals) and the original data (Z, dat).

```r
med_vals <- apply(Z, 2, median)
Znew <- matrix(med_vals, nrow = 1)
h_true <- dat$HFun(Znew)

h_est1 <- ComputePostmeanHnew(fitkm, Znew = Znew, method = "approx")
h_est2 <- ComputePostmeanHnew(fitkm, Znew = Znew, method = "exact")
set.seed(111)
samps3 <- SamplePred(fitkm, Znew = Znew, Xnew = cbind(0))

h_est_compare <- data.frame(
  method = c("truth", 1:3),
  post_mean = c(h_true, h_est1$postmean, h_est2$postmean, mean(samps3)),
  post_sd = c(NA, sqrt(h_est1$postvar), sqrt(h_est2$postvar), sd(samps3))
)
h_est_compare
```

--------------------------------

### Fit BKMR Model using kmbayes

Source: https://jenfb.github.io/bkmr/SimData1

Fits a Bayesian Kernel Machine Regression (BKMR) model using the `kmbayes` function, which employs a Markov chain Monte Carlo (MCMC) algorithm. Key arguments include the number of MCMC iterations (`iter`), verbosity of output (`verbose`), and whether to perform variable selection (`varsel`).

```r
set.seed(1000)
fitkm <- kmbayes(Y, Z = scale_expos, X = covar, iter = 10000, 
                 verbose = FALSE, varsel = TRUE)
```

--------------------------------

### R: Fit BKMR with Different Numbers of Knots

Source: https://jenfb.github.io/bkmr/SimData1

This R code demonstrates how to fit a BKMR model using a Gaussian predictive process approach with varying numbers of knots (100 and 50). It utilizes the `fields::cover.design` function to generate knot locations and `kmbayes` for model fitting. The `set.seed` function is used for reproducibility.

```r
set.seed(1000)
knots100 <- fields::cover.design(scale_expos, nd = 100)$design
fitkm_knots100 <- kmbayes(Y, scale_expos, covar, iter = 10000, 
                          verbose = FALSE, varsel = TRUE,
                          knots = knots100)

set.seed(1000)
knots50 <- fields::cover.design(scale_expos, nd = 50)$design
fitkm_knots50 <- kmbayes(Y, scale_expos, covar, iter = 10000, 
                             verbose = FALSE, varsel = TRUE,
                             knots = knots50)
```

--------------------------------

### Simulate Correlated Data for BKMR Analysis

Source: https://jenfb.github.io/bkmr/overview

Generates simulated data for BKMR analysis where exposure variables (Z) are correlated. Uses the 'SimData' function with 'Zgen = "corr"' to introduce correlation. The output includes a correlation matrix of the simulated predictors.

```r
set.seed(111)
d2 <- SimData(n = 100, M = 4, Zgen = "corr", sigsq.true = 2.2)
round(cor(d2$Z), 2)
```

--------------------------------

### Compute Single-Predictor Interaction Summaries

Source: https://jenfb.github.io/bkmr/overview

Computes specific 'interaction' parameters by comparing single-predictor health risks under different fixed quantile scenarios for other predictors. Uses the 'SingVarIntSummaries' function, requiring a fitted model ('fitkm'), response ('y'), predictor matrices ('Z', 'X'), and specified quantile differences ('qs.diff', 'qs.fixed').

```r
risks.int <- SingVarIntSummaries(fit = fitkm, y = y, Z = Z, X = X, 
                                 qs.diff = c(0.25, 0.75), 
                                 qs.fixed = c(0.25, 0.75),
                                 method = "exact")
risks.int
```

--------------------------------

### Prepare Bivariate Cross-section Exposure Pairs

Source: https://jenfb.github.io/bkmr/SimData1

Generates a data frame of exposure pairs for visualizing bivariate exposure-response functions. This involves creating combinations of two exposure variables, ensuring that each pair is unique and that the first exposure index is less than the second.

```r
expos.pairs <- subset(data.frame(expand.grid(expos1 = c(1,2,4,5,7), expos2 = c(1,2,4,5,7))), expos1 < expos2)
```

--------------------------------

### Plot Single-Predictor Health Risks using ggplot

Source: https://jenfb.github.io/bkmr/overview

Visualizes the single-predictor health risks, showing the estimated risk and its confidence interval for each predictor at different fixed quantile levels. Uses 'ggplot2' and requires 'risks.singvar' data frame. The plot is flipped for better readability of predictor names.

```r
ggplot(risks.singvar, aes(variable, est, ymin = est - 1.96*sd, 
                          ymax = est + 1.96*sd, col = q.fixed)) +
    geom_pointrange(position = position_dodge(width = 0.75)) +
  coord_flip()
```

--------------------------------

### Visualize Single Exposure Risk Summaries (R)

Source: https://jenfb.github.io/bkmr/SimData1

Generates a plot to visualize the single-exposure risk summaries calculated by `SingVarRiskSummaries`. It uses `ggplot2` to create a point range plot, showing the estimated risk and its confidence interval for different variables and fixed exposure quantiles. The plot is designed to be easily interpretable for identifying the contribution of individual exposures.

```r
ggplot(risks.singvar, aes(variable, est, ymin = est - 1.96*sd, ymax = est + 1.96*sd, 
                          col = q.fixed)) + 
    geom_pointrange(position = position_dodge(width = 0.75)) + coord_flip()
```

--------------------------------

### Plot Exposure-Response at Fixed Quantiles in R

Source: https://jenfb.github.io/bkmr/SimData1

Visualizes the exposure-response function of a single exposure (`z1`) with the second exposure (`z2`) fixed at different quantiles. It uses `geom_smooth` to show the relationship and `facet_grid` to separate by quantile.

```r
ggplot(pred.resp.bivar.levels, aes(z1, est)) +
    geom_smooth(aes(col = quantile), stat = "identity") +
    facet_grid(variable2 ~ variable1) +
    ggtitle("h(expos1 | quantiles of expos2)") +
    xlab("expos1")
```

--------------------------------

### Extract Posterior Inclusion Probabilities (PIP) from BKMR Fit

Source: https://jenfb.github.io/bkmr/overview

Estimates the posterior inclusion probability (PIP) for each exposure variable in a fitted BKMR model. This function is used when variable selection is enabled during model fitting (varsel = TRUE). It returns a data frame with variables and their corresponding PIPs.

```r
ExtractPIPs(fitkm)
```

--------------------------------

### Calculate Overall Risk Summaries with BKMR

Source: https://jenfb.github.io/bkmr/overview

Computes summary statistics of the predictor-response function to assess the overall effect of predictors. It compares the function's value when all predictors are at specified quantiles versus when they are at a fixed quantile (defaulting to the 50th percentile). Requires fitted BKMR model (`fit`), outcome (`y`), and predictors (`Z`, `X`).

```r
risks.overall <- OverallRiskSummaries(fit = fitkm, y = y, Z = Z, X = X, 
                                      qs = seq(0.25, 0.75, by = 0.05), 
                                      q.fixed = 0.5, method = "exact")
risks.overall
```

--------------------------------

### R: Plot Estimated Health Risks Comparison

Source: https://jenfb.github.io/bkmr/SimData1

This R code generates plots to compare estimated health risks (h) from a full BKMR model against those estimated using Gaussian predictive processes with 100 and 50 knots. It visualizes the agreement between the different model approximations by plotting estimated values against each other with a 1:1 reference line.

```r
hi_est <- ComputePostmeanHnew(fitkm)$postmean
hi_knot100 <- ComputePostmeanHnew(fitkm_knots100)$postmean
hi_knot50 <- ComputePostmeanHnew(fitkm_knots50)$postmean

par(mfrow = c(1, 2))
plot(hi_est, hi_knot100, xlab = "Full model", 
     ylab = "100 knots", main = "Estimated h")
abline(0, 1, col = "red")
plot(hi_est, hi_knot50, xlab = "Full model", 
     ylab = "50 knots", main = "Estimated h")
abline(0, 1, col = "red")
```

--------------------------------

### PredictorResponseBivarLevels Function for BKMR

Source: https://jenfb.github.io/bkmr/overview

This R function, `PredictorResponseBivarLevels`, refines the bivariate predictor-response analysis by allowing specific quantiles (`qs`) of one predictor to be fixed while examining the relationship with another. It takes the output from `PredictorResponseBivar` as input.

```r
pred.resp.bivar.levels <- PredictorResponseBivarLevels(
  pred.resp.df = pred.resp.bivar, 
                                                       
  Z = Z, qs = c(0.1, 0.5, 0.9))
```

--------------------------------

### Investigate Prior Distributions for RM/RM Parameters using R

Source: https://jenfb.github.io/bkmr/overview

The `InvestigatePrior` function fits univariate Kernel Machine Regression (KMR) models with a Gaussian kernel to assess how different rmrm parameter values influence the exposure-response function. It takes outcome data (y), exposure matrix (Z), and covariate matrix (X) as input. The function allows specification of rmrm values via `r.seq` or derived values through `q.seq`, outputting model fits for analysis.

```r
priorfits <- InvestigatePrior(y = y, Z = Z, X = X, 
                              q.seq = c(2, 1/2, 1/4, 1/16))

```

--------------------------------

### Plot Univariate Exposure-Response Functions

Source: https://jenfb.github.io/bkmr/SimData1

Visualizes the univariate exposure-response function by plotting the relationship between a single exposure and the outcome. Other exposures are fixed at a specified quantile (defaulting to the median). Uses the `PredictorResponseUnivar` function and `ggplot2` for plotting.

```r
pred.resp.univar <- PredictorResponseUnivar(fit = fitkm)
```

```r
ggplot(pred.resp.univar, aes(z, est, ymin = est - 1.96*se, 
                             ymax = est + 1.96*se)) + 
    geom_smooth(stat = "identity") + 
    facet_wrap(~variable, ncol = 4) + 
    xlab("expos") + 
    ylab("h(expos)")
```

--------------------------------

### R: Visualize Univariate Exposure-Response Relationships

Source: https://jenfb.github.io/bkmr/SimData1

This R code generates plots to visualize the estimated univariate exposure-response functions from a BKMR model against the true functions. It fixes all but one exposure at their median values to isolate the effect of each individual exposure. The output is a faceted plot comparing estimated (blue) and true (red) relationships.

```r
meds <- apply(expos, 2, median)
df_preds_true <- data.frame()
for (i in 1:ncol(expos)) {
    new_grid <- as.list(meds)
    new_grid[[i]] <- seq(min(expos[, i]), max(expos[, i]), length.out = 50)
    new_grid <- expand.grid(new_grid)
    pred <- with(new_grid, htrue_fun(z1, z2, z4, z5, z7))
    pred_center <- pred - mean(pred)
    df0 <- data.frame(variable = paste0("z", i), 
                      expos = new_grid[, i],
                      z = (new_grid[, i] - mean(expos[, i]))/sd(expos[, i]),
                      est = pred_center)
    df_preds_true <- rbind(df_preds_true, df0)
}
df_preds_true$se <- NA

ggplot(pred.resp.univar, aes(z, est, ymin = est - 1.96*se, 
                             ymax = est + 1.96*se)) + 
    geom_smooth(stat = "identity") + 
    geom_line(data = df_preds_true, col = "red", lty = 2) +
    facet_wrap(~variable, ncol = 4) +
    xlab("expos") +
    ylab("h(expos)")
```

--------------------------------

### Visualize Single Exposure Interaction Summaries (R)

Source: https://jenfb.github.io/bkmr/SimData1

Visualizes the single-exposure interaction summaries computed by `SingVarIntSummaries` using `ggplot2`. The plot displays the estimated interaction effects and their confidence intervals, with a horizontal line at zero to indicate no interaction. This visualization aids in understanding how changes in one exposure's quantile affect health risks when other exposures are held at different quantile levels.

```r
ggplot(risks.int, aes(variable, est, ymin = est - 1.96*sd, ymax = est + 1.96*sd)) + 
    geom_pointrange(position = position_dodge(width = 0.75)) +
    geom_hline(yintercept = 0, lty = 2, col = "brown") +
    coord_flip()
```

--------------------------------

### Plot Three-way Interaction Surfaces in R

Source: https://jenfb.github.io/bkmr/SimData1

Visualizes the results of the three-way interaction analysis using ggplot2. It creates facet-wrapped raster plots showing the bivariate exposure-response surface for `z1` and `z2` at different fixed levels of `z5`.

```r
ggplot(pred.bivar.by.expos5, aes(z1, z2, fill = est)) +
    geom_raster() +
    facet_wrap(~expos5) +
    scale_fill_gradientn(colours = c("#0000FFFF","#FFFFFFFF","#FF0000FF")) +
    xlab(unique(pred.bivar.by.expos5$variable1)) +
    ylab(unique(pred.bivar.by.expos5$variable2))
```

--------------------------------

### Calculate Single Exposure Risk Summaries (R)

Source: https://jenfb.github.io/bkmr/SimData1

Computes summaries of health risks associated with changes in a single exposure variable, while holding other exposures fixed at specified quantiles. It requires a fitted BKMR model (`fitkm`), a vector of quantiles to compare (`qs.diff`), and a vector of quantiles for fixed exposures (`q.fixed`). The output is a tibble containing estimates and standard deviations for different exposure scenarios.

```r
risks.singvar <- SingVarRiskSummaries(
    fit = fitkm, qs.diff = c(0.25, 0.75), 
    q.fixed = c(0.25, 0.50, 0.75))

subset(risks.singvar, variable %in% c("z1", "z2", "z4", "z5", "z7"))
```

--------------------------------

### Scale Exposure Variables for BKMR

Source: https://jenfb.github.io/bkmr/SimData1

Applies z-scoring to the exposure variables matrix (`expos`). Scaling ensures that all exposure variables have a mean of 0 and a standard deviation of 1, which is often recommended for kernel-based methods like BKMR.

```R
scale_expos <- scale(expos)
```

--------------------------------

### Hierarchical vs. Component-wise Variable Selection in BKMR

Source: https://jenfb.github.io/bkmr/overview

This snippet demonstrates fitting a BKMR model using both component-wise and hierarchical variable selection. Hierarchical selection groups correlated variables (z1z1 and z3z3 into Group 1), allowing for group-level and within-group variable selection. Component-wise selection evaluates each variable independently. The 'groups' argument in `kmbayes` specifies the group assignment for each variable in the Z matrix.

```R
set.seed(111)
fitkm_corr <- kmbayes(y = d2$y, Z = d2$Z, X = d2$X, iter = 10000, varsel = TRUE, verbose = FALSE)
fitkm_hier <- kmbayes(y = d2$y, Z = d2$Z, X = d2$X, iter = 10000, varsel = TRUE, 
                      groups = c(1,2,1,3), verbose = FALSE)
```

--------------------------------

### PredictorResponseBivar Function for BKMR

Source: https://jenfb.github.io/bkmr/overview

The `PredictorResponseBivar` R function calculates the bivariate predictor-response function from a fitted BKMR model. It allows for the visualization of the joint effect of two predictors on the outcome. The `min.plot.dist` argument controls the minimum distance for plotting.

```r
pred.resp.bivar <- PredictorResponseBivar(fit = fitkm, min.plot.dist = 1)
```

--------------------------------

### Predict Bivariate Exposure-Response at Quantiles in R

Source: https://jenfb.github.io/bkmr/SimData1

Estimates the bivariate exposure-response function while fixing the second exposure at specified quantiles. This is useful for investigating the impact of one exposure while varying another. The `qs` argument defines the quantiles for the second exposure.

```r
pred.resp.bivar.levels <- PredictorResponseBivarLevels(
    pred.resp.bivar, scale_expos, qs = c(0.10, 0.5, 0.90))
```

--------------------------------

### Compute Single-Predictor Health Risks

Source: https://jenfb.github.io/bkmr/overview

Calculates the contribution of individual predictors to the response by comparing risk at different quantiles (e.g., 25th vs. 75th percentile) while holding other predictors constant. Uses the 'SingVarRiskSummaries' function, which requires a fitted model ('fitkm'), response variable ('y'), predictor matrices ('Z', 'X'), and specified quantiles ('qs.diff', 'q.fixed').

```r
risks.singvar <- SingVarRiskSummaries(fit = fitkm, y = y, Z = Z, X = X, 
                                      qs.diff = c(0.25, 0.75), 
                                      q.fixed = c(0.25, 0.50, 0.75),
                                      method = "exact")
risks.singvar
```

--------------------------------

### Calculate Single Exposure Interaction Summaries (R)

Source: https://jenfb.github.io/bkmr/SimData1

Calculates summaries of interaction effects between single exposures by comparing health risks when other exposures are fixed at different quantiles. This function, `SingVarIntSummaries`, takes a fitted BKMR model (`fitkm`), quantiles for comparison (`qs.diff`), and quantiles for fixed exposures (`qs.fixed`). The output is a tibble detailing the interaction estimates and their standard deviations for each variable.

```r
risks.int <- SingVarIntSummaries(fit = fitkm, 
qs.diff = c(0.25, 0.75), qs.fixed = c(0.25, 0.75))
risks.int
```

--------------------------------

### PredictorResponseUnivar Function for BKMR

Source: https://jenfb.github.io/bkmr/overview

This R function, `PredictorResponseUnivar`, is used within the BKMR framework to calculate the univariate predictor-response function. It takes a fitted BKMR model object as input and returns a data frame suitable for plotting the relationship between individual predictors and the outcome.

```r
pred.resp.univar <- PredictorResponseUnivar(fit = fitkm)
```

--------------------------------

### Plot Prior Fits to Visualize RM/RM Parameter Effects using R

Source: https://jenfb.github.io/bkmr/overview

The `PlotPriorFits` function visualizes the results generated by `InvestigatePrior`, illustrating the impact of rmrm parameters on the exposure-response relationship. It requires the outcome data (y), exposure matrix (Z), covariate matrix (X), and the fits object from `InvestigatePrior`. The plots show correlations based on exposure 'distance' and estimated exposure responses for different rmrm values across exposure variables.

```r
PlotPriorFits(y = y, Z = Z, X = X, 
              fits = priorfits)

```

--------------------------------

### Calculate Overall Cumulative Risk Summaries in R

Source: https://jenfb.github.io/bkmr/SimData1

Computes overall cumulative risk summaries by comparing the exposure-response function across a range of quantiles for all exposures against a fixed quantile (default is median). The `OverallRiskSummaries` function calculates the estimated effect (`est`) and its standard deviation (`sd`).

```r
risks.overall <- OverallRiskSummaries(fit = fitkm, qs = seq(0.25, 0.75, by = 0.05), q.fixed = 0.5)
```

--------------------------------

### Plot Univariate Predictor-Response Function with ggplot2

Source: https://jenfb.github.io/bkmr/overview

Generates a plot of the univariate predictor-response function using ggplot2. It visualizes the relationship between a single exposure (z) and the outcome (est) with confidence intervals, faceted by exposure variable. Assumes `pred.resp.univar` object is available from `PredictorResponseUnivar`.

```r
library(ggplot2)
ggplot(pred.resp.univar, aes(z, est, ymin = est - 1.96*se, ymax = est + 1.96*se)) + 
    geom_smooth(stat = "identity") + 
    facet_wrap(~ variable) + 
  ylab("h(z)")
```

--------------------------------

### Plot Bivariate Exposure-Response Surface in R

Source: https://jenfb.github.io/bkmr/SimData1

Generates a raster plot of the bivariate exposure-response surface using ggplot2. It facets the plot by variable pairs and uses a gradient fill to represent the estimated effect. Points outside the specified `min.plot.dist` are grayed out.

```r
ggplot(pred.resp.bivar, aes(z1, z2, fill = est)) +
    geom_raster() +
    facet_grid(variable2 ~ variable1) +
    scale_fill_gradientn(colours = c("#0000FFFF","#FFFFFFFF","#FF0000FF")) +
    xlab("expos1") +
    ylab("expos2") +
    ggtitle("h(expos1, expos2)")
```

--------------------------------

### Estimate Three-way Interaction Bivariate Surfaces in R

Source: https://jenfb.github.io/bkmr/SimData1

Calculates bivariate exposure-response functions for pairs of exposures (z7, z1) while fixing a third exposure (z5) at different probability levels (quantiles). This helps visualize three-way interactions. Multiple calls to `PredictorResponseBivar` are used for different probabilities.

```r
pred.bivar.10 <- PredictorResponseBivar(
    fit = fitkm, min.plot.dist = 1, z.pairs = rbind(c(7,1)), 
    whichz3 = 5, prob = 0.10)
pred.bivar.50 <- PredictorResponseBivar(
    fit = fitkm, min.plot.dist = 1, z.pairs = rbind(c(7,1)), 
    whichz3 = 5, prob = 0.50)
pred.bivar.90 <- PredictorResponseBivar(
    fit = fitkm, min.plot.dist = 1, z.pairs = rbind(c(7,1)), 
    whichz3 = 5, prob = 0.90)

pred.bivar.10$expos5 <- "10% of z5"
pred.bivar.50$expos5 <- "50% of z5"
pred.bivar.90$expos5 <- "90% of z5"

pred.bivar.by.expos5 <- rbind(pred.bivar.10, pred.bivar.50, pred.bivar.90)
pred.bivar.by.expos5$expos5 <- as.factor(pred.bivar.by.expos5$expos5)
```

--------------------------------

### Plot Bivariate Predictor-Response Function with ggplot2

Source: https://jenfb.github.io/bkmr/overview

Creates a plot of the bivariate predictor-response function using ggplot2, visualizing the relationship between two exposures (z1, z2) and the outcome (est) as a heatmap. It uses `geom_raster` and allows faceting by exposure pairs. Assumes `pred.resp.bivar` object is available from `PredictorResponseBivar`.

```r
ggplot(pred.resp.bivar, aes(z1, z2, fill = est)) + 
    geom_raster() + 
    facet_grid(variable2 ~ variable1) + 
  scale_fill_gradientn(colours=c("#0000FFFF","#FFFFFFFF","#FF0000FF")) + 
  xlab("expos1") + 
  ylab("expos2") + 
  ggtitle("h(expos1, expos2)")
```

--------------------------------

### R: Calculate Root Mean Square Error (rMSE) for BKMR Approximations

Source: https://jenfb.github.io/bkmr/SimData1

This R code calculates the root mean square error (rMSE) to quantify the difference between the health risks estimated by a full BKMR model and those estimated by approximations using 100 and 50 knots. The rMSE values are presented in a data frame, indicating the accuracy of the approximations.

```r
rMSE <- data.frame(
    knots50 = sqrt(mean((hi_knot50 - hi_est)^2)), 
    knots100 = sqrt(mean((hi_knot100 - hi_est)^2))
)
rMSE
```

--------------------------------

### R: Evaluate BKMR Overall Fit

Source: https://jenfb.github.io/bkmr/SimData1

This R code defines a function to calculate the true exposure-response function and then plots the estimated health risks against the true health risks, with a reference line at y=x. It helps evaluate how well the BKMR model captures the overall exposure-response relationship.

```r
htrue_fun <- function(z1, z2, z4, z5, z7) {
    alpha0 <- 2
    alpha1 <- 20
    TKt <- 0.5
    K1 <- 1.5
    K2 <- 3
    K4 <- 4.5
    K5 <- 1
    K7 <- 0.5
    R00 <- 1
    lambda <- 2
    R0 <- R00 + lambda*z7/(K7 + z7)
    num <- TKt + z1/K1 + z2/K2
    alpha0 + alpha1*num/(1 + num + z4/K4 + z5/K5)*R0
}

hi_true <- htrue_fun(z1 = expos[, "z1"], z2 = expos[, "z2"], z4 = expos[, "z4"], z5 = expos[, "z5"], z7 = expos[, "z7"])
hi_est <- ComputePostmeanHnew(fitkm)$postmean

plot(hi_true, hi_est, xlab = "True h", ylab = "Estimated h")
abline(0, 1, col = "red")
```

--------------------------------

### Plot Overall Cumulative Risk Summaries in R

Source: https://jenfb.github.io/bkmr/SimData1

Generates a plot of the overall cumulative risk summaries using ggplot2. It displays the estimated risk (`est`) against the quantile, with error bars representing the 95% confidence interval (est +/- 1.96*sd). A horizontal line at y=0 is included for reference.

```r
ggplot(risks.overall, aes(quantile, est, ymin = est - 1.96*sd, 
                          ymax = est + 1.96*sd)) +
geom_hline(yintercept = 0, lty = 2, col = "brown") +
geom_pointrange()
```

--------------------------------

### Plot Bivariate Predictor-Response Function Levels with ggplot2

Source: https://jenfb.github.io/bkmr/overview

Generates plots showing the predictor-response function for one exposure while varying quantiles of a second exposure. This is useful for understanding conditional relationships when direct bivariate plots are complex. It uses `geom_smooth` with `stat = "identity"` and faceting. Assumes `pred.resp.bivar.levels` object is available from `PredictorResponseBivarLevels`.

```r
ggplot(pred.resp.bivar.levels, aes(z1, est)) + 
    geom_smooth(aes(col = quantile), stat = "identity") + 
    facet_grid(variable2 ~ variable1) + 
    ggtitle("h(expos1 | quantiles of expos2)") + 
  xlab("expos1")
```

--------------------------------

### Rename Exposure and Covariate Variables

Source: https://jenfb.github.io/bkmr/SimData1

Renames the columns of the exposure matrix to 'z1' through 'z7' and the covariate matrix to 'x'. This is done for consistency with common notation in BKMR literature.

```R
colnames(covar) <- "x"
colnames(expos) <- paste0("z", 1:ncol(expos))
```

--------------------------------

### Predict Bivariate Exposure-Response Function in R

Source: https://jenfb.github.io/bkmr/SimData1

Estimates the bivariate exposure-response function, useful for creating contour or image plots. It allows specifying a minimum plotting distance to focus on areas near observed data points. The output is a data frame suitable for plotting.

```r
pred.resp.bivar <- PredictorResponseBivar(fit = fitkm, 
                                          min.plot.dist = 0.5,
                                          z.pairs = expos.pairs)
```