#' Periodic kernel
#'
#' @param X covariate (dimension Q x N; i.e., covariates x samples)
#' @param sigma scalar parameter
#' @param rho scalar bandwidth parameter
#' @param period period length (in days)
#' @param jitter small scalar to add to off-diagonal of gram matrix
#' (for numerical underflow issues)
#' @return Gram Matrix (N x N) (e.g., the Kernel evaluated at
#' each pair of points)
#' @export
<- (X, =1, rho=1, period=24, =0){
<- (((X)))
G <- ^2 * (-2*((*/period)^2)/(rho^2)) + *(())
}
#' Trivially expands the dimensions of the posterior estimates objects in a MAP estimate basset result
#'
#' @param fit a bassetfit object
#' @return NULL
#' @export
#' @examples
#' fit <- fix_MAP_dims(fit)
<- (fit) {
(((fit$fit$Eta)) == 2) {
(fit$fit$Eta) <- ((fit$fit$Eta),1)
(fit$fit$Lambda) <- ((fit$fit$Lambda),1)
(fit$fit$Sigma) <- ((fit$fit$Sigma),1)
}
(fit)
}
#' Pulls a list of fitted basset models from the designated output directory
#'
#' @param tax_level taxonomic level at which to agglomerate data
#' @param DLM if TRUE, looks for DLM model fits instead of GP model fits
#' @param MAP use MAP estimate model output instead of full posterior output
#' @return NULL
#' @export
#' @examples
#' model_list <- get_fitted_model_list(tax_level = "ASV", DLM = TRUE, MAP = FALSE)
<- (tax_level = "ASV", DLM = , MAP = ) {
(DLM) {
pattern_str <- "*_labraduckfit.rds"
regexpr_str <- "_labraduckfit.rds"
} {
pattern_str <- "*_bassetfit.rds"
regexpr_str <- "_bassetfit.rds"
}
(MAP) {
level_dir <- ("output","model_fits",(tax_level,"_MAP"))
} {
level_dir <- ("output","model_fits",tax_level)
}
fitted_models <- (path=level_dir, =pattern_str, full.names=, recursive=)
# pull out some useful summary information: a list of fitted hosts, models, and dimensions
hosts <- ((fitted_models, (x) { idx <- (regexpr_str, x); ((x, idx-3, idx-1)) } ))
(MAP) {
path <- ("output","model_fits",(tax_level,"_MAP/",hosts[1],regexpr_str))
} {
path <- ("output","model_fits",(tax_level,"/",hosts[1],regexpr_str))
}
fit <- (path)
((hosts=hosts,
pattern_str=pattern_str,
regexpr_str=regexpr_str,
model_list=fitted_models,
=fit$fit$,
n_samples=fit$fit$iter))
}
#' Set up a basic ALR prior
#'
#' @param D number of features including reference (where the ALR will represent D-1)
#' @param total_variance scale of the log variance
#' @return list containing inverse Wishart parameters degrees of freedom and scale matrix
#' @export
#' @examples
#' params <- get_Xi(D=100, log_var_scale=1)
#' Sigma <- matrixsampling::rinvwishart(1, params$upsilon, params$Xi)[,,1]
<- (, total_variance=1) {
upsilon <- -1+10 # specify low certainty/concentration
GG <- ((-1), -1) # log contrast for ALR with last taxon as reference
Xi <- GG%*%(()*total_variance)%*%(GG) # take diag as covariance over log abundances
Xi <- Xi*(upsilon--1)
((upsilon=upsilon, Xi=Xi))
}
#' Generate empirical estimates for variance components to use in the Gaussian process
#'
#' @param data a phyloseq object
#' @return list containing median log abundance taxon and mean empirical total variance
#' across hosts using this ALR representation
#' @import phyloseq
#' @import driver
#' @export
<- () {
metadata <- sample_data()
sample_status <- metadata$sample_status # replicate labels
# pick an intermediate abundance ALR reference taxon
counts <- otu_table()@.Data
alr_ref <- (((counts, 2, )) == ((counts)/2))
# apply the ALR we'll model in
log_ratios <- alr(counts + 0.5, d=alr_ref) # samples x taxa
# get mean within-host ALR total variance
host_var <- ()
(host (metadata$sname)) {
retain_sids <- metadata[metadata$sname == host,]$sample_id
((retain_sids) > 3) {
lr <- log_ratios((log_ratios) %in% retain_sids),]
lr <- (lr, =)
sample_var <- ((lr))
host_var <- (host_var, (sample_var))
}
}
mean_total_var <- (host_var)
# return squared exponential kernel variance as 45% (90%/2) the total empirical ALR variance and the periodic
# kernel variance as 5% (10%/2) the total empirical ALR variance
((total_variance=mean_total_var, alr_ref=alr_ref))
}
#' Build the design matrix
#'
#' @param metadata metadata data.frame associated with the data phyloseq object
#' @param metadata.diet metadata data.frame parsed from separate diet and climate data
#' @return matrix with observation time (days), diet PCs, and climate (rainfall and max temp) covariates
#' @export
<- (metadata, metadata.diet = ) {
baseline_date <- metadata$collection_date[1]
X <- (metadata$collection_date, (x) ((x, baseline_date, ="days"))) + 1
# render X a row matrix
(X) <- (1, (X))
(!(metadata.diet)) {
# just a few diet PCs for now
X <- (X, (metadata.diet,("diet_PC1","diet_PC2","diet_PC3"),=F]))
X <- (X, (metadata.diet,("rain_monthly","tempmax_monthly")]))
}
# do samples line up?
(((metadata$sample_id == metadata.diet$sample_id) == ) > 0) {
("Metadata and metadata-diet samples don't line up!\n")
}
# standardize diet and climate predictors and impute NAs
((X) > 1) {
(i 2:(X)) {
X[i,] <- (X[i,])
((((X[i,]))) > 0) {
x <- 1:(X[i,])
y <- X[i,]
X[i,] <- (x = x, y = y, xout = x, method = "linear")$y
}
}
}
(X)
}
#' Build the design matrix for prediction by interpolating within the
#'
#' @param X existing design matrix
#' @return matrix with observation time (days), diet PCs, and climate (rainfall and max temp) covariates interpolated from
#' earliest to latest observation date
#' @export
<- (X) {
last_time_point <- (X[1,])
X_predict <- (1:last_time_point)
n_features <- (X)-1 # these are the features to (linearly) interpolate
(n_features > 0) {
(f 1:n_features) {
X_predict <- (X_predict, ((x = X[1,], y = X[f+1,], xout = 1:last_time_point, method = "linear")$y))
}
}
(X_predict)
}
#' Define bandwidth of squared exponential kernel
#'
#' @param min_correlation minimum correlation to assume between (within-host) samples
#' @param days_to_baseline days at which squared exponential kernel decays to baseline correlation of ~0.1
#' @return bandwidth parameter for squared exponential kernel
#' @export
<- (min_correlation = 0.1, days_to_baseline = 90) {
# back-calculate the squared exponential bandwidth parameter by finding a bandwidth that gives
# a desired minimum correlation at the number of days specified by days_to_baseline
rho <- (-days_to_baseline^2/(2*(min_correlation))) # back calculate the decay
(rho)
}
#' Define a kernel (function) over samples
#'
#' @param kernel_scale total variance for the composite kernel
#' @param proportions proportion variance to attribute to each of 3 kernels (see details)
#' @param rho bandwidth for SE kernel
#' @details Composite kernel is built from (1) squared exponential kernel (base autocorrelation component),
#' (2) seasonal kernel (periodic), and (3) diet and climate component (another squared exponential)
#' @return list containing kernel function and bandwidth parameter
#' @import fido
#' @export
<- (kernel_scale, , min_correlation = 0.1, days_to_baseline = 90) {
rho <- (min_correlation = min_correlation, days_to_baseline = days_to_baseline)
# back-calculate the squared exponential bandwidth parameter by finding a bandwidth that gives
# a desired minimum correlation at the number of days specified by SE_days_to_baseline
<- (X) {
<- 1e-08
<- ()/(()) # just in case we pass something that isn't a composition
part.1 <- kernel_scale * [1]
part.2 <- kernel_scale * [2]
part.3 <- kernel_scale * [3]
SE(X[1,,=F], = (part.1), rho = rho, = ) +
(X[1,,=F], = (part.2), rho = 1, period = 365, = ) +
SE(X[2:6,,=F], = (part.3), rho = 1, = ) # a few diet PCs
}
((rho = rho, = ))
}
#' Utility function to calculate RMSE
#'
#' @param y1 count vector
#' @param y2 count vector
#' @param pc optional pseudocount parameter
#' @return root mean squared error of log counts
#' @export
<- (y1, y2, pc = 0.5) {
(((1:(y1), (yy) {
((y1[yy] + pc) - (y2[yy] + pc))^2
}))/(y1))
}
#' Utility function to calculate fold change (as error)
#'
#' @param y1 count vector
#' @param y2 count vector
#' @param pc optional pseudocount parameter
#' @return average fold difference between vectors
#' @export
<- (y1, y2, pc = 0.5) {
((1:(y1), (yy) {
a <- (y1[yy] + pc, y2[yy] + pc)
b <- (y1[yy] + pc, y2[yy] + pc)
a/b
})) - 1
}
#' Fit a Gaussian process to a single host series using basset
#'
#' @param data a phyloseq object
#' @param host host short name (e.g. ACA)
#' @param taxa_covariance list of prior covariance parameters over taxa
#' @param sample_covariance kernel function
#' @param rho bandwidth for SE kernel
#' @param tax_level taxonomic level at which to agglomerate data
#' @param alr_ref index of reference ALR coordinate
#' @param n_samples number of posterior samples to draw
#' @param MAP compute MAP estimate only (as single posterior sample)
#' @param holdout_proportion if non-zero, proportion of host's sample to use as a test set
#' @param return_model if TRUE, returns the fitted model instead of saving it
#' @param scramble if TRUE, host samples are scrambled in time; this serve the diagnostic purpose of
#' identifying whether some of the correlation we see between individual hosts' dynamics
#' are the result of systematic differences in abundances (etc.)
#' @details Fitted model and metadata saved to designated model output directory.
#' @return NULL
#' @import phyloseq
#' @import fido
#' @export
#' @examples
#' tax_level <- "ASV"
#' data <- load_data(tax_level = tax_level)
#' params <- formalize_parameters(data)
#' sample_covariance <- get_Gamma(kernel_scale = 2, proportions = c(1, 0, 0), min_correlation = 0.1, days_to_baseline = 90)
#' taxa_covariance <- get_Xi(phyloseq::ntaxa(data), total_variance = 1)
#' fit_GP(data, host = "GAB", taxa_covariance = taxa_covariance, sample_covariance = sample_covariance, tax_level = tax_level, alr_ref = params$alr_ref, MAP = TRUE)
<- (, host, taxa_covariance, sample_covariance, tax_level = "ASV", alr_ref = , n_samples = 100, MAP = , holdout_proportion = 0, return_model = , scramble = ) {
(MAP) {
(("Fitting fido::basset model (MAP) to host ",host,"\n"))
(holdout_proportion > 0) {
("Leave-one-out predictions only work for full posterior estimates (currently).\n")
}
} {
(("Fitting fido::basset model to host ",host,"\n"))
}
# global assign is a hack seemingly necessary for this phyloseq::subset_samples function call
host <<- host
host_data <- subset_samples(, sname == host)
# encode observations as differences from baseline in units of days
host_metadata <- sample_data(host_data)
# read diet and climate covariate data
data.diet <- ("input/ps_w_covs.RDS")
data.name_mapping <- ("input/host_subject_id_to_sname_key.csv")
data.name_mapping <- (data.name_mapping,("sname","host_subject_id2")])
host.num <<- (data.name_mapping[data.name_mapping$sname == host,]$host_subject_id2)
data.diet <- subset_samples(data.diet, host == host.num)
metadata.diet <- sample_data(data.diet)
# pull out the count table
Y <- otu_table(host_data)@.Data
X <- (host_metadata, metadata.diet) # covariates x samples
Y <- (Y) # taxa x samples
# get dimensions
<- (Y)
N <- (Y)
# scramble, if required
(scramble) {
("Scrambling host observations...\n")
(i 1:) {
per_taxon_col_assignments <- (1:N)
Y[i,] <- Y[i,per_taxon_col_assignments]
}
}
# fido uses the D^th element as the ALR reference by default
# if we'd like to use a different reference, do some row shuffling in Y to put the reference at the end
(!(alr_ref)) {
Y <- Y[c((1:,alr_ref),alr_ref),]
}
# strip off sequence variant labels
(Y) <-
(Y) <-
(holdout_proportion > 0) {
holdout_n <- (N*holdout_proportion)
holdout_idx <- (1:N)[1:holdout_n]
X_test <- X,holdout_idx]
X <- X,(1:N, holdout_idx)]
Y_test <- Y,holdout_idx]
Y <- Y,(1:N, holdout_idx)]
}
# define the prior over baselines
alr_ys <- driver::alr(((Y) + 0.5))
alr_means <- (alr_ys)
Theta <- (X) (alr_means, -1, (X))
(MAP) {
n_samples <- 0
}
fit <- fido::basset(Y, X, taxa_covariance$upsilon, Theta, sample_covariance$, taxa_covariance$Xi,
n_samples = n_samples, ret_mean = MAP,
b2 = 0.98, step_size = 0.004, eps_f = 1e-11, eps_g = 1e-05,
max_iter = 10000L, optim_method = "adam")
(MAP) {
# fill out dimensions; some function (including fido predictive functions) expect arrays
# not matrices here
fit <- (fit)
}
(holdout_proportion > 0) {
predicted <- (fit, X_test, response = "Y", iter = fit$iter) # taxa x holdout samples x posterior samples
err1 <- ((1:(predicted)[2], (i) {
# for each held out sample
ref <- Y_test,i]
(1:(predicted)[3], (j) {
# for each posterior sample
(ref, predicted,i,j])
})
}))
err2 <- ((1:(predicted)[2], (i) {
# for each held out sample
ref <- Y_test,i]
(1:(predicted)[3], (j) {
# for each posterior sample
(ref, predicted,i,j])
})
}))
(( = err1, = err2))
}
# pack up results
fit_obj <- (Y = Y, alr_ys = alr_ys, alr_ref = alr_ref, X = X, fit = fit)
# a hack: for later prediction, these will need to be available in the workspace for Gamma
fit_obj$kernelparams$rho <- sample_covariance$rho
(return_model) {
(fit_obj)
}
# save results
(MAP) {
save_dir <- (("output","model_fits",(tax_level,"_MAP")))
} {
save_dir <- (("output","model_fits",tax_level))
}
(fit_obj, (save_dir,(host,"_bassetfit.rds")))
}
#' Fit a dynamic linear model to a single host series using labraduck
#'
#' @param data a phyloseq object
#' @param host host short name (e.g. ACA)
#' @param taxa_covariance list of prior covariance parameters over taxa
#' @param var_scale combined scale of the total variance components Sigma and Gamma (default 1)
#' @param tax_level taxonomic level at which to agglomerate data
#' @param alr_ref index of reference ALR coordinate
#' @param n_samples number of posterior samples to draw
#' @param MAP compute MAP estimate only (as single posterior sample)
#' @param use_covariates if TRUE uses available rain, temperature, and diet data as covariates in the model
#' @param split_diet if TRUE, splits diet PC1 into positive and negative subsets to treat as two 'environments'
#' @details Fitted model and metadata saved to designated model output directory.
#' @return NULL
#' @import phyloseq
#' @import fido
#' @export
#' @examples
#' tax_level <- "ASV"
#' data <- load_data(tax_level = tax_level)
#' params <- formalize_parameters(data)
#' taxa_covariance <- get_Xi(phyloseq::ntaxa(data), total_variance = 1)
#' fit_DLM(data, host = "GAB", taxa_covariance = taxa_covariance, tax_level = tax_level, alr_ref = params$alr_ref, MAP = TRUE)
<- (, host, taxa_covariance, var_scale = 1, tax_level = "ASV", alr_ref = ,
n_samples = 100, MAP = , use_covariates = , split_diet = ) {
(MAP) {
(("Fitting fido::labraduck model (MAP) to host ",host,"\n"))
} {
(("Fitting fido::labraduck model to host ",host,"\n"))
}
# global assign is a hack seemingly necessary for this phyloseq::subset_samples function call
host <<- host
host_data <- subset_samples(, sname == host)
# encode observations as differences from baseline in units of days
host_metadata <- sample_data(host_data)
# read diet and climate covariate data
data.diet <- ("input/ps_w_covs.RDS")
data.name_mapping <- ("input/host_subject_id_to_sname_key.csv")
data.name_mapping <- (data.name_mapping,("sname","host_subject_id2")])
host.num <<- (data.name_mapping[data.name_mapping$sname == host,]$host_subject_id2)
data.diet <- subset_samples(data.diet, host == host.num)
metadata.diet <- sample_data(data.diet)
(split_diet) {
use_covariates =
env_idx1 <- (metadata.diet$diet_PC1 < 0)
env_idx2 <- (metadata.diet$diet_PC1 >= 0)
env_idxs <- (env_idx1, env_idx2)
} {
env_idxs <- (1:(metadata.diet$diet_PC1))
}
fits <- ()
(idx_it 1:(env_idxs)) {
idx <- env_idxs[idx_it]]
days <- host_metadata$collection_date[idx]
day0 <- (days)
days <- (((days, (x) (x, day0, = "days")))) + 1
# pull out the count table
Y <- (otu_table(host_data)@.Data) # taxa x samples
Y <- Y,idx]
# strip off sequence variant labels
(Y) <-
(Y) <-
# clr.Y <- clr_array(Y + 0.5, parts = 1)
T <- (days)
# Build the pseudo-covariate matrix
(use_covariates) {
F <- (0, 9, T)
F[1,] <- 1
(i 1:(days)) {
# # F[2,data$days[i]] <- data$season[i] # season actually worsens the fit here
F[2,days[i]] <- (((metadata.diet$rain_monthly + 0.1)))[i]
F[3,days[i]] <- ((metadata.diet$tempmax_monthly))[i]
F[4,days[i]] <- ((metadata.diet$diet_PC1))[i]
F[5,days[i]] <- ((metadata.diet$diet_PC2))[i]
F[6,days[i]] <- ((metadata.diet$diet_PC3))[i]
F[7,days[i]] <- ((metadata.diet$diet_PC4))[i]
F[8,days[i]] <- ((metadata.diet$diet_PC5))[i]
F[9,days[i]] <- ((metadata.diet$diet_PC6))[i]
}
# Some covariates as NA; impute these with the mean (zero for these scaled vars)
(((F)) > 0) {
na_fill <- ((F), arr.ind = T)
(na_idx 1:(na_fill)) {
i <- na_fill[na_idx,1]
j <- na_fill[na_idx,2]
F[i,j] <- 0
}
}
} {
F <- (1, 1, T)
}
Q <- (F)
<- (Y)
W <- (Q)
# scale covariate-inclusive and covariate-exclusive models to have similar total variance
W <- W/(W)
G <- (Q)
# fido uses the D^th element as the ALR reference by default
# if we'd like to use a different reference, do some row shuffling in Y to put the reference at the end
(!(alr_ref)) {
Y <- Y[c((1:,alr_ref),alr_ref),]
}
# define the prior over baselines
C0 <- W
alr_ys <- driver::alr(((Y) + 0.5))
alr_means <- (alr_ys)
M0 <- (0, Q, -1)
M0[1,] <- alr_means
(MAP) {
n_samples <- 0
}
# I'm giving W about 1/2 the scale of gamma
# My thinking here is that in the gLV models we've been simulating from, we've seen that most of
# the dynamism in the model results from /environmental interactions/, not innate volatility.
# In my mind, giving W a smaller scale than gamma is consistent with relatively more volatility
# presumed to result from interactions with the environment.
fit <- labraduck(Y = Y, upsilon = taxa_covariance$upsilon, Xi = taxa_covariance$Xi, F = F, G = G, W = W, M0 = M0, C0 = C0,
observations = days, gamma_scale = (var_scale * 2/3), W_scale = (var_scale * 1/3), apply_smoother = MAP,
n_samples = n_samples, ret_mean = MAP)
fits[idx_it]] <- fit
}
(fit_it 1:(fits)) {
fit <- fits[fit_it]]
# pack up results
fit_obj <- (Y = Y, alr_ys = alr_ys, alr_ref = alr_ref, fit = fit)
# save results
(MAP) {
save_dir <- (("output","model_fits",(tax_level,"_MAP")))
} {
save_dir <- (("output","model_fits",tax_level))
}
(fit_obj, (save_dir,(host,"_labraduckfit_",fit_it,".rds")))
}
}
#' Perform k-fold cross-validation for a given host and kernel configuration choice
#'
#' @param data a phyloseq object
#' @param host host short name (e.g. ACA)
#' @param sample_covariance composite kernel function
#' @param rho scalar bandwidth parameter
#' @param tax_level taxonomic level at which to agglomerate data
#' @param holdout_proportion proportion of host's sample to use as a test set
#' @param iterations number of CV iterations to perform; this should be <= n_samples for this host
#' @details average error (currently RMSE of log counts)
#' @return NULL
#' @import phyloseq
#' @export
#' @examples
#' tax_level <- "ASV"
#' data <- load_data(tax_level = tax_level)
#' sample_covariance <- get_Gamma(kernel_scale = 2, proportions = c(1, 0, 0), min_correlation = 0.1, days_to_baseline = 90)
#' taxa_covariance <- get_Xi(phyloseq::ntaxa(data), total_variance = 1)
#' perform_cv(data, host = "GAB", taxa_covariance = taxa_covariance, sample_covariance = sample_covariance, tax_level = tax_level)
<- (, host, taxa_covariance, sample_covariance, tax_level = "ASV", holdout_proportion = 0.2) {
params <- ()
((, host, taxa_covariance, sample_covariance, rho, tax_level, alr_ref = params$alr_ref, n_samples = 100,
MAP = , holdout_proportion = holdout_proportion))
}
