diff --git a/.Rbuildignore b/.Rbuildignore index e13c405..250f717 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -1,3 +1,6 @@ ^.*\.Rproj$ ^\.Rproj\.user$ ^LICENSE\.md$ + +^_pkgdown\.yml$ +^pkgdown$ \ No newline at end of file diff --git a/.github/workflows/R-CMD-check.yaml b/.github/workflows/R-CMD-check.yaml index 56c5fe7..0c0de53 100644 --- a/.github/workflows/R-CMD-check.yaml +++ b/.github/workflows/R-CMD-check.yaml @@ -23,6 +23,7 @@ jobs: env: R_REMOTES_NO_ERRORS_FROM_WARNINGS: true + _R_CHECK_FORCE_SUGGESTS_: false RSPM: ${{ matrix.config.rspm }} GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} diff --git a/.github/workflows/pkgdown.yaml b/.github/workflows/pkgdown.yaml new file mode 100644 index 0000000..a7276e8 --- /dev/null +++ b/.github/workflows/pkgdown.yaml @@ -0,0 +1,48 @@ +# Workflow derived from https://github.com/r-lib/actions/tree/v2/examples +# Need help debugging build failures? Start at https://github.com/r-lib/actions#where-to-find-help +on: + push: + branches: [main, master] + pull_request: + branches: [main, master] + release: + types: [published] + workflow_dispatch: + +name: pkgdown + +jobs: + pkgdown: + runs-on: ubuntu-latest + # Only restrict concurrency for non-PR jobs + concurrency: + group: pkgdown-${{ github.event_name != 'pull_request' || github.run_id }} + env: + GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} + permissions: + contents: write + steps: + - uses: actions/checkout@v4 + + - uses: r-lib/actions/setup-pandoc@v2 + + - uses: r-lib/actions/setup-r@v2 + with: + use-public-rspm: true + + - uses: r-lib/actions/setup-r-dependencies@v2 + with: + extra-packages: any::pkgdown, local::. + needs: website + + - name: Build site + run: pkgdown::build_site_github_pages(new_process = FALSE, install = FALSE) + shell: Rscript {0} + + - name: Deploy to GitHub pages 🚀 + if: github.event_name != 'pull_request' + uses: JamesIves/github-pages-deploy-action@v4.5.0 + with: + clean: false + branch: gh-pages + folder: docs diff --git a/DESCRIPTION b/DESCRIPTION index c165bf2..8071925 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,31 +1,36 @@ Package: radEmu Title: Using Relative Abundance Data to Estimate of Multiplicative Differences in Mean Absolute Abundance -Version: 1.1.0.0 +Version: 1.2.0.0 Authors@R: as.person(c( "David Clausen [aut, cre]", "Amy Willis [aut]", "Sarah Teichman [aut]" )) Description: A differential abundance method for the analysis of microbiome data. radEmu estimates fold-differences in the abundance of taxa across samples relative to "typical" fold-differences. Notably, it does not require pseudocounts, nor choosing a denominator taxon. +URL: https://github.com/statdivlab/radEmu, https://statdivlab.github.io/radEmu/ License: MIT + file LICENSE Encoding: UTF-8 LazyData: true Roxygen: list(markdown = TRUE) -RoxygenNote: 7.3.1 +RoxygenNote: 7.3.2 Depends: MASS, Matrix, methods, + magrittr, + rlang, R (>= 2.10) Suggests: testthat (>= 3.0.0), numDeriv, phyloseq, + TreeSummarizedExperiment, + SummarizedExperiment, knitr, - magrittr, - dplyr, - ggplot2, + dplyr, + ggplot2, stringr, + parallel, rmarkdown Config/testthat/edition: 3 VignetteBuilder: knitr diff --git a/NAMESPACE b/NAMESPACE index d4519a2..145b1b9 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,12 +1,15 @@ # Generated by roxygen2: do not edit by hand +S3method(plot,emuFit) S3method(print,emuFit) export(emuFit) export(emuFit_micro) export(score_test) import(MASS) import(Matrix) +importFrom(magrittr,"%>%") importFrom(methods,as) +importFrom(rlang,.data) importFrom(stats,cov) importFrom(stats,median) importFrom(stats,model.matrix) @@ -14,7 +17,9 @@ importFrom(stats,optim) importFrom(stats,pchisq) importFrom(stats,qnorm) importFrom(stats,rbinom) +importFrom(stats,reorder) importFrom(stats,rnbinom) importFrom(stats,rnorm) importFrom(stats,rpois) +importFrom(stats,setNames) importFrom(stats,weighted.mean) diff --git a/R/d_l_aug_dB.R b/R/d_l_aug_dB.R new file mode 100644 index 0000000..8b469f9 --- /dev/null +++ b/R/d_l_aug_dB.R @@ -0,0 +1,41 @@ + +d_l_aug_dB <- function(Bj, + Bj_constr_no_k_constr, + z, + Y, + X, + Bk_constr_no_j_j_constr, + k_constr, + j, + j_constr, + constraint_fn, + constraint_grad_fn){ + #update parameter determined by null constraint + Bk_constr_j_constr <- + constraint_fn(c(Bk_constr_no_j_j_constr, + Bj[k_constr])) + #create object to store updated two columns of B (j and j_constr) in + Bj_j_constr <- cbind(Bj,0) + + #update object with elements of B_j_constr not subject to null constraint + Bj_j_constr[-k_constr,2] <- Bj_constr_no_k_constr + + #update object with element of B_j_constr subject to null constraint + Bj_j_constr[k_constr,2] <- Bk_constr_j_constr + + #construct log means for Yj and Yj_constr + log_means <- X%*%Bj_j_constr + matrix(z,ncol = 1)%*%matrix(1,nrow = 1, ncol = 2) + + #compute ll + deriv_j <- t(X)%*%(Y[,j,drop = FALSE] - exp(log_means[,1])) + deriv_j_constr <- t(X)%*%(Y[,j_constr,drop = FALSE] - exp(log_means[,2])) + constr_grad <- constraint_grad_fn(c(Bk_constr_no_j_j_constr, + Bj[k_constr])) + constr_grad <- constr_grad[length(constr_grad)] + deriv_j[k_constr] <- deriv_j[k_constr] + deriv_j_constr[k_constr]*constr_grad + + + #return ll + return(c(deriv_j,deriv_j_constr[-k_constr])) + +} \ No newline at end of file diff --git a/R/emuFit.R b/R/emuFit.R index 51a982e..e45d692 100644 --- a/R/emuFit.R +++ b/R/emuFit.R @@ -4,12 +4,17 @@ #' @param X an n x p matrix or dataframe of covariates (optional) #' @param formula a one-sided formula specifying the form of the mean model to be fit #' @param data an n x p data frame containing variables given in \code{formula} -#' @param cluster a numeric vector giving cluster membership for each row of Y to -#' be used in computing GEE test statistics. Default is NULL, in which case rows of -#' Y are treated as independent. +#' @param assay_name a string containing the desired assay name within a `TreeSummarizedExperiment` object. +#' This is only required if Y is a `TreeSummarizedExperiment` object, otherwise this argument does nothing +#' and can be ignored. +#' @param cluster a vector giving cluster membership for each row of Y to be used in computing +#' GEE test statistics. Default is NULL, in which case rows of Y are treated as independent. #' @param penalize logical: should Firth penalty be used in fitting model? Default is TRUE. #' @param B starting value of coefficient matrix (p x J). If not provided, #' B will be initiated as a zero matrix. +#' @param B_null_list list of starting values of coefficient matrix (p x J) for null estimation. This should either +#' be a list with the same length as \code{test_kj}. If you only want to provide starting values for some tests, +#' include the other elements of the list as \code{NULL}. #' @param fitted_model a fitted model produced by a separate call to emuFit; to #' be provided if score tests are to be run without refitting the full unrestricted model. #' Default is NULL. @@ -31,16 +36,22 @@ #' @param use_both_cov logical: should score tests be run using information and #' empirical score covariance evaluated both under the null and full models? #' Used in simulations +#' @param match_row_names logical: Make sure rows on covariate data and response data correspond to +#' the same sample by comparing row names and subsetting/reordering if necessary. #' @param constraint_fn function g defining a constraint on rows of B; g(B_k) = 0 #' for rows k = 1, ..., p of B. Default function is a smoothed median (minimizer of -#' pseudohuber loss). +#' pseudohuber loss). If a number is provided a single category constraint will be used +#' with the provided category as a reference category. #' @param constraint_grad_fn derivative of constraint_fn with respect to its #' arguments (i.e., elements of a row of B) #' @param constraint_param If pseudohuber centering is used (this is the default), #' parameter controlling relative weighting of elements closer and further from center. #' (Limit as \code{constraint_param} approaches infinity is the mean; as this parameter approaches zero, #' the minimizer of the pseudo-Huber loss approaches the median.) -#' @param verbose provide updates as model is being fitted? Defaults to TRUE. +#' @param verbose provide updates as model is being fitted? Defaults to FALSE. If user sets verbose = TRUE, +#' then key messages about algorithm progress will be displayed. If user sets verbose = "development", +#' then key messages and technical messages about convergence will be displayed. Most users who want status +#' updates should set verbose = TRUE. #' @param tolerance tolerance for stopping criterion in full model fitting; once #' no element of B is updated by more than this value in a single step, we exit #' optimization. Defaults to 1e-3. @@ -71,20 +82,42 @@ #' @param c1 numeric: parameter for Armijo line search. Default is 1e-4. #' @param trackB logical: should values of B be recorded across optimization #' iterations and be returned? Primarily used for debugging. Default is FALSE. +#' @param return_nullB logical: should values of B under null hypothesis be returned. Primarily used for debugging. Default is FALSE. +#' @param return_score_components logical: should components of score statistic be returned? Primarily used for debugging. Default is FALSE. #' @param return_both_score_pvals logical: should score p-values be returned using both #' information matrix computed from full model fit and from null model fits? Default is #' FALSE. This parameter is used for simulations - in any applied analysis, type of #' p-value to be used should be chosen before conducting tests. +#' @param remove_zero_comparison_pvals Should score p-values be replaced with NA for zero-comparison parameters? These parameters occur +#' for categorical covariates with three or more levels, and represent parameters that compare a covariate level to the reference level for +#' a category in which the comparison level and reference level both have 0 counts in all samples. These parameters can have misleadingly +#' small p-values and are not thought to have scientifically interesting signals. We recommend removing them before analyzing data further. +#' If TRUE, all zero-comparison parameter p-values will be set to NA. If FALSE no zero-comparison parameter p-values will be set to NA. +#' If a value between 0 and 1, all zero-comparison p-values below the value will be set to NA. +#' Default is \code{0.01}. +#' @param unobserved_taxon_error logical: should an error be thrown if Y includes taxa that have 0 counts for all samples? Default is TRUE. #' #' @return A list containing elements 'coef', 'B', 'penalized', 'Y_augmented', -#' 'I', and 'Dy'. Parameter estimates by -#' covariate and outcome category (e.g., taxon for microbiome data), as well as -#' optionally confidence intervals and p-values, are contained in 'coef'. 'B' -#' contains parameter estimates in matrix format (rows indexing covariates and -#' columns indexing outcome category / taxon). 'penalized' is equal to TRUE -#' if Firth penalty is used in estimation (default) and FALSE otherwise. 'I' and -#' 'Dy' contain an information matrix and empirical score covariance matrix -#' computed under the full model. +#' 'z_hat', 'I', 'Dy', and 'score_test_hyperparams' if score tests are run. +#' Parameter estimates by covariate and outcome category (e.g., taxon for microbiome data), +#' as well as optionally confidence intervals and p-values, are contained in 'coef'. +#' Any robust score statistics and score test p-values are also included in 'coef'. +#' If there are any zero-comparison parameters in the model, a column 'zero_comparison' +#' is also included, which is TRUE for any parameters that compare the level of a categorical +#' covariate to a reference level for a category with only zero counts for both the comparison +#' level and the reference level. This check is currently implemented for an arbitrary design matrix +#' generated using the \code{formula} and \code{data} arguments, and for a design matrix with no more +#' than one categorical covariate if the design matrix \code{X} is input directly. +#' 'B' contains parameter estimates in matrix format (rows indexing covariates and +#' columns indexing outcome category / taxon). +#' 'penalized' is equal to TRUE f Firth penalty is used in estimation (default) and FALSE otherwise. +#' 'z_hat' returns the nuisance parameters calculated in Equation 7 of the radEmu manuscript, +#' corresponding to either 'Y_augmented' or 'Y' if the 'penalized' is equal to TRUE +#' or FALSE, respectively. +#' I' and 'Dy' contain an information matrix and empirical score covariance matrix computed under the +#' full model. +#' 'score_test_hyperparams' contains parameters and hyperparameters related to estimation under the null, +#' including whether or not the algorithm converged, which can be helpful for debugging. #' #' @importFrom stats cov median model.matrix optim pchisq qnorm weighted.mean #' @import Matrix @@ -96,9 +129,11 @@ emuFit <- function(Y, X = NULL, formula = NULL, data = NULL, + assay_name = NULL, cluster = NULL, penalize = TRUE, B = NULL, + B_null_list = NULL, fitted_model = NULL, refit = TRUE, test_kj = NULL, @@ -109,6 +144,7 @@ emuFit <- function(Y, use_fullmodel_info = FALSE, use_fullmodel_cov = FALSE, use_both_cov = FALSE, + match_row_names = TRUE, constraint_fn = pseudohuber_center, constraint_grad_fn = dpseudohuber_center_dx, constraint_param = 0.1, @@ -126,14 +162,20 @@ emuFit <- function(Y, inner_maxit = 25, max_step = 1, trackB = FALSE, - return_both_score_pvals = FALSE - - -) { + return_nullB = FALSE, + return_score_components = FALSE, + return_both_score_pvals = FALSE, + remove_zero_comparison_pvals = 0.01, + unobserved_taxon_error = TRUE) { # Record call call <- match.call(expand.dots = FALSE) + # confirm that input to verbose is valid + if (!(verbose %in% c(FALSE, TRUE, "development"))) { + stop('The argument "verbose" must be set to one of TRUE, FALSE, or "development".') + } + # check if Y is a phyloseq object if ("phyloseq" %in% class(Y)) { if (requireNamespace("phyloseq", quietly = TRUE)) { @@ -142,12 +184,32 @@ emuFit <- function(Y, } else { data <- data.frame(phyloseq::sample_data(Y)) X <- model.matrix(formula, data) + taxa_are_rows <- Y@otu_table@taxa_are_rows Y <- as.matrix(phyloseq::otu_table(Y)) + if (taxa_are_rows) { + Y <- t(Y) + } } } else { stop("You are trying to use a `phyloseq` data object or `phyloseq` helper function without having the `phyloseq` package installed. Please either install the package or use a standard data frame.") } - } else if ("data.frame" %in% class(Y)) { + + # check if Y is a TreeSummarizedExperiment object + } else if ("TreeSummarizedExperiment" %in% class(Y)) { + if (requireNamespace("SummarizedExperiment", quietly = TRUE)) { + if (is.null(assay_name) | is.null(formula)) { + stop("If Y is a `TreeSummarizedExperiment` object, make sure to include the assay_name and formula arguments.") + } + data <- as.data.frame(SummarizedExperiment::colData(Y)) + X <- model.matrix(formula, data) + Y <- as.data.frame(t(SummarizedExperiment::assay(Y, assay_name))) + } else { + stop("You are trying to use a `TreeSummarizedExperiment` data object or `TreeSummarizedExperiment` helper function without having the `SummarizedExperiment` package installed. Please either install the package or use a standard data frame.") + } + } + + # convert Y from a data.frame object to a matrix, even if it was extracted directly from `phyloseq` or `TreeSummarizedExperiment` + if ("data.frame" %in% class(Y)) { Y <- as.matrix(Y) if (!is.numeric(Y)) { stop("Y is a data frame that cannot be coerced to a numeric matrix. Please fix and try again.") @@ -168,20 +230,70 @@ covariates in formula must be provided.") } } + # check that if X and Y match in the row names + if (is.null(rownames(X)) || is.null(rownames(Y))){ + if (nrow(X) == nrow(Y)){ + if(match_row_names){ + if(is.null(rownames(X))){ + message("Row names are missing from the covariate matrix X. We will assume the rows are in the same order as in the response matrix Y. You are responsible for ensuring the order of your observations is the same in both matrices.") + } else { + message("Row names are missing from the response matrix Y. We will assume the rows are in the same order as in the covariate matrix X. You are responsible for ensuring the order of your observations is the same in both matrices.") + } + } + } else { + if(is.null(rownames(X))){ + stop("Row names are missing from the covariate matrix X, and the number of rows does not match the number of rows in the response matrix Y. Please resolve this issue before refitting the model.") + } else { + stop("Row names are missing from the response matrix Y, and the number of rows does not match the number of rows in the covariate matrix X. Please resolve this issue before refitting the model.") + } + } + } else{ + if(match_row_names){ + names_X <- rownames(X) + names_Y <- rownames(Y) + + #Checking if any row names are duplicated + if (any(duplicated(names_X))) stop("Covariate matrix X has duplicated row names. Please ensure all row names are unique before refitting the model.") + if (any(duplicated(names_Y))) stop("Response matrix Y has duplicated row names. Please ensure all row names are unique before refitting the model.") + + # Find common row names + common_names <- intersect(names_X, names_Y) + + if (length(common_names) < length(names_X) || length(common_names) < length(names_Y)) { + warning(sprintf("According to the rownames, there are observations that are missing either in the covariate matrix (X) and/or the response matrix (Y). We will subset to common rows only, resulting in %d samples.", length(common_names))) + + X <- X[common_names, , drop = FALSE] + Y <- Y[common_names, , drop = FALSE] + + } else if(all.equal(rownames(Y), rownames(X)) != TRUE){ + message("There is a different row ordering between the covariate matrix (X) and the response matrix (Y). Covariate data will be reordered to match response data.") + X <- X[rownames(Y), , drop = FALSE] + } + } else { + if(nrow(X) != nrow(Y)){ + stop("The number of rows does not match between the covariate matrix (X) and the response matrix (Y), and subsetting/matching by row name has been disabled. Please resolve this issue before refitting the model.") + } + } + } + if (min(rowSums(Y))==0) { stop("Some rows of Y consist entirely of zeroes, meaning that some samples have no observations. These samples must be excluded before fitting model.") } + if (unobserved_taxon_error) { + if (min(colSums(Y)) == 0) { + stop("Some columns of Y consist entirely of zeroes, meaning that some categories have zero counts for all samples. These + categories must be excluded before fitting the model.") + } + } + #check that cluster is correctly type if provided if(!is.null(cluster)){ - if(!is.numeric(cluster)){ - stop("If provided, argument 'cluster' must be a numeric vector.") - } if(length(cluster)!=nrow(Y)){ - stop("If provided, argument 'cluster' must be a numeric vector with + stop("If provided as a vector, argument 'cluster' must have length equal to n (the number of rows in Y).") - } + } if(length(unique(cluster)) == nrow(Y)){ warning("Number of unique values in 'cluster' equal to number of rows of Y; ignoring argument 'cluster'.") @@ -189,16 +301,56 @@ ignoring argument 'cluster'.") } } + # check that B_null_list is the correct length if provided + if (!is.null(B_null_list)) { + if (length(B_null_list) != nrow(test_kj)) { + warning("Length of 'B_null_list' is different than the number of tests specified in 'test_kj'. Ignoring object 'B_null_list'.") + B_null_list <- NULL + } + } + + # check for valid argument remove_zero_comparison_pvals + if (remove_zero_comparison_pvals != TRUE & remove_zero_comparison_pvals != FALSE) { + if (!(is.numeric(remove_zero_comparison_pvals) & remove_zero_comparison_pvals <= 1 & + remove_zero_comparison_pvals >= 0)) { + stop("Please set `remove_zero_comparison_pvals` to either TRUE, FALSE, or a numeric value between 0 and 1.") + } + } + + if (length(constraint_fn) == 1 & is.numeric(constraint_fn)) { + constraint_cat <- constraint_fn + constraint_fn <- function(x) {x[constraint_cat]} + constraint_grad_fn <- function(x) { + grad <- rep(0, length(x)) + grad[constraint_cat] <- 1 + return(grad) + } + constraint_param <- NA + } + n <- nrow(Y) J <- ncol(Y) p <- ncol(X) - X_cup <- X_cup_from_X(X,J) + if (is.null(colnames(X))) { + if (p > 1) { + colnames(X) <- c("Intercept", paste0("covariate_", 1:(ncol(X) - 1))) + } else { + colnames(X) <- "Intercept" + } + } + if (is.null(colnames(Y))) { + colnames(Y) <- paste0("category_", 1:ncol(Y)) + } + # check for zero-comparison parameters + zero_comparison_res <- zero_comparison_check(X = X, Y = Y) + + X_cup <- X_cup_from_X(X,J) if (is.logical(all.equal(constraint_fn, pseudohuber_center))) { if (all.equal(constraint_fn, pseudohuber_center)) { - if (verbose) message("Centering rows of B with pseudo-Huber smoothed median with smoothing parameter ", constraint_param, ".") + if (verbose %in% c(TRUE, "development")) message("Centering rows of B with pseudo-Huber smoothed median with smoothing parameter ", constraint_param, ".") stopifnot(!is.na(constraint_param)) @@ -238,10 +390,14 @@ and the corresponding gradient function to constraint_grad_fn.") maxit = maxit, max_step = max_step, tolerance = tolerance, - verbose = verbose, + verbose = (verbose == "development"), j_ref = j_ref) Y_test <- fitted_model$Y_augmented fitted_B <- fitted_model$B + converged_estimates <- fitted_model$convergence + if (!converged_estimates) { + warning("Optimization to estimate parameters did not converge. Try running again with a larger value of 'maxit'.") + } } else { fitted_model <- @@ -252,7 +408,8 @@ and the corresponding gradient function to constraint_grad_fn.") maxit = maxit, max_stepsize = max_step, tolerance = tolerance, - j_ref = j_ref) + j_ref = j_ref, + verbose = (verbose == "development")) fitted_B <- fitted_model Y_test <- Y } @@ -263,29 +420,54 @@ and the corresponding gradient function to constraint_grad_fn.") } if (!is.null(fitted_model)) { fitted_B <- fitted_model$B - Y_test <- Y_augmented <- fitted_model$Y_augmented - penalize <- fitted_model$penalized + if (penalize != fitted_model$penalized) { + stop("Your argument to `penalize` does not match the `penalize` argument within your `fitted_model` object. Please use the `penalize` argument that matches the `penalized` return object within your `fitted_model`.") + } + if (penalize) { + Y_test <- Y_augmented <- fitted_model$Y_augmented + } else { + Y_augmented <- NULL + Y_test <- Y + } if (!is.null(B)) { warning("B and fitted_model provided to emuFit; B ignored in favor of fitted_model.") } } else { fitted_B <- B + if (penalize) { + X_cup <- X_cup_from_X(X, J) + G <- get_G_for_augmentations(X, J, n, X_cup) + Y_test <- Y_augmented <- Y + + get_augmentations(X = X, G = G, Y = Y, B = fitted_B) + } else { + Y_augmented <- NULL + Y_test <- Y + } } } - - coefficients <- expand.grid(1:J, 2:p) - coefficients <- data.frame(j = coefficients[ , 1], - k = coefficients[ ,2]) - coefficients$estimate <- do.call(c, lapply(2:p, function(k) fitted_B[k,])) - coefficients$lower <- NA - coefficients$upper <- NA - coefficients$pval <- coefficients$score_stat <- NA - - + if (p == 1) { + message("You are running an intercept-only model. In this model the intercept beta_0^j is an unidentifiable combination of intercept for category j and the detection efficiency of category j. Therefore this parameter is not interpretable.") + + coefficients <- expand.grid(1:J, 1) + coefficients <- data.frame(j = coefficients[ , 1], + k = coefficients[ ,2]) + coefficients$estimate <- fitted_B[1, ] + coefficients$lower <- NA + coefficients$upper <- NA + coefficients$pval <- coefficients$score_stat <- NA + } else { + coefficients <- expand.grid(1:J, 2:p) + coefficients <- data.frame(j = coefficients[ , 1], + k = coefficients[ ,2]) + coefficients$estimate <- do.call(c, lapply(2:p, function(k) fitted_B[k,])) + coefficients$lower <- NA + coefficients$upper <- NA + coefficients$pval <- coefficients$score_stat <- NA + } if (compute_cis) { - if (verbose) { + if (verbose %in% c(TRUE, "development")) { message("Performing Wald tests and constructing CIs.") } @@ -297,7 +479,7 @@ and the corresponding gradient function to constraint_grad_fn.") constraint_fn = constraint_fn, constraint_grad_fn = constraint_grad_fn, nominal_coverage = 1 - alpha, - verbose = verbose, + verbose = (verbose == "development"), j_ref = j_ref, cluster = cluster) @@ -346,29 +528,52 @@ and the corresponding gradient function to constraint_grad_fn.") if (run_score_tests) { + score_test_hyperparams <- data.frame(u = rep(NA, nrow(test_kj)), + rho = NA, + tau = NA, + inner_maxit = NA, + gap = NA, + converged = NA) + + if (return_score_components) { + score_components <- vector(mode = "list", length = nrow(test_kj)) + } + if (return_both_score_pvals) { colnames(coefficients)[colnames(coefficients) == "pval"] <- "score_pval_full_info" colnames(coefficients)[colnames(coefficients) == "score_stat"] <- "score_stat_full_info" - coefficients$score_pval_null_info <- numeric(nrow(coefficients)) - coefficients$score_stat_null_info <- numeric(nrow(coefficients)) + coefficients$score_pval_null_info <- NA + coefficients$score_stat_null_info <- NA if (!use_fullmodel_info) { stop("If return_both_score_pvals = TRUE, use_fullmodel_info must be TRUE as well.") } } - + if (return_nullB) { + nullB_list <- vector(mode = "list", length = nrow(test_kj)) + } for(test_ind in 1:nrow(test_kj)) { - if (verbose) { + if (verbose %in% c(TRUE, "development")) { print(paste("Running score test ", test_ind, " of ", nrow(test_kj)," (row of B k = ", test_kj$k[test_ind], "; column of B j = ", test_kj$j[test_ind],").",sep = "")) } + B_to_use <- fitted_B + if (!is.null(B_null_list)) { + if (!is.null(B_null_list[[test_ind]])) { + B_to_use <- B_null_list[[test_ind]] + if (!(nrow(B_to_use) == nrow(fitted_B) & ncol(B_to_use) == ncol(fitted_B))) { + warning("'B_null_list' contains objects that are not the correct dimension for 'B'. The 'B_null_list' argument will be ignored.") + B_to_use <- fitted_B + } + } + } - test_result <- score_test(B = fitted_B, #B (MPLE) + test_result <- score_test(B = B_to_use, #B (MPLE or starting value if provided) Y = Y_test, #Y (with augmentations) X = X, #design matrix X_cup = X_cup, @@ -387,63 +592,78 @@ and the corresponding gradient function to constraint_grad_fn.") maxit = maxit, inner_maxit = inner_maxit, ntries = ntries, - verbose = verbose, + verbose = (verbose == "development"), trackB = trackB, I_inv = I_inv, Dy = Dy, return_both_score_pvals = return_both_score_pvals, cluster = cluster) - which_row <- which((as.numeric(coefficients$k) == as.numeric(test_kj$k[test_ind]))& - (as.numeric(coefficients$j) == as.numeric(test_kj$j[test_ind]))) - - if (!return_both_score_pvals) { - coefficients[which_row ,c("pval","score_stat")] <- - c(test_result$pval,test_result$score_stat) - } else { - coefficients[which_row ,c("score_pval_full_info","score_stat_full_info", - "score_pval_null_info","score_stat_null_info")] <- - c(test_result$pval,test_result$score_stat, - test_result$pval_null_info,test_result$score_stat_null_info) + if (return_score_components & !(is.null(test_result))) { + score_components[[test_ind]] <- test_result$score_pieces } - if (use_both_cov) { + if (is.null(test_result)) { + if (return_nullB) { + nullB_list[[test_ind]] <- NA + } + } else { - #adjustment factor from guo GEE paper (https://doi.org/10.1002/sim.2161) - alt_score_stat <- get_score_stat(Y = Y_test, - X_cup = X_cup, - X = X, - B = test_result$null_B, - k_constr = test_kj$k[test_ind], - j_constr = test_kj$j[test_ind], - constraint_grad_fn = constraint_grad_fn, - indexes_to_remove = (j_ref - 1)*p + 1:p, - j_ref = j_ref, - J = J, - n = n, - p = p, - I_inv=I_inv, - Dy = just_wald_things$Dy, - cluster = cluster) + score_test_hyperparams[test_ind, ] <- + c(test_result$u, test_result$rho, test_result$tau, test_result$inner_maxit, + test_result$gap, test_result$convergence) + if (return_nullB) { + null_B <- test_result$null_B + for (k in 1:p) { + null_B[k, ] <- null_B[k, ] - constraint_fn(null_B[k, ]) + } + nullB_list[[test_ind]] <- null_B + } which_row <- which((as.numeric(coefficients$k) == as.numeric(test_kj$k[test_ind]))& (as.numeric(coefficients$j) == as.numeric(test_kj$j[test_ind]))) - coefficients[which_row, c("score_fullcov_p")] <- pchisq(alt_score_stat,1, - lower.tail = FALSE) - } - + + if (!return_both_score_pvals) { + coefficients[which_row ,c("pval","score_stat")] <- + c(test_result$pval,test_result$score_stat) + } else { + coefficients[which_row ,c("score_pval_full_info","score_stat_full_info", + "score_pval_null_info","score_stat_null_info")] <- + c(test_result$pval,test_result$score_stat, + test_result$pval_null_info,test_result$score_stat_null_info) + } + + if (use_both_cov) { + + #adjustment factor from guo GEE paper (https://doi.org/10.1002/sim.2161) + alt_score_stat <- get_score_stat(Y = Y_test, + X_cup = X_cup, + X = X, + B = test_result$null_B, + k_constr = test_kj$k[test_ind], + j_constr = test_kj$j[test_ind], + constraint_grad_fn = constraint_grad_fn, + indexes_to_remove = (j_ref - 1)*p + 1:p, + j_ref = j_ref, + J = J, + n = n, + p = p, + I_inv=I_inv, + Dy = just_wald_things$Dy, + cluster = cluster)$score_stat + + + which_row <- which((as.numeric(coefficients$k) == as.numeric(test_kj$k[test_ind]))& + (as.numeric(coefficients$j) == as.numeric(test_kj$j[test_ind]))) + coefficients[which_row, c("score_fullcov_p")] <- pchisq(alt_score_stat,1, + lower.tail = FALSE) + } + } } } - if (is.null(colnames(X))) { - colnames(X) <- c("Intercept", paste0("covariate_", 1:(ncol(X) - 1))) - } - if (is.null(colnames(Y))) { - colnames(Y) <- paste0("category_", 1:ncol(Y)) - } - if (!is.null(colnames(X))) { if (length(unique(colnames(X))) == ncol(X)) { k_to_covariates <- data.frame(k = 1:p, @@ -472,62 +692,56 @@ and the corresponding gradient function to constraint_grad_fn.") } if (!compute_cis) { - coefficients <- + coef_df <- + cbind(data.frame(covariate = coefficients$covariate, + category = coefficients$category, + category_num = coefficients$j), + coefficients[ , c("estimate","lower","upper")]) + } else { + coef_df <- cbind(data.frame(covariate = coefficients$covariate, category = coefficients$category, category_num = coefficients$j), - coefficients[ , c("estimate","lower","upper","score_stat","pval")]) + coefficients[ , c("estimate","se", "lower","upper")]) + } + if (use_both_cov) { + coef_df <- cbind(coef_df, coefficients[ , c("score_stat","pval","score_fullcov_p")]) } else { - if (!return_wald_p) { - coefficients <- - cbind(data.frame(covariate = coefficients$covariate, - category = coefficients$category, - category_num = coefficients$j), - coefficients[ , c("estimate","se","lower","upper","score_stat","pval")]) + if (return_both_score_pvals) { + coef_df <- cbind(coef_df, coefficients[ , c("score_stat_full_info", + "score_pval_full_info", + "score_stat_null_info", + "score_pval_null_info")]) } else { - if (use_both_cov) { - coefficients <- - cbind(data.frame(covariate = coefficients$covariate, - category = coefficients$category, - category_num = coefficients$j), - coefficients[ , c("estimate","se","lower","upper","score_stat","pval", - "wald_p","score_fullcov_p")]) - } else { - if (return_both_score_pvals) { - coefficients <- - cbind(data.frame(covariate = coefficients$covariate, - category = coefficients$category, - category_num = coefficients$j), - coefficients[ , c("estimate","se","lower","upper", - "score_stat_full_info", - "score_pval_full_info", - "score_stat_null_info", - "score_pval_null_info", - "wald_p")]) - - } else { - coefficients <- - cbind(data.frame(covariate = coefficients$covariate, - category = coefficients$category, - category_num = coefficients$j), - coefficients[ , c("estimate","se","lower","upper","score_stat","pval","wald_p")]) - } - } + coef_df <- cbind(coef_df, coefficients[ , c("score_stat","pval")]) } } + if (return_wald_p) { + coef_df$wald_p <- coefficients[, "wald_p"] + } + coefficients <- coef_df if (penalize) { - Y_augmented <- fitted_model$Y_augmented - } else { if (!is.null(fitted_model)) { Y_augmented <- fitted_model$Y_augmented - } else { - Y_augmented <- NULL } + } else { + # set Y_augmented to NUll because without penalty there is no Y augmentation + Y_augmented <- NULL } if (!is.null(fitted_model)) { - B <- fitted_model$B + if (penalize) { + B <- fitted_model$B + } else { + B <- fitted_model + } + } + + if (penalize) { + z_hat <- log(rowSums(Y_augmented)) - log(rowSums(exp(X %*% B))) + } else { + z_hat <- log(rowSums(Y)) - log(rowSums(exp(X %*% B))) } if (is.null(just_wald_things)) { @@ -546,14 +760,54 @@ and the corresponding gradient function to constraint_grad_fn.") rownames(B) <- colnames(X) } - return(structure(list("call" = call, - "coef" = coefficients, - "B" = B, - "penalized" = penalize, - "Y_augmented" = Y_augmented, - "I" = I, - "Dy" = Dy, - "cluster" = cluster), class = "emuFit")) + if (!is.null(zero_comparison_res)) { + coefficients <- dplyr::full_join(coefficients, zero_comparison_res, + by = c("covariate", "category")) + coefficients$zero_comparison[is.na(coefficients$zero_comparison)] <- FALSE + + if (remove_zero_comparison_pvals == TRUE | is.numeric(remove_zero_comparison_pvals)) { + pval_cols <- which(grepl("pval", names(coefficients))) + for (col in pval_cols) { + if (remove_zero_comparison_pvals == TRUE) { + coefficients[coefficients$zero_comparison, col] <- NA + } else { + ind <- ifelse(is.na(coefficients[, col]), FALSE, + coefficients[, col] <= remove_zero_comparison_pvals & coefficients$zero_comparison) + coefficients[ind, col] <- NA + } + } + } + } + + results <- list("call" = call, + "coef" = coefficients, + "B" = B, + "penalized" = penalize, + "Y_augmented" = Y_augmented, + "z_hat" = z_hat, + "I" = I, + "Dy" = Dy, + "cluster" = cluster) + + if (refit & penalize) { + results$estimation_converged <- converged_estimates + } + if (run_score_tests & return_nullB) { + results$null_B <- nullB_list + } + if (run_score_tests) { + results$score_test_hyperparams <- score_test_hyperparams + if (sum(score_test_hyperparams$converged != "converged") > 0) { + unconverged_test_kj <- test_kj[which(score_test_hyperparams$converged != "converged"), ] + results$null_estimation_unconverged <- unconverged_test_kj + warning("Optimization for estimation under the null for robust score tests failed to converge for some tests. See 'null_estimation_unconverged' within the returned emuFit object for which tests are affected by this.") + } + } + if (run_score_tests & return_score_components) { + results$score_components <- score_components + } + + return(structure(results, class = "emuFit")) } diff --git a/R/emuFit_micro.R b/R/emuFit_micro.R index 46670fb..a326536 100644 --- a/R/emuFit_micro.R +++ b/R/emuFit_micro.R @@ -67,6 +67,14 @@ emuFit_micro <- Y_start[i,] <- Y_start[i,] - mean(Y_start[i,]) } B <- matrix(nrow = p,ncol = J) + + ##Checking if design matrix is rank-deficient and soln_mat can be found + if (qr(t(X)%*%X)$rank < ncol(X)){ + stop("Design matrix X inputted for the model is rank-deficient, preventing proper model fitting. + This might be due to multicollinearity, overparameterization, or redundant factor levels included in covariates. + Consider removing highly correlated covariates or adjusting factor levels to ensure a full-rank design. \n") + } + soln_mat <- qr.solve(t(X)%*%X,t(X)) for(j in 1:J){ B[,j] <- as.numeric(as.matrix(soln_mat%*%Y_start[,j,drop = FALSE])) @@ -125,6 +133,28 @@ emuFit_micro <- B_diff <- Inf while(!converged){ old_B <- B + if(!is.null(j_ref)){ + shifty <- optim(rep(0,p),function(x) lil_ll(x, + B = B, + p = p, + X = X, + Y = Y, + J = J, + j_ref = j_ref), + method = "BFGS") + + shift <- shifty$par + + # print(signif(shifty$par,2)) + + for(k in 1:p){ + B[k,-j_ref] <- B[k,-j_ref] + shifty$par[k] + } + # + # } + + z <- update_z(X = X, Y = Y, B = B) + } for(j in loop_js){ # print(j) # llj <- function(x){ diff --git a/R/emuFit_micro_penalized.R b/R/emuFit_micro_penalized.R index f654f50..13a6d4b 100644 --- a/R/emuFit_micro_penalized.R +++ b/R/emuFit_micro_penalized.R @@ -127,15 +127,18 @@ maintained only for testing purposes.") if(B_diff < tolerance){ converged <- TRUE + actually_converged <- TRUE } if(counter>maxit){ converged <- TRUE + actually_converged <- FALSE } counter <- counter + 1 } return(list("Y_augmented" = Y_augmented, - "B" = fitted_model)) + "B" = fitted_model, + "convergence" = actually_converged)) } diff --git a/R/fit_null_repar.R b/R/fit_null_repar.R new file mode 100644 index 0000000..b13a9a9 --- /dev/null +++ b/R/fit_null_repar.R @@ -0,0 +1,231 @@ +#' fits model with B_kj constrained to equal g(B_k) for constraint fn g +#' +#' @param B description +#' @param Y Y (with augmentations) +#' @param X design matrix +#' @param k_constr row index of B to constrain +#' @param j_constr col index of B to constrain +#' @param j_ref column index of convenience constraint +#' @param constraint_fn constraint function +#' @param constraint_grad_fn gradient of constraint fn +#' @param tolerance convergence tolerance. Once no element of B changes by +#' this amount or more in an iteration, convergence is declared. +#' @param maxit maximum iterations +#' @param verbose shout at you? +#' @param trackB track value of beta across iterations and return? +#' @param starting_stepsize stepsize at which backtracking line search begins. Default is 0.5. +#' @param do_shift perform optimization of row-wise location of B in between rounds of coordinate descent? +#' Default is TRUE. +#' +#' @returns A list containing elements `B`, `k_constr`, `j_constr`, `niter` +#' `converged`, `dll_dB`, and `Bs`. `B` is a matrix containing parameter estimates +#' under the null (obtained by maximum likelihood on augmented observations Y), +#' 'k_constr', and 'j_constr' give row and column indexes of the parameter +#' fixed to be equal to the constraint function g() under the null. `niter` is a +#' scalar giving total number of outer iterations used to fit the null model. +#' 'converged' is an indicator that convergence was reached (TRUE) -- if optimization +#' stopped because the iteration limit was reached, `convergence` will have value FALSE. +#' `dll_dB` is the derivative of the Poisson log likelihood with respect to the +#' parameters in the reparametrized null model, and +#' 'Bs' is a data frame containing values of B by iteration if trackB was set +#' equal to TRUE (otherwise it contains a NULL value). +#' +fit_null_repar <- function(B, + Y, + X, + k_constr, + j_constr, + j_ref, + constraint_fn, + constraint_grad_fn, + tolerance = 1e-3, + maxit = 100, + verbose = FALSE, + trackB = FALSE, + starting_stepsize = 0.5, + method = "fisher", + do_shift = TRUE, + shift_tolerance = 1e-3 +) { + + J <- ncol(Y) + n <- nrow(Y) + p <- ncol(X) + + if(j_ref == j_constr){ + stop("Taxon j chosen for convenience constraint may not be null taxon.") + } + + if(j_ref != J){ + stop("fit_null_repar has not yet been tested with j_ref != J.") + } + + #change id. constr. by subtracting approp. col of B from others + for(k in 1:p) { + B[k,] <- B[k,] - B[k,j_ref] + } + + # impose null constraint -- only valid for permutation-invariant constraint functions + B[k_constr,j_constr] <- constraint_fn(B[k_constr,-j_constr]) + + #update z + z <- update_z(X = X, Y = Y, B = B) + + if(trackB){ + Bs <- vector(mode = "list",length = maxit + 1) + Bs[[1]] <- B + } else{ + Bs <- NULL + } + + #evaluate log likelihood + log_means <- X%*%B + matrix(z, ncol = 1)%*%matrix(1,ncol = J,nrow = 1) + ll <- sum(Y*log_means - exp(log_means)) + lls <- ll + + prev_ll <- -Inf + Bs <- list(B) + + #loop over j except j_constr and j_ref + loop_j <- setdiff(1:J,c(j_constr,j_ref)) + + #initiate prev_B at an arbitrary value that will not trigger convergence + stop_iteration <- FALSE + + iter <- 1 + while(!stop_iteration){ + #store previous value of B + prev_B <- B + + if(do_shift){ + lil_ll_rep <- function(shifts,B,X,Y){ + for(k in 1:p){ + B[k,-j_ref] <- B[k,-j_ref] + shifts[k] + } + B[k_constr,j_constr] <- constraint_fn(B[k_constr,-j_constr]) + z <- update_z(X = X, Y = Y, B = B) + log_means <- X%*%B + matrix(z, ncol = 1)%*%matrix(1,ncol = J,nrow = 1) + + return(-sum(Y*log_means - exp(log_means))) + } + + shifty <- optim(rep(0,p),function(x) lil_ll_rep(x,B = B, X = X, Y = Y), + method = "BFGS") + + shift <- shifty$par + + # if(max(abs(shift))1){ + old_dl <- deriv_long + } + deriv_long <- do.call(c,lapply(1:p,function(k) deriv_mat[k,])) + + # + # if(iter ==1){ + # plot(asinh(deriv_long),col = "red",pch = 20) + # } else{ + # points(asinh(old_dl), col = "grey",pch = 20) + # points(asinh(deriv_long),col = "red",pch = 20) + # } + } + iter <- iter + 1 + + if(iter > maxit){ + stop_iteration <- TRUE + converged <- FALSE + } + + if(verbose){ + message("Maximum absolute difference in B: ",signif(max(abs(B - prev_B)),3)) + } + # hist(as.numeric(log(abs(B - prev_B))/log(10)), + # breaks = 100,xlab = "Log 10 absolute difference") + if(max(abs(B - prev_B))% +#' @importFrom rlang .data +#' @importFrom stats reorder +#' @importFrom stats setNames +#' +#' @return Object of class \code{ggplot}. Plot of \code{radEmu} model fit with 95% confidence intervals. +#' +#' @examples +#' data(wirbel_sample) +#' data(wirbel_otu) +#' +#' subset_studies <- which(wirbel_sample$Study %in% c("FR-CRC", "US-CRC", "AT-CRC")) +#' +#' chosen_genera <- c("Eubacterium", "Faecalibacterium", "Fusobacterium", "Porphyromonas") +#' +#' mOTU_names <- colnames(wirbel_otu) +#' mOTU_name_df <- data.frame(name = mOTU_names) %>% +#' dplyr::mutate(base_name = stringr::str_remove(mOTU_names, "unknown ") %>% +#' stringr::str_remove("uncultured ")) %>% +#' dplyr::mutate(genus_name = stringr::word(base_name, 1)) +#' +#' restricted_mOTU_names <- mOTU_name_df %>% +#' dplyr::filter(genus_name %in% chosen_genera) %>% +#' dplyr::pull(name) +#' +#' small_Y <- wirbel_otu[subset_studies, restricted_mOTU_names] +#' category_to_rm <- which(colSums(small_Y) == 0) +#' +#' small_sample <- wirbel_sample[subset_studies, ] +#' +#' ch_fit <- emuFit(formula = ~ Group + Study, +#' data = small_sample, +#' Y = small_Y, +#' run_score_tests = FALSE) +#' +#' plot_key <- list(p1 = c("Control" = "GroupCTR"), +#' p2 = c("FR-Control" = "StudyFR-CRC", +#' "US-Control" = "StudyUS-CRC")) +#' +#' out <- plot(x = ch_fit, +#' plot_key = plot_key, +#' display_taxon_names = FALSE) +#' +#' out$plots$p1 +#' out$plots$p2 +#' @export + +plot.emuFit <- function(x, + plot_key = NULL, + title = NULL, + taxon_names = NULL, + display_taxon_names = TRUE, + data_only = FALSE, ...) { + mod <- x + + # confirm that values in plot_key are actually in the coefficient table + if (!all(unlist(plot_key) %in% unique(mod$coef$covariate))) { + stop("At least one of the coefficient names included in the plot key is not in + the covariate column of the coefficient table of the radEmu() output.") + } + + # confirm that no coefficient is listed multiple times + if (length(unlist(plot_key)) != length(unique(unlist(plot_key)))) { + stop("One of the coefficient names is included in the plot key multiple times, + however each coefficient can only be included in one plot.") + } + + # determine which levels are not included in the plot key + remaining_variables <- setdiff(unique(mod$coef$covariate), unlist(plot_key)) + + # construct a list of variables that should be plotted together + plot_key_default <- append(plot_key, + lapply(as.list(remaining_variables), + function(vec){ + setNames(vec, vec) + })) + + plot_key_default <- lapply(plot_key_default, function(vec){ + if (is.null(names(vec))) { + setNames(vec, vec) + } else { + setNames(vec, ifelse(names(vec) != "", names(vec), vec)) + } + + }) + + # for each covariate, determine which it should be included in + covariate_groups <- sapply(plot_key_default, function(key){ + mod$coef$covariate %in% key + }) + + # add a variable to identify which plot each coefficient will be plotted on + mod$coef$plot_key <- apply(covariate_groups, 1, which) + + # rename variable levels as given in the plot key + new_variable_names <- names(unlist(unname(plot_key_default))) + old_variable_names <- unname(unlist(unname(plot_key_default))) + mod$coef <- mod$coef %>% + dplyr::mutate(covariate = factor(covariate, + levels = old_variable_names, + labels = new_variable_names)) + + # now separate output coefficients into a list for each plot + coef_list <- mod$coef %>% + split(.$plot_key) + + # match coefficient names using user-provided "taxon_names" data frame + coef_list_renamed <- lapply(coef_list, function(coef_subset){ + if (!is.null(taxon_names)) { + dplyr::left_join(coef_subset, taxon_names, by = "category") %>% + dplyr::arrange(covariate, estimate) %>% + dplyr::mutate(cat_small = reorder(cat_small, estimate)) + } else { + coef_subset %>% + dplyr::mutate(cat_small = category) %>% + dplyr::arrange(covariate, estimate) %>% + dplyr::mutate(cat_small = factor(cat_small, + levels = unique(cat_small), + labels = unique(cat_small))) + } + }) + + # return plots + if (!data_only) { + + p_list <- lapply(coef_list_renamed, function(coef_subset){ + p <- ggplot2::ggplot(coef_subset) + + ggplot2::geom_point(ggplot2::aes(x = estimate, + y = as.character(cat_small), + color = covariate, + group = covariate), + position = ggplot2::position_dodge2(width = 0.5), + size = 2) + + ggplot2::geom_errorbar(ggplot2::aes(y = cat_small, + xmin = lower, + xmax = upper, + color = covariate, + group = covariate), + position = ggplot2::position_dodge2(width = 0.5), + width = 0.5) + + ggplot2::geom_vline(xintercept = 0, alpha = 0.5) + + ggplot2::theme_bw() + + ggplot2::labs(title = title) + + ggplot2::guides(color = ggplot2::guide_legend(title = "Comparison")) + + ggplot2::theme(legend.position = "bottom") + + ggplot2::labs(y = "Category") + + ggplot2::labs(x = "Estimate") + + if (!display_taxon_names) { + p <- p + + ggplot2::theme(axis.text.y = ggplot2::element_blank(), + axis.ticks.y = ggplot2::element_blank()) + } + + p + }) + + p_list <- setNames(p_list, paste0("p", 1:length(p_list))) + + return(list(plots = invisible(p_list), + data = dplyr::bind_rows(coef_list_renamed))) + + } else { + + # return data + return(list(plots = NULL, + data = dplyr::bind_rows(coef_list_renamed))) + } +} + + diff --git a/R/profile_ll.R b/R/profile_ll.R new file mode 100644 index 0000000..75e9d60 --- /dev/null +++ b/R/profile_ll.R @@ -0,0 +1,7 @@ + +#compute log likelihood after profiling out z +profile_ll <- function(X,Y,B){ + z <- update_z(Y = Y, X = X, B = B) + log_means <- X%*%B + matrix(z,ncol = 1)%*%matrix(1,nrow = 1, ncol = ncol(Y)) + return(sum(Y*log_means - exp(log_means))) +} \ No newline at end of file diff --git a/R/score_test.R b/R/score_test.R index 7a586d5..d17d211 100644 --- a/R/score_test.R +++ b/R/score_test.R @@ -1,6 +1,7 @@ #' Run robust score test -#' @param B value of coefficient matrix (p x J) returned by full model fit +#' @param B value of coefficient matrix (p x J) returned by full model fit or value of coefficient +#' matrix to start null estimation at given as input to emuFit #' @param Y an n x J matrix or dataframe of *augmented* nonnegative observations (i.e., #' observations Y plus augmentations from last iteration of maximum penalized likelihood estimation #' for full model) @@ -70,13 +71,16 @@ #' Y are treated as independent. #' #' @return A list containing elements 'score_stat', 'pval', 'log_pval','niter', -#' 'convergence', 'gap', 'u', 'rho', 'null_B', and 'Bs'. 'score_stat' gives the +#' 'convergence', 'gap', 'u', 'rho', 'tau', 'inner_maxit', 'null_B', and 'Bs'. 'score_stat' gives the #' value of the robust score statistic for H_0: B_{k_constr,j_constr} = g(B_{k_constr}). #' 'pval' and 'log_pval' are the p-value (on natural and log scales) corresponding to #' the score statistic (log_pval may be useful when the p-value is very close to zero). #' 'gap' is the final value of g(B_{k_constr}) - B_{k_constr, j_constr} obtained in #' optimization under the null. 'u' and 'rho' are final values of augmented -#' Lagrangian parameters returned by null fitting algorithm. 'null_B' is the value of +#' Lagrangian parameters returned by null fitting algorithm. 'tau' is the final value of 'tau' that +#' is used to update the 'rho' values and 'inner_maxit' is the final maximum number of iterations for +#' the inner optimization loop in optimization under the null, in which B and z parameter values are +#' maximized for specific 'u' and 'rho' parameters. 'null_B' is the value of #' B returned but the null fitting algorithm. 'Bs' is by default NULL; if trackB = TRUE, #' 'Bs is a data frame containing values of B by outcome category, covariate, and #' iteration. @@ -123,54 +127,54 @@ score_test <- function(B, #B (MPLE) good_enough_fit <- FALSE while(!accept_try){ #fit under null - constrained_fit <- try(fit_null(B = B, #B (MPLE) - Y = Y, #Y (with augmentations) - X = X, #design matrix - X_cup = X_cup, - k_constr = k_constr, #row index of B to constrain - j_constr = j_constr, #col index of B to constrain - constraint_fn = constraint_fn, #constraint function - constraint_grad_fn = constraint_grad_fn, #gradient of constraint fn - # constraint_hess_fn = constraint_hess_fn, - rho_init = rho_init, - tau = tau, - kappa = kappa, - B_tol = B_tol, - inner_tol = inner_tol, - constraint_tol = constraint_tol, - j_ref = j_ref, - c1 = c1, - maxit = maxit, - inner_maxit = inner_maxit, - verbose = verbose, - trackB = trackB - # I = I, - # Dy = Dy - )) - if(inherits(constrained_fit,"try-error")){ - accept_try <- FALSE - } else{ - if((abs(constrained_fit$gap) <= constraint_tol) & - (constrained_fit$niter < maxit)){ - accept_try <- TRUE - good_enough_fit <- TRUE + constrained_fit <- try(fit_null(B = B, #B (MPLE) + Y = Y, #Y (with augmentations) + X = X, #design matrix + X_cup = X_cup, + k_constr = k_constr, #row index of B to constrain + j_constr = j_constr, #col index of B to constrain + constraint_fn = constraint_fn, #constraint function + constraint_grad_fn = constraint_grad_fn, #gradient of constraint fn + # constraint_hess_fn = constraint_hess_fn, + rho_init = rho_init, + tau = tau, + kappa = kappa, + B_tol = B_tol, + inner_tol = inner_tol, + constraint_tol = constraint_tol, + j_ref = j_ref, + c1 = c1, + maxit = maxit, + inner_maxit = inner_maxit, + verbose = verbose, + trackB = trackB + # I = I, + # Dy = Dy + )) + if(inherits(constrained_fit,"try-error")){ + accept_try <- FALSE } else{ - tau <- tau^(3/4) - inner_maxit <- 2*inner_maxit - message("Constrained optimization failed to converge within iteration limit; + if((abs(constrained_fit$gap) <= constraint_tol) & + (constrained_fit$niter < maxit)){ + accept_try <- TRUE + good_enough_fit <- TRUE + } else{ + tau <- tau^(3/4) + inner_maxit <- 2*inner_maxit + message("Constrained optimization failed to converge within iteration limit; retrying with smaller penalty scaling parameter tau and larger inner_maxit.") + } + } + tries_so_far <- tries_so_far + 1 + if(tries_so_far == ntries){ + accept_try <- TRUE } - } - tries_so_far <- tries_so_far + 1 - if(tries_so_far == ntries){ - accept_try <- TRUE - } } if(!good_enough_fit){ warning("Optimization for null fit with k = ",k_constr," and j = ",j_constr," failed to converge across ", ntries, ifelse(ntries>1," attempts."," attempt.")) -} + } B <- constrained_fit$B z <- update_z(Y,X,B) p <- ncol(X) @@ -182,7 +186,7 @@ retrying with smaller penalty scaling parameter tau and larger inner_maxit.") #indexes in long format corresponding to the j_constr-th col of B #get score stat indexes_to_remove <- (j_ref - 1)*p + 1:p - score_stat <- + score_res <- try( get_score_stat(Y = Y, X_cup = X_cup, X = X, @@ -197,12 +201,21 @@ retrying with smaller penalty scaling parameter tau and larger inner_maxit.") p = p, I_inv = I_inv, Dy = Dy, - cluster = cluster) + cluster = cluster)) + if (inherits(score_res, "try-error")) { + score_stat <- score_res + } else { + score_stat <- score_res$score_stat + } if(!return_both_score_pvals){ #typically we want only one score p-value #(using only one version of information matrix) - + if (inherits(score_stat, "try-error")) { + warning("score statistic for test of k = ", k_constr, " and j = ", j_constr, " cannot be computed, likely because the information matrix is computationally singular.") + score_stat <- NA + } return(list("score_stat" = score_stat, + "score_pieces" = score_res, "pval" = pchisq(score_stat,1,lower.tail = FALSE), "log_pval" = pchisq(score_stat,1,lower.tail = FALSE, log.p = TRUE), "niter" = constrained_fit$niter, @@ -211,13 +224,15 @@ retrying with smaller penalty scaling parameter tau and larger inner_maxit.") "gap" = constrained_fit$gap, "u" = constrained_fit$u, "rho" = constrained_fit$rho, + "tau" = tau, + "inner_maxit" = inner_maxit, "null_B" = constrained_fit$B, # "score_stats" = constrained_fit$score_stats, "Bs" = constrained_fit$Bs)) } else{ #for simulations -- if we want to return both the score p-value using #information from full model fit and from null model - score_stat_with_null_info <- + score_res_with_null_info <- get_score_stat(Y = Y, X_cup = X_cup, X = X, @@ -232,13 +247,23 @@ retrying with smaller penalty scaling parameter tau and larger inner_maxit.") p = p, I_inv = NULL, Dy = Dy) - + if (inherits(score_res_with_null_info, "try-error")) { + score_stat_with_null_info <- score_res_with_null_info + } else { + score_stat_with_null_info <- score_res_with_null_info$score_stat + } score_stat_with_null_info <- score_stat_with_null_info + if (inherits(score_stat_with_null_info, "try-error")) { + warning("one of the score statistics for test of k = ", k_constr, " and j = ", j_constr, " cannot be computed, likely because the information matrix is computationally singular.") + score_stat_with_null_info <- NA + } return(list("score_stat" = score_stat, + "score_pieces" = score_res, "pval" = pchisq(score_stat,1,lower.tail = FALSE), "log_pval" = pchisq(score_stat,1,lower.tail = FALSE, log.p = TRUE), "score_stat_null_info" = score_stat_with_null_info, + "score_pieces_null_info" = score_res_with_null_info, "pval_null_info" = pchisq(score_stat_with_null_info,1,lower.tail = FALSE), "log_pval_null_info" = pchisq(score_stat_with_null_info,1,lower.tail = FALSE,log.p = TRUE), "niter" = constrained_fit$niter, @@ -247,6 +272,8 @@ retrying with smaller penalty scaling parameter tau and larger inner_maxit.") "gap" = constrained_fit$gap, "u" = constrained_fit$u, "rho" = constrained_fit$rho, + "tau" = tau, + "inner_maxit" = inner_maxit, "null_B" = constrained_fit$B, # "score_stats" = constrained_fit$score_stats, "Bs" = constrained_fit$Bs)) diff --git a/R/simulate_data.R b/R/simulate_data.R index 01a4f10..f363375 100644 --- a/R/simulate_data.R +++ b/R/simulate_data.R @@ -1,45 +1,103 @@ +#' Data simulation function +#' +#' Function to simulate data for simulations in Clausen & Willis (2024) and for the cluster vignette +#' +#' @param n Number of samples +#' @param J Number of categories +#' @param b0 Intercept parameter vector +#' @param b1 Covariate paramter vector +#' @param distn Distribution to simulate from, either "Poisson" or "ZINB" +#' @param zinb_size Size parameter for negative binomial draw for ZINB data +#' @param zinb_zero_prop Proportion of zeros for ZINB data +#' @param mean_z Parameter controlling the mean of the sample-specific effects. +#' @param X Optional design matrix, this must have two columns and n rows. +#' @param B Optional B matrix, if p is not equal to 2 +#' @param cluster Optional cluster vector, this must have n elements. +#' +#' @return \code{Y}. A \code{n times J} dimension matrix of simulated response counts. +#' #' @importFrom stats rnorm rbinom rpois rnbinom -#function to simulate data for simulations in Clausen & Willis (2024) +#' simulate_data <- function(n, J, - b0, - b1, + b0 = NULL, + b1 = NULL, distn, zinb_size = NULL, zinb_zero_prop = NULL, - mean_count_before_ZI) { + mean_z, + X = NULL, + B = NULL, + cluster = NULL) { + + if (is.null(X)) { + X <- cbind(1,rep(c(0,1),each = n/2)) + } + if (!is.null(b0) & !is.null(b1)) { + B <- rbind(b0,b1) + if (nrow(B) != ncol(X)) { + stop("You've input b0 and b1 but your X matrix does not have 2 columns. Please use the B argument when your design matrix does not have 2 columns.") + } + } else if (is.null(b0) | is.null(b1)) { + if (is.null(B)) { + stop("Please input either parameter vectors b0 and b1, or parameter matrix B.") + } + } - X <- cbind(1,rep(c(0,1),each = n/2)) - B <- rbind(b0,b1) log_means <- do.call(cbind, lapply(1:J, - function(j) X%*%B[,j,drop = FALSE])) - - row_means <- rowSums(exp(log_means))/J - - z <- sapply(row_means,function(x) log(mean_count_before_ZI) - log(x) + stats::rnorm(1)) + function(j) X %*% B[,j,drop = FALSE])) + + z <- stats::rnorm(n = n, mean = mean_z, sd = 1) + Y <- matrix(0, ncol = J, nrow = n) - + for(i in 1:n){ log_means[i,] <- log_means[i,] + z[i] } - - for(i in 1:n){ - accepted <- FALSE - while(!accepted){ - for(j in 1:J){ - if(distn == "Poisson"){ - Y[i,j] <- stats::rpois(1,lambda = exp(log_means[i,j])) - } - if(distn == "ZINB"){ - Y[i,j] <- stats::rnbinom(1,mu = exp(log_means[i,j]), size= zinb_size)* - (1 - stats::rbinom(1,1,prob = zinb_zero_prop)) - } - if(sum(Y[i,])>0){ - accepted <- TRUE + + if (is.null(cluster)) { + for(i in 1:n){ + accepted <- FALSE + while(!accepted){ + for(j in 1:J){ + if(distn == "Poisson"){ + Y[i,j] <- stats::rpois(1,lambda = exp(log_means[i,j])) + } + if(distn == "ZINB"){ + Y[i,j] <- stats::rnbinom(1,mu = exp(log_means[i,j]), size= zinb_size)* + (1 - stats::rbinom(1,1,prob = zinb_zero_prop)) + } + if(sum(Y[i,])>0){ + accepted <- TRUE + } + } } + } + } else { + cluster_effs <- lapply(1:length(unique(cluster)), function(i) log(matrix(stats::rexp(2*J), nrow = 2))) + for(i in 1:n){ + accepted <- FALSE + while(!accepted){ + for(j in 1:J){ + if(distn == "Poisson"){ + Y[i,j] <- stats::rpois(1,lambda = exp(log_means[i,j] + + cluster_effs[[cluster[i]]][, j])) + } + if(distn == "ZINB"){ + Y[i,j] <- stats::rnbinom(1, + mu = exp(log_means[i,j] + + cluster_effs[[cluster[i]]][, j]), + size= zinb_size)* + (1 - stats::rbinom(1,1,prob = zinb_zero_prop)) + } + if(sum(Y[i,])>0){ + accepted <- TRUE + } + } } } } + return(Y) } diff --git a/R/zero_comparison_check.R b/R/zero_comparison_check.R new file mode 100644 index 0000000..502ad96 --- /dev/null +++ b/R/zero_comparison_check.R @@ -0,0 +1,120 @@ +# function to check for zero comparison parameters +zero_comparison_check <- function(X, Y) { + + # X is a matrix from model.matrix and has "assign" attribute + if ("assign" %in% names(attributes(X))) { + + col_assign <- attr(X, "assign") + col_shared <- sapply(col_assign, function(x) {sum(col_assign == x)}) + col_kept <- col_assign[col_shared > 1] + + if (max(col_shared) > 1) { + + zero_comp_dat <- data.frame(covariate = NULL, category = NULL, + zero_comparison = NULL) + + cat_covs <- unique(col_kept) + for (col in cat_covs) { + base_X <- X[, col_assign == col] + + # get indices for each group + n_groups <- ncol(base_X) + 1 + group_ind <- vector(mode = "list", length = n_groups) + group_ind[[1]] <- which(rowSums(base_X) == 0) + + for (i in 1:ncol(base_X)) { + group_ind[[i + 1]] <- which(base_X[, i] == 1) + } + + # get Y sums for each group + J <- ncol(Y) + group_counts <- matrix(NA, nrow = n_groups, ncol = J) + for (i in 1:n_groups) { + group_counts[i, ] <- colSums(Y[group_ind[[i]], ]) + } + + # get matrix that is (p - 1) x J that gives whether or not parameter is zero-comparison + zero_comp <- matrix(FALSE, nrow = n_groups - 1, J) + for (i in 1:(n_groups - 1)) { + zero_comp[i, ] <- (group_counts[1, ] == 0) * (group_counts[i + 1, ] == 0) == 1 + } + + # if there are any zero-comparison parameters + if (sum(zero_comp) > 0) { + cov <- colnames(base_X) + cat <- colnames(Y) + new_zero_comp_dat <- data.frame(covariate = rep(cov, each = length(cat)), + category = rep(cat, length(cov))) + new_zero_comp_dat$zero_comparison <- as.vector(t(zero_comp)) + zero_comp_dat <- rbind(zero_comp_dat, new_zero_comp_dat) + } + } + if (nrow(zero_comp_dat) > 0) { + return(zero_comp_dat) + } + } + } else { + # X is a matrix (not from model matrix), need to check manually + # this will only identify a singular categorical covariate + + # remove intercept + base_X <- X[, -1, drop = FALSE] + + # remove any columns of X that include values other than 0 and 1 + non_cat_cols <- which(apply(base_X, 2, function(x) {sum(!(x %in% c(0, 1)))} > 0)) + if (length(non_cat_cols) > 0 & length(non_cat_cols) == ncol(base_X)) { + return(NULL) + } else { + if (length(non_cat_cols) > 0) { + base_X <- base_X[, -non_cat_cols, drop = FALSE] + } + } + + # check if there are more than 1 columns left + if (ncol(base_X) > 1) { + + # check if there is a single categorical covariate + # i.e. all rows add up to 0 or 1 + if (all(rowSums(base_X) %in% c(0, 1))) { + + # get indices for each group + n_groups <- ncol(base_X) + 1 + group_ind <- vector(mode = "list", length = n_groups) + group_ind[[1]] <- which(rowSums(base_X) == 0) + + for (i in 1:ncol(base_X)) { + group_ind[[i + 1]] <- which(base_X[, i] == 1) + } + + # get Y sums for each group + J <- ncol(Y) + group_counts <- matrix(NA, nrow = n_groups, ncol = J) + for (i in 1:n_groups) { + group_counts[i, ] <- colSums(Y[group_ind[[i]], ]) + } + + # get matrix that is (p - 1) x J that gives whether or not parameter is zero-comparison + zero_comp <- matrix(NA, nrow = n_groups - 1, J) + for (i in 1:(n_groups - 1)) { + zero_comp[i, ] <- (group_counts[1, ] == 0) * (group_counts[i + 1, ] == 0) == 1 + } + + # if there are any zero-comparison parameters + if (sum(zero_comp) > 0) { + cov <- colnames(base_X) + cat <- colnames(Y) + zero_comp_dat <- data.frame(covariate = rep(cov, each = length(cat)), + category = rep(cat, length(cov))) + zero_comp_dat$zero_comparison <- as.vector(t(zero_comp)) + return(zero_comp_dat) + } + + } + + } + } + + # if there are no zero-comparison parameters, return NULL + return(NULL) + +} \ No newline at end of file diff --git a/README.md b/README.md index ecdccf4..6fc45f8 100644 --- a/README.md +++ b/README.md @@ -1,5 +1,4 @@ -# radEmu - +# radEmu [![R-CMD-check](https://github.com/statdivlab/radEmu/workflows/R-CMD-check/badge.svg)](https://github.com/statdivlab/radEmu/actions) @@ -48,7 +47,7 @@ We are currently only releasing `radEmu` via GitHub. If you'd like us to conside ## Use -The vignette demonstrates example usage of the main functions. Please [file an issue](https://github.com/statdivlab/radEmu/issues) if you have a request for a tutorial that is not currently included. The following code shows the easy-to-use syntax if your data is in a `phyloseq` object: +The vignettes demonstrate example usage of the main functions. Please [file an issue](https://github.com/statdivlab/radEmu/issues) if you have a request for a tutorial that is not currently included. The following code shows the easy-to-use syntax if your data is in a `phyloseq` object: ``` r ch_fit <- emuFit(formula = ~ Group + Study + Gender + Sampling, @@ -62,7 +61,9 @@ all_fit <- emuFit(formula = ~ Group + Study + Gender + Sampling, data = my_covariates_df, Y = my_abundances_df) ``` +## Documentation +We additionally have a `pkgdown` [website](https://statdivlab.github.io/radEmu/) that contains pre-built versions of our function [documentation](https://statdivlab.github.io/radEmu/reference/index.html) and our vignettes (an introductory [vignette](https://statdivlab.github.io/radEmu/articles/intro_radEmu.html), an introductory [vignette](https://statdivlab.github.io/radEmu/articles/intro_radEmu_with_phyloseq.html) that uses `phyloseq` data, a [vignette](https://statdivlab.github.io/radEmu/articles/parallel_radEmu.html) for running `radEmu` tests in parallel for more efficient computation, and a [vignette](https://statdivlab.github.io/radEmu/articles/radEmu_clustered_data.html) for running `radEmu` with clustered data). ## Citation diff --git a/_pkgdown.yml b/_pkgdown.yml new file mode 100644 index 0000000..ac7975b --- /dev/null +++ b/_pkgdown.yml @@ -0,0 +1,4 @@ +url: https://statdivlab.github.io/radEmu/ +template: + bootstrap: 5 + diff --git a/man/emuFit.Rd b/man/emuFit.Rd index 8cd0d04..5a0222b 100644 --- a/man/emuFit.Rd +++ b/man/emuFit.Rd @@ -9,9 +9,11 @@ emuFit( X = NULL, formula = NULL, data = NULL, + assay_name = NULL, cluster = NULL, penalize = TRUE, B = NULL, + B_null_list = NULL, fitted_model = NULL, refit = TRUE, test_kj = NULL, @@ -22,6 +24,7 @@ emuFit( use_fullmodel_info = FALSE, use_fullmodel_cov = FALSE, use_both_cov = FALSE, + match_row_names = TRUE, constraint_fn = pseudohuber_center, constraint_grad_fn = dpseudohuber_center_dx, constraint_param = 0.1, @@ -39,7 +42,11 @@ emuFit( inner_maxit = 25, max_step = 1, trackB = FALSE, - return_both_score_pvals = FALSE + return_nullB = FALSE, + return_score_components = FALSE, + return_both_score_pvals = FALSE, + remove_zero_comparison_pvals = 0.01, + unobserved_taxon_error = TRUE ) } \arguments{ @@ -51,15 +58,22 @@ emuFit( \item{data}{an n x p data frame containing variables given in \code{formula}} -\item{cluster}{a numeric vector giving cluster membership for each row of Y to -be used in computing GEE test statistics. Default is NULL, in which case rows of -Y are treated as independent.} +\item{assay_name}{a string containing the desired assay name within a \code{TreeSummarizedExperiment} object. +This is only required if Y is a \code{TreeSummarizedExperiment} object, otherwise this argument does nothing +and can be ignored.} + +\item{cluster}{a vector giving cluster membership for each row of Y to be used in computing +GEE test statistics. Default is NULL, in which case rows of Y are treated as independent.} \item{penalize}{logical: should Firth penalty be used in fitting model? Default is TRUE.} \item{B}{starting value of coefficient matrix (p x J). If not provided, B will be initiated as a zero matrix.} +\item{B_null_list}{list of starting values of coefficient matrix (p x J) for null estimation. This should either +be a list with the same length as \code{test_kj}. If you only want to provide starting values for some tests, +include the other elements of the list as \code{NULL}.} + \item{fitted_model}{a fitted model produced by a separate call to emuFit; to be provided if score tests are to be run without refitting the full unrestricted model. Default is NULL.} @@ -91,9 +105,13 @@ recomputed for each null model fit for score testing.} empirical score covariance evaluated both under the null and full models? Used in simulations} +\item{match_row_names}{logical: Make sure rows on covariate data and response data correspond to +the same sample by comparing row names and subsetting/reordering if necessary.} + \item{constraint_fn}{function g defining a constraint on rows of B; g(B_k) = 0 for rows k = 1, ..., p of B. Default function is a smoothed median (minimizer of -pseudohuber loss).} +pseudohuber loss). If a number is provided a single category constraint will be used +with the provided category as a reference category.} \item{constraint_grad_fn}{derivative of constraint_fn with respect to its arguments (i.e., elements of a row of B)} @@ -103,7 +121,10 @@ parameter controlling relative weighting of elements closer and further from cen (Limit as \code{constraint_param} approaches infinity is the mean; as this parameter approaches zero, the minimizer of the pseudo-Huber loss approaches the median.)} -\item{verbose}{provide updates as model is being fitted? Defaults to TRUE.} +\item{verbose}{provide updates as model is being fitted? Defaults to FALSE. If user sets verbose = TRUE, +then key messages about algorithm progress will be displayed. If user sets verbose = "development", +then key messages and technical messages about convergence will be displayed. Most users who want status +updates should set verbose = TRUE.} \item{tolerance}{tolerance for stopping criterion in full model fitting; once no element of B is updated by more than this value in a single step, we exit @@ -148,21 +169,47 @@ will be rescaled if a step in any parameter exceeds this value. Defaults to 0.5. \item{trackB}{logical: should values of B be recorded across optimization iterations and be returned? Primarily used for debugging. Default is FALSE.} +\item{return_nullB}{logical: should values of B under null hypothesis be returned. Primarily used for debugging. Default is FALSE.} + +\item{return_score_components}{logical: should components of score statistic be returned? Primarily used for debugging. Default is FALSE.} + \item{return_both_score_pvals}{logical: should score p-values be returned using both information matrix computed from full model fit and from null model fits? Default is FALSE. This parameter is used for simulations - in any applied analysis, type of p-value to be used should be chosen before conducting tests.} + +\item{remove_zero_comparison_pvals}{Should score p-values be replaced with NA for zero-comparison parameters? These parameters occur +for categorical covariates with three or more levels, and represent parameters that compare a covariate level to the reference level for +a category in which the comparison level and reference level both have 0 counts in all samples. These parameters can have misleadingly +small p-values and are not thought to have scientifically interesting signals. We recommend removing them before analyzing data further. +If TRUE, all zero-comparison parameter p-values will be set to NA. If FALSE no zero-comparison parameter p-values will be set to NA. +If a value between 0 and 1, all zero-comparison p-values below the value will be set to NA. +Default is \code{0.01}.} + +\item{unobserved_taxon_error}{logical: should an error be thrown if Y includes taxa that have 0 counts for all samples? Default is TRUE.} } \value{ A list containing elements 'coef', 'B', 'penalized', 'Y_augmented', -'I', and 'Dy'. Parameter estimates by -covariate and outcome category (e.g., taxon for microbiome data), as well as -optionally confidence intervals and p-values, are contained in 'coef'. 'B' -contains parameter estimates in matrix format (rows indexing covariates and -columns indexing outcome category / taxon). 'penalized' is equal to TRUE -if Firth penalty is used in estimation (default) and FALSE otherwise. 'I' and -'Dy' contain an information matrix and empirical score covariance matrix -computed under the full model. +'z_hat', 'I', 'Dy', and 'score_test_hyperparams' if score tests are run. +Parameter estimates by covariate and outcome category (e.g., taxon for microbiome data), +as well as optionally confidence intervals and p-values, are contained in 'coef'. +Any robust score statistics and score test p-values are also included in 'coef'. +If there are any zero-comparison parameters in the model, a column 'zero_comparison' +is also included, which is TRUE for any parameters that compare the level of a categorical +covariate to a reference level for a category with only zero counts for both the comparison +level and the reference level. This check is currently implemented for an arbitrary design matrix +generated using the \code{formula} and \code{data} arguments, and for a design matrix with no more +than one categorical covariate if the design matrix \code{X} is input directly. +'B' contains parameter estimates in matrix format (rows indexing covariates and +columns indexing outcome category / taxon). +'penalized' is equal to TRUE f Firth penalty is used in estimation (default) and FALSE otherwise. +'z_hat' returns the nuisance parameters calculated in Equation 7 of the radEmu manuscript, +corresponding to either 'Y_augmented' or 'Y' if the 'penalized' is equal to TRUE +or FALSE, respectively. +I' and 'Dy' contain an information matrix and empirical score covariance matrix computed under the +full model. +'score_test_hyperparams' contains parameters and hyperparameters related to estimation under the null, +including whether or not the algorithm converged, which can be helpful for debugging. } \description{ Fit radEmu model diff --git a/docs/radEmu_hex.png b/man/figures/logo.png similarity index 100% rename from docs/radEmu_hex.png rename to man/figures/logo.png diff --git a/man/plot.emuFit.Rd b/man/plot.emuFit.Rd new file mode 100644 index 0000000..673e857 --- /dev/null +++ b/man/plot.emuFit.Rd @@ -0,0 +1,76 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/plot_emuFit.R +\name{plot.emuFit} +\alias{plot.emuFit} +\title{Plotting function} +\usage{ +\method{plot}{emuFit}( + x, + plot_key = NULL, + title = NULL, + taxon_names = NULL, + display_taxon_names = TRUE, + data_only = FALSE, + ... +) +} +\arguments{ +\item{x}{Output from emuFit()} + +\item{plot_key}{(Optional) Default \code{NULL}. List of named vectors containing names in the "covariate" column of the \code{coef} output of the radEmu model object. If you wish for multiple covariate values to be plotted on the same plot, then those variables should be included in the same named vector. By default, each column of the design matrix receives its own plot.} + +\item{title}{(Optional). Default \code{NULL}. Character string. The main title for the graphic.} + +\item{taxon_names}{(Optional). Default \code{NULL}. Data frame. If \code{NULL}, keep taxon names as listed in radEmu model. Otherwise, users can input a data frame with two columns: one labelled "category" with the same levels as in the radEmu output and another labelled "cat_small" with the preferred labels.} + +\item{display_taxon_names}{(Optional). Default \code{TRUE}. Boolean. If \code{FALSE}, remove sample names from the plot.} + +\item{data_only}{(Optional). Default \code{FALSE}. Boolean. If \code{TRUE}, only returns data frame.} + +\item{...}{There are no optional parameters at this time.} +} +\value{ +Object of class \code{ggplot}. Plot of \code{radEmu} model fit with 95\% confidence intervals. +} +\description{ +Plotting function +} +\examples{ +data(wirbel_sample) +data(wirbel_otu) + +subset_studies <- which(wirbel_sample$Study \%in\% c("FR-CRC", "US-CRC", "AT-CRC")) + +chosen_genera <- c("Eubacterium", "Faecalibacterium", "Fusobacterium", "Porphyromonas") + +mOTU_names <- colnames(wirbel_otu) +mOTU_name_df <- data.frame(name = mOTU_names) \%>\% + dplyr::mutate(base_name = stringr::str_remove(mOTU_names, "unknown ") \%>\% + stringr::str_remove("uncultured ")) \%>\% + dplyr::mutate(genus_name = stringr::word(base_name, 1)) + +restricted_mOTU_names <- mOTU_name_df \%>\% + dplyr::filter(genus_name \%in\% chosen_genera) \%>\% + dplyr::pull(name) + +small_Y <- wirbel_otu[subset_studies, restricted_mOTU_names] +category_to_rm <- which(colSums(small_Y) == 0) + +small_sample <- wirbel_sample[subset_studies, ] + +ch_fit <- emuFit(formula = ~ Group + Study, + data = small_sample, + Y = small_Y, + run_score_tests = FALSE) + +plot_key <- list(p1 = c("Control" = "GroupCTR"), + p2 = c("FR-Control" = "StudyFR-CRC", + "US-Control" = "StudyUS-CRC")) + +out <- plot(x = ch_fit, + plot_key = plot_key, + display_taxon_names = FALSE) + +out$plots$p1 +out$plots$p2 +} diff --git a/man/score_test.Rd b/man/score_test.Rd index 91921bd..5f9af0f 100644 --- a/man/score_test.Rd +++ b/man/score_test.Rd @@ -33,7 +33,8 @@ score_test( ) } \arguments{ -\item{B}{value of coefficient matrix (p x J) returned by full model fit} +\item{B}{value of coefficient matrix (p x J) returned by full model fit or value of coefficient +matrix to start null estimation at given as input to emuFit} \item{Y}{an n x J matrix or dataframe of \emph{augmented} nonnegative observations (i.e., observations Y plus augmentations from last iteration of maximum penalized likelihood estimation @@ -115,13 +116,16 @@ Y are treated as independent.} } \value{ A list containing elements 'score_stat', 'pval', 'log_pval','niter', -'convergence', 'gap', 'u', 'rho', 'null_B', and 'Bs'. 'score_stat' gives the +'convergence', 'gap', 'u', 'rho', 'tau', 'inner_maxit', 'null_B', and 'Bs'. 'score_stat' gives the value of the robust score statistic for H_0: B_{k_constr,j_constr} = g(B_{k_constr}). 'pval' and 'log_pval' are the p-value (on natural and log scales) corresponding to the score statistic (log_pval may be useful when the p-value is very close to zero). 'gap' is the final value of g(B_{k_constr}) - B_{k_constr, j_constr} obtained in optimization under the null. 'u' and 'rho' are final values of augmented -Lagrangian parameters returned by null fitting algorithm. 'null_B' is the value of +Lagrangian parameters returned by null fitting algorithm. 'tau' is the final value of 'tau' that +is used to update the 'rho' values and 'inner_maxit' is the final maximum number of iterations for +the inner optimization loop in optimization under the null, in which B and z parameter values are +maximized for specific 'u' and 'rho' parameters. 'null_B' is the value of B returned but the null fitting algorithm. 'Bs' is by default NULL; if trackB = TRUE, 'Bs is a data frame containing values of B by outcome category, covariate, and iteration. diff --git a/man/simulate_data.Rd b/man/simulate_data.Rd new file mode 100644 index 0000000..9f7ee00 --- /dev/null +++ b/man/simulate_data.Rd @@ -0,0 +1,49 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/simulate_data.R +\name{simulate_data} +\alias{simulate_data} +\title{Data simulation function} +\usage{ +simulate_data( + n, + J, + b0 = NULL, + b1 = NULL, + distn, + zinb_size = NULL, + zinb_zero_prop = NULL, + mean_z, + X = NULL, + B = NULL, + cluster = NULL +) +} +\arguments{ +\item{n}{Number of samples} + +\item{J}{Number of categories} + +\item{b0}{Intercept parameter vector} + +\item{b1}{Covariate paramter vector} + +\item{distn}{Distribution to simulate from, either "Poisson" or "ZINB"} + +\item{zinb_size}{Size parameter for negative binomial draw for ZINB data} + +\item{zinb_zero_prop}{Proportion of zeros for ZINB data} + +\item{mean_z}{Parameter controlling the mean of the sample-specific effects.} + +\item{X}{Optional design matrix, this must have two columns and n rows.} + +\item{B}{Optional B matrix, if p is not equal to 2} + +\item{cluster}{Optional cluster vector, this must have n elements.} +} +\value{ +\code{Y}. A \code{n times J} dimension matrix of simulated response counts. +} +\description{ +Function to simulate data for simulations in Clausen & Willis (2024) and for the cluster vignette +} diff --git a/pkgdown/favicon/apple-touch-icon-120x120.png b/pkgdown/favicon/apple-touch-icon-120x120.png new file mode 100644 index 0000000..d41e6b1 Binary files /dev/null and b/pkgdown/favicon/apple-touch-icon-120x120.png differ diff --git a/pkgdown/favicon/apple-touch-icon-152x152.png b/pkgdown/favicon/apple-touch-icon-152x152.png new file mode 100644 index 0000000..d8f7096 Binary files /dev/null and b/pkgdown/favicon/apple-touch-icon-152x152.png differ diff --git a/pkgdown/favicon/apple-touch-icon-180x180.png b/pkgdown/favicon/apple-touch-icon-180x180.png new file mode 100644 index 0000000..c5c9256 Binary files /dev/null and b/pkgdown/favicon/apple-touch-icon-180x180.png differ diff --git a/pkgdown/favicon/apple-touch-icon-60x60.png b/pkgdown/favicon/apple-touch-icon-60x60.png new file mode 100644 index 0000000..0c95625 Binary files /dev/null and b/pkgdown/favicon/apple-touch-icon-60x60.png differ diff --git a/pkgdown/favicon/apple-touch-icon-76x76.png b/pkgdown/favicon/apple-touch-icon-76x76.png new file mode 100644 index 0000000..7369776 Binary files /dev/null and b/pkgdown/favicon/apple-touch-icon-76x76.png differ diff --git a/pkgdown/favicon/apple-touch-icon.png b/pkgdown/favicon/apple-touch-icon.png new file mode 100644 index 0000000..b2bade4 Binary files /dev/null and b/pkgdown/favicon/apple-touch-icon.png differ diff --git a/pkgdown/favicon/favicon-16x16.png b/pkgdown/favicon/favicon-16x16.png new file mode 100644 index 0000000..4ee3184 Binary files /dev/null and b/pkgdown/favicon/favicon-16x16.png differ diff --git a/pkgdown/favicon/favicon-32x32.png b/pkgdown/favicon/favicon-32x32.png new file mode 100644 index 0000000..63e4086 Binary files /dev/null and b/pkgdown/favicon/favicon-32x32.png differ diff --git a/pkgdown/favicon/favicon.ico b/pkgdown/favicon/favicon.ico new file mode 100644 index 0000000..b5f2da9 Binary files /dev/null and b/pkgdown/favicon/favicon.ico differ diff --git a/tests/testthat/test-augmentation-failures.R b/tests/testthat/test-augmentation-failures.R index 960da78..43df471 100644 --- a/tests/testthat/test-augmentation-failures.R +++ b/tests/testthat/test-augmentation-failures.R @@ -1,11 +1,17 @@ test_that("confirm Matrix Csparse_transpose issue is not happening", { - Y <- structure(c(1087, 3541, 0, 2432, 0, 18538, 1158, 2282, 625, 0, - 3759, 0, 0, 4658, 0, 3719, 0, 0, 7531, 48316, 0, 0, 20273, 0, - 1227, 0, 1471, 0, 479, 10602, 3115, 3286, 0, 1969, 0, 3045, 0, - 8018, 0, 1622, 1307, 34117, 9338, 0, 0, 7909, 0, 0), dim = c(12L, 4L)) + set.seed(1) X <- structure(c(rep(1, 18), rep(0, 6)), dim = c(12L, 2L)) - covariates <- structure(list(group = c(0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1)), class = "data.frame", row.names = c(NA, -12L)) + + Y <- radEmu:::simulate_data(n = 12L, J = 4L, + b0 = runif(4L, min = 0, max = 4), + b1 = runif(4L, min = 0, max = 4), + X = X, + distn = "Poisson", + mean_z = 5) + + covariates <- structure(list(group = c(0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1)), + class = "data.frame", row.names = c(NA, -12L)) # devtools::load_all() # info <- methods::as(info, "symmetricMatrix") fitted_model <- emuFit(Y = Y, @@ -24,16 +30,16 @@ test_that("confirm Matrix Csparse_transpose issue is not happening", { ### check data frame inputs ok fitted_model_df <- emuFit(Y = as.data.frame(Y)[5:8, ], - X = as.data.frame(X)[5:8, ], - formula = ~group, - data = covariates, - verbose = FALSE, - B_null_tol = 1e-2, - tolerance = 0.01, - tau = 2, - run_score_test = TRUE, - return_wald_p = TRUE) - expect_true("emuFit" %in% class(fitted_model_df)) + X = as.data.frame(X)[5:8, ], + formula = ~group, + data = covariates, + verbose = FALSE, + B_null_tol = 1e-2, + tolerance = 0.01, + tau = 2, + run_score_test = TRUE, + return_wald_p = TRUE) + expect_true("emuFit" %in% class(fitted_model_df)) }) diff --git a/tests/testthat/test-cluster.R b/tests/testthat/test-cluster.R index e01a88b..1393c3c 100644 --- a/tests/testthat/test-cluster.R +++ b/tests/testthat/test-cluster.R @@ -1,28 +1,66 @@ -set.seed(11) -J <- 6 -p <- 2 -n <- 12 -X <- cbind(1,rnorm(n)) -z <- rnorm(n) +5 -b0 <- rnorm(J) -b1 <- seq(1,5,length.out = J) -b1 <- b1 - mean(b1) -b <- rbind(b0,b1) -Y <- matrix(NA,ncol = J, nrow = n) +test_that("clusters work as I want", { + + set.seed(100) + n <- 64 + J <- 50 + + cage_num <- rep(c(1:16), 4) + treatment <- (cage_num <= 8) + XX <- data.frame(treatment) + + Y <- radEmu:::simulate_data(n = n, + J = J, + X = cbind(1, treatment), + b0 = runif(J, min = 0, max = 4), + b1 = runif(J, min = 0, max = 4), + distn = "ZINB", + zinb_size = 10, + zinb_zero_prop = 0.3, + mean_z = 5) + + # check that cluster argument works as a numeric vector + ef_num <- emuFit(formula = ~ treatment, + data = XX, + Y = Y, + cluster=cage_num, + run_score_tests=FALSE) #### very fast + expect_equal(ef_num$coef %>% class, "data.frame") + + # check that cluster argument works as character vector and gives + # equivalent results to numeric vector + cage_char <- rep(c(LETTERS[1:16]), 4) + ef_char <- emuFit(formula = ~ treatment, + data = XX, + Y = Y, + cluster=cage_char, + run_score_tests=FALSE) + expect_equal(ef_num$coef, ef_char$coef) + + # check that cluster argument works as factor and gives equivalent results + # to numeric vector + cage_fact <- cage_char %>% as.factor + ef_fact <- emuFit(formula = ~ treatment, + data = XX, + Y = Y, + cluster=cage_fact, + run_score_tests=FALSE) + expect_equal(ef_num$coef, ef_fact$coef) +}) -for(i in 1:n){ - for(j in 1:J){ - temp_mean <- exp(X[i,,drop = FALSE]%*%b[,j,drop = FALSE] + z[i]) - Y[i,j] <- rnbinom(1, mu= temp_mean,size = 2)*rbinom(1,1,0.8) - } -} + +set.seed(11) +X <- cbind(1,rnorm(12)) +Y <- radEmu:::simulate_data(n = 12, + J = 6, + X = X, + b0 = rnorm(6), + b1 = seq(1,5,length.out = 6) - mean(seq(1,5,length.out = 6)), + distn = "ZINB", + zinb_size = 2, + zinb_zero_prop = 0.8, + mean_z = 5) covariates <- data.frame(group = X[,2]) -b0 <- rnorm(J) -b1 <- seq(1,5,length.out = J) -b1 <- b1 - mean(b1) -b1[3:4] <- 0 -b <- rbind(b0,b1) test_that("GEE with cluster covariance gives plausible type 1 error ",{ @@ -39,24 +77,21 @@ test_that("GEE with cluster covariance gives plausible type 1 error ",{ pval = numeric(2*nsim))[-(1:(2*nsim)),] results_noGEE <- results for(sim in 1:nsim){ - print(sim) - X <- cbind(1,rnorm(n)) - covariates <- data.frame(group = X[,2]) - Y <- matrix(NA,ncol = J, nrow = n) - - cluster_effs <- lapply(1:4, - function(i) - log(matrix(rexp(2*J),nrow= 2))) + # print(sim) - for(i in 1:n){ - Y[i,] <- 0 - while(sum(Y[i,])==0){ - for(j in 1:J){ - temp_mean <- exp(X[i,,drop = FALSE]%*%(b[,j,drop = FALSE] + - cluster_effs[[ cluster[i] ]][,j]) + z[i]) - Y[i,j] <- rnbinom(1, mu= temp_mean,size = 5)*rbinom(1,1,0.8) - }} - } + X <- cbind(1, rnorm(12)) + Y <- radEmu:::simulate_data(n = 12, + J = 6, + X = X, + b0 = rnorm(6), + b1 = seq(1,5,length.out = 6) - + mean(seq(1,5,length.out = 6)), + distn = "ZINB", + zinb_size = 5, + zinb_zero_prop = 0.8, + mean_z = 5, + cluster = cluster) + covariates <- data.frame(group = X[,2]) # expect_silent({ fitted_model_cluster <- emuFit(Y = Y, diff --git a/tests/testthat/test-emuFit.R b/tests/testthat/test-emuFit.R index 9cb0024..ae7f9d8 100644 --- a/tests/testthat/test-emuFit.R +++ b/tests/testthat/test-emuFit.R @@ -1,60 +1,26 @@ set.seed(11) J <- 6 -p <- 2 n <- 12 X <- cbind(1,rnorm(n)) -z <- rnorm(n) +5 b0 <- rnorm(J) -b1 <- seq(1,5,length.out = J) -b1 <- b1 - mean(b1) -b <- rbind(b0,b1) -Y <- matrix(NA,ncol = J, nrow = n) +b1 <- seq(1,5,length.out = J) - + mean(seq(1,5,length.out = J)) +b <- rbind(b0, b1) +Y <- radEmu:::simulate_data(n = n, + J = J, + X = X, + b0 = b0, + b1 = b1, + distn = "ZINB", + zinb_size = 2, + zinb_zero_prop = 0.2, + mean_z = 5) -for(i in 1:n){ - for(j in 1:J){ - temp_mean <- exp(X[i,,drop = FALSE]%*%b[,j,drop = FALSE] + z[i]) - Y[i,j] <- rnbinom(1, mu= temp_mean,size = 2)*rbinom(1,1,0.8) - } -} - -# Y <- structure(c(534337, 0, 0, 0, 376, 41, 19, 103, 0, 0, 85, 0, 42794, -# 0, 0, 0, 95, 0, 0, 15, 0, 0, 0, 26, 0, 149, 0, 0, 0, 0, 0, 211, -# 0, 0, 0, 0, 0, 103, 0, 0, 0, 1372, 83, 337, 0, 0, 0, 0, 0, 53, -# 0, 0, 0, 0, 259, 0, 0, 0, 14, 0, 0, 0, 0, 193, 0, 0, 0, 0, 0, -# 0, 402, 0), dim = c(12L, 6L)) -# X <- structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, -2.49421928123597, -# -0.579053917775617, -0.974555155010523, 0.237710056670222, 0.41240637454179, -# 1.12994912631468, 0.706485861932659, 0.588878125500377, 0.0834756145662259, -# 1.99483775157368, 0.227951737778031, -1.03963299361785), dim = c(12L, -# 2L)) +#To ensure the messages about lack of row names do not show in the tests +rownames(X) <- paste0("Sample_",1:12) +rownames(Y) <- paste0("Sample_",1:12) -# -# Y <- structure(c(1748, 8286, 4096, 4289, 1122, 30007, 5087, 3841, -# 3059, 3105, 80, 32, 0, 20, 13, 0, 0, 30, 41, 54, 0, 124134, 43569, -# 122134, 15785, 99540, 0, 41104, 0, 0, 0, 0, 0, 572, 0, 1497, -# 0, 0, 314, 0, 0, 0, 0, 416, 920, 0, 0, 1931, 1279, 0, 0, 0, 2, -# 21, 49, 41, 0, 89, 0, 85, 1287, 1716, 0, 0, 1354, 8783, 3040, -# 6271, 2274, 0, 26431, 5186, 4147, 0, 0, 6450, 0, 0, 1483, 0, -# 0, 0, 0, 5936, 0, 0, 0, 33557, 11459, 0, 0, 4065, 0, 5391, 6721, -# 8997, 9225, 13951, 4061, 3871, 0, 0, 0, 0, 0, 0, 3954, 1338, -# 886, 426, 0, 0, 0, 496, 0, 709, 508, 840, 680, 0, 0, 4529, 2885, -# 0, 0, 13382, 11802, 0, 2144, 2622, 92214, 18326, 6183, 11737, -# 0, 0, 12808, 7604, 4684, 9348, 0, 3564, 2250, 0, 0, 20486, 0, -# 3442, 5133, 5103, 299927, 34273, 17407, 0, 17896, 149402, 42592, -# 0, 0, 25395, 1875, 21975, 1685, 0, 1654, 28670, 9331, 0, 4765, -# 0, 0, 0, 76158, 0, 85495, 247396, 51659, 93587, 0, 84277, 0, -# 0, 217, 0, 0, 0, 0, 0, 0, 0, 0, 76950, 0, 19911, 33042, 43690, -# 68014, 43885, 6086, 0), dim = c(20L, 10L)) -# X <- structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, -# 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, -# 1, 1), dim = c(20L, 2L)) covariates <- data.frame(group = X[,2]) -# b <- structure(c(0.0231122313161752, -4.5, 2.34881909126062, -3.5, -# -0.949623561962149, -2.5, -0.23942645176718, 0.556978794649417, -# 0.681191947270542, 0.499524350059682, -1.14509952873781, 0.5, -# 0.258373890958553, 1.5, 0.163812326430595, 2.5, 0.491510413255902, -# 3.5, -1.63267035802524, 4.5), dim = c(2L, 10L), dimnames = list( -# c("b0", "b1"), NULL)) test_that("emuFit takes formulas and actually fits a model", { @@ -219,12 +185,6 @@ test_that("emuFit takes cluster argument without breaking ",{ }) -b0 <- rnorm(J) -b1 <- seq(1,5,length.out = J) -b1 <- b1 - mean(b1) -b1[3:4] <- 0 -b <- rbind(b0,b1) - test_that("GEE with cluster covariance gives plausible type 1 error ",{ skip("Skipping -- test requires fitting models to 100 simulated datasets.") set.seed(44022) @@ -238,24 +198,24 @@ test_that("GEE with cluster covariance gives plausible type 1 error ",{ results_noGEE <- results for(sim in 1:nsim){ print(sim) - X <- cbind(1,rnorm(n)) + + X <- cbind(1,rnorm(12)) covariates <- data.frame(group = X[,2]) - Y <- matrix(NA,ncol = J, nrow = n) - - cluster_effs <- lapply(1:4, - function(i) - log(matrix(rexp(2*J),nrow= 2))) - - for(i in 1:n){ - Y[i,] <- 0 - while(sum(Y[i,])==0){ - for(j in 1:J){ - temp_mean <- exp(X[i,,drop = FALSE]%*%(b[,j,drop = FALSE] + - cluster_effs[[ cluster[i] ]][,j]) + z[i]) - Y[i,j] <- rnbinom(1, mu= temp_mean,size = 5)*rbinom(1,1,0.8) - }} - } - + + b1 <- seq(1,5,length.out = 6) + b1 <- b1 - mean(b1) + b1[3:4] <- 0 + + Y <- radEmu:::simulate_data(n = 12, J = 6, + X = X, + b0 = rnorm(6), + b1 = b1, + distn = "ZINB", + zinb_size = 2, + zinb_zero_prop = 0.2, + mean_z = 5, + cluster = cluster) + # expect_silent({ fitted_model_cluster <- emuFit(Y = Y, X = X, @@ -327,7 +287,7 @@ test_that("GEE with cluster covariance gives plausible type 1 error ",{ # b1 <- b1 - mean(b1) # b1[5] <- pseudohuber_center(b1[-5],0.1) # -# Y <- simulate_data(n=10, J=10, b0=rnorm(10), distn="Poisson", b1=b1, mean_count_before_ZI=500) +# Y <- radEmu:::simulate_data(n=10, J=10, b0=rnorm(10), distn="Poisson", b1=b1, mean_z=500) # set.seed(894334) # n <- 100 @@ -386,3 +346,257 @@ test_that("GEE with cluster covariance gives plausible type 1 error ",{ # expect_true(cor(fitted_model$B[2,],b1)>0.95) # # }) + +test_that("emuFit runs without penalty", { + + expect_silent({ + fitted_model <- emuFit(Y = Y, + X = X, + penalize = FALSE, + formula = ~group, + data = covariates, + verbose = FALSE, + B_null_tol = 1e-2, + tolerance = 0.01, + tau = 2, + return_wald_p = FALSE, + compute_cis = TRUE, + run_score_tests = TRUE, + use_fullmodel_info = FALSE, + use_fullmodel_cov = FALSE, + return_both_score_pvals = FALSE) + }) +}) + +test_that("emuFit runs with just intercept model", { + + expect_message({ + fitted_model <- emuFit(Y = Y, + formula = ~1, + data = covariates, + verbose = FALSE, + B_null_tol = 1e-2, + tolerance = 0.01, + tau = 2, + return_wald_p = FALSE, + compute_cis = TRUE, + run_score_tests = TRUE, + use_fullmodel_info = FALSE, + use_fullmodel_cov = FALSE, + return_both_score_pvals = FALSE) + }) + + expect_message({ + fitted_model1 <- emuFit(Y = Y, + X = X[, 1, drop = FALSE], + verbose = FALSE, + B_null_tol = 1e-2, + tolerance = 0.01, + tau = 2, + return_wald_p = FALSE, + compute_cis = TRUE, + run_score_tests = TRUE, + use_fullmodel_info = FALSE, + use_fullmodel_cov = FALSE, + return_both_score_pvals = FALSE) + }) + + expect_equal(fitted_model$coef[, 2:9], fitted_model1$coef[, 2:9]) + +}) + +test_that("emuFit has 'score_test_hyperparams' object and throws warnings when convergence isn't hit", { + # check that warning is returned when estimation under the alternative doesn't converge + expect_warning({ + fitted_model <- emuFit(Y = Y, + X = cbind(X, rnorm(nrow(X))), + verbose = FALSE, + B_null_tol = 1e-2, + tolerance = 0.01, + tau = 2, + return_wald_p = FALSE, + compute_cis = FALSE, + run_score_tests = FALSE, + maxit = 1) + }) + expect_false(fitted_model$estimation_converged) + + # check that warning is returned when estimation under the null doesn't converge + suppressWarnings({ + fitted_model <- emuFit(Y = Y, + X = cbind(X, rnorm(nrow(X))), + verbose = FALSE, + B_null_tol = 1e-2, + tolerance = 0.01, + tau = 2, + return_wald_p = FALSE, + compute_cis = FALSE, + run_score_tests = TRUE, + test_kj = data.frame(k = 1, j = 1:2), + maxit = 1, + inner_maxit = 1) + }) + + # check that fitted model contains score_test_hyperparams object + expect_true("score_test_hyperparams" %in% names(fitted_model)) + + # check that fitted model contains data frame of unconverged test_kj + expect_type(fitted_model$null_estimation_unconverged, "list") +}) + +test_that("test that B_null_list object can be used and throws appropriate warnings when used incorrectly", { + expect_warning({ + fitted_model <- emuFit(Y = Y, + X = X, + B_null_list = list(b), + verbose = FALSE, + B_null_tol = 1e-2, + tolerance = 0.01, + tau = 2, + return_wald_p = FALSE, + compute_cis = FALSE, + run_score_tests = TRUE, + test_kj = data.frame(k = 1, j = 1:2)) + }) + + expect_silent({ + fitted_model <- emuFit(Y = Y, + X = X, + B_null_list = list(NULL, b), + verbose = FALSE, + B_null_tol = 1e-2, + tolerance = 0.01, + tau = 2, + return_wald_p = FALSE, + compute_cis = FALSE, + run_score_tests = TRUE, + test_kj = data.frame(k = 1, j = 1:2)) + }) + + expect_warning({ + fitted_model <- emuFit(Y = Y, + X = X, + B_null_list = list(NULL, b[, -3]), + verbose = FALSE, + B_null_tol = 1e-2, + tolerance = 0.01, + tau = 2, + return_wald_p = FALSE, + compute_cis = FALSE, + run_score_tests = TRUE, + test_kj = data.frame(k = 1, j = 1:2)) + }) + +}) + +test_that("emuFit reorders X and X and Y rownames don't match", { + dat1 <- data.frame(group = c(covariates$group[12], covariates$group[1:11])) + rownames(dat1) <- paste0("sample", c(12, 1:11)) + dat2 <- covariates + rownames(dat2) <- paste0("sample", 1:12) + rownames(Y) <- paste0("sample", 1:12) + + expect_message({ + fitted_model1 <- emuFit(formula = ~ group, + data = dat1, + Y = Y, + compute_cis = FALSE, + run_score_tests = FALSE) + }) + + expect_silent({ + fitted_model2 <- emuFit(formula = ~ group, + data = dat2, + Y = Y, + compute_cis = FALSE, + run_score_tests = FALSE) + }) + + expect_true(all.equal(fitted_model1$coef, fitted_model2$coef)) +}) + +test_that("emuFit throws error when there is a category with all zero counts", { + + Y_zero <- Y + Y_zero[, 1] <- 0 + + expect_error({ + fitted_model <- emuFit(Y = Y_zero, + X = X, + formula = ~group, + data = covariates, + verbose = FALSE, + B_null_tol = 1e-2, + tolerance = 0.01, + tau = 2, + return_wald_p = FALSE, + compute_cis = FALSE, + run_score_tests = FALSE) + }) +}) + +test_that("Confirm zi is different when penalty is applied or not", { + + fit_penT <- emuFit(Y = Y, + X = X, + formula = ~group, + data = covariates, + verbose = FALSE, + penalize = TRUE, + B_null_tol = 1e-2, + tolerance = 0.01, + tau = 2, + return_wald_p = FALSE, + compute_cis = FALSE, + run_score_tests = FALSE) + + fit_penF <- emuFit(Y = Y, + X = X, + formula = ~group, + data = covariates, + verbose = FALSE, + penalize = FALSE, + B_null_tol = 1e-2, + tolerance = 0.01, + tau = 2, + return_wald_p = FALSE, + compute_cis = FALSE, + run_score_tests = FALSE) + + expect_false(isTRUE(all.equal(fit_penT$z_hat, + fit_penF$z_hat))) +}) + +test_that("Single category constraint works", { + + emuRes <- emuFit(Y = Y, + X = X, + constraint_fn = 3, + run_score_tests = FALSE) + expect_true(emuRes$B[2, 3] == 0) + +}) + +test_that("emuFit works with fitted objects passed in", { + emuRes <- emuFit(Y = Y, + X = X, + run_score_tests = FALSE) + # can run emuFit with fitted model + expect_silent({ + emuRes2 <- emuFit(Y = Y, X = X, fitted_model = emuRes, refit = FALSE, + compute_cis = FALSE, test_kj = data.frame(k = 2, j = 1)) + }) + # get error if have penalize arguments that don't match + expect_error({ + emuRes2 <- emuFit(Y = Y, X = X, fitted_model = emuRes, refit = FALSE, + compute_cis = FALSE, test_kj = data.frame(k = 2, j = 1), + penalize = FALSE) + }) + # can run emuFit with only B + # can run emuFit with fitted model + expect_silent({ + emuRes2 <- emuFit(Y = Y, X = X, B = emuRes$B, refit = FALSE, + compute_cis = FALSE, test_kj = data.frame(k = 2, j = 1)) + }) + +}) \ No newline at end of file diff --git a/tests/testthat/test-emuFit_TreeSummarizedExperiment.R b/tests/testthat/test-emuFit_TreeSummarizedExperiment.R new file mode 100644 index 0000000..28379db --- /dev/null +++ b/tests/testthat/test-emuFit_TreeSummarizedExperiment.R @@ -0,0 +1,48 @@ +library(radEmu) + +test_that("emuFit works with a TreeSummarizedExperiment object", { + if (requireNamespace("TreeSummarizedExperiment", quietly = TRUE)) { + # make TreeSummarizedExperiment object from data + data(wirbel_sample) + data(wirbel_otu) + data(wirbel_taxonomy) + + sub_samples <- wirbel_sample$Country == "FRA" & wirbel_sample$Gender == "F" + sub_taxa <- colSums(wirbel_otu[sub_samples, , drop = FALSE]) > 0 + + wirbel_sample_sub <- wirbel_sample[sub_samples, ] + wirbel_otu_sub <- wirbel_otu[sub_samples, sub_taxa] + wirbel_taxonomy_sub <- wirbel_taxonomy[sub_taxa, ] + + Y <- TreeSummarizedExperiment::TreeSummarizedExperiment( + assays = list("counts" = t(wirbel_otu_sub)), + rowData = wirbel_taxonomy_sub, + colData = wirbel_sample_sub) + + # fit model using TreeSummarizedExperiment + fit <- emuFit(Y = Y, + formula = ~ Group, + assay_name = "counts", + run_score_tests = FALSE, tolerance = 0.01) + + # fit model using data.frames directly + fit2 <- emuFit(Y = wirbel_otu_sub, + X = model.matrix(object = ~ Group, data = wirbel_sample_sub), + run_score_tests = FALSE, tolerance = 0.01) + + # confirm the results match when data are extracted from TreeSummarizedExperiment + # or when data.frames are used directly + expect_true(all.equal(fit$coef, fit2$coef)) + + # confirm an error is returned when assay_name is not provided by Y is + # a TreeSummarizedExperiment object + expect_error( + fit <- emuFit(Y = Y, + formula = ~ Group, + run_score_tests = FALSE, tolerance = 0.01) + ) + + } else { + expect_error(stop("You don't have TreeSummarizedExperiment installed.")) + } +}) diff --git a/tests/testthat/test-emuFit_micro.R b/tests/testthat/test-emuFit_micro.R index 2865915..9e9de90 100644 --- a/tests/testthat/test-emuFit_micro.R +++ b/tests/testthat/test-emuFit_micro.R @@ -2,58 +2,50 @@ test_that("ML fit to simple example give reasonable output", { set.seed(4323) X <- cbind(1,rep(c(0,1),each = 20)) - z <- rnorm(40) +5 J <- 10 - p <- 2 n <- 40 - b0 <- rnorm(J) - b1 <- seq(1,10,length.out = J) - b <- rbind(b0,b1) - Y <- matrix(NA,ncol = J, nrow = 40) - - for(i in 1:40){ - for(j in 1:J){ - temp_mean <- exp(X[i,,drop = FALSE]%*%b[,j,drop = FALSE] + z[i]) - Y[i,j] <- rpois(1, lambda = temp_mean) - } - } - ml_fit <- emuFit_micro(X, - Y, - constraint_fn = function(x) mean(x), - maxit = 200, - tolerance = 1e-3, - verbose = FALSE) + b1 <- 1:J - mean(1:J) + Y <- radEmu:::simulate_data(n = n, + J = J, + X = X, + b0 = rnorm(J), + b1 = b1, + distn = "Poisson", + mean_z = 10) + + ml_fit <- emuFit_micro(X = X, + Y = Y, + constraint_fn = function(x) mean(x), + maxit = 200, + tolerance = 1e-3, + verbose = FALSE) # plot(b1-mean(b1),ml_fit[2,]) # abline(a = 0,b = 1,lty =2,col = "red") - expect_true(max(abs(ml_fit[2,] - (b1 - mean(b1))))<.1) + expect_true(max(abs(ml_fit[2,] - (b1 - mean(b1)))) < 1) # increased to 1 }) test_that("With or without 'working_constraint' we get same results", { set.seed(4323) X <- cbind(1,rep(c(0,1),each = 20)) - z <- rnorm(40) +5 J <- 10 - p <- 2 n <- 40 - b0 <- rnorm(J) - b1 <- seq(1,10,length.out = J) - b <- rbind(b0,b1) - Y <- matrix(NA,ncol = J, nrow = 40) - - for(i in 1:40){ - for(j in 1:J){ - temp_mean <- exp(X[i,,drop = FALSE]%*%b[,j,drop = FALSE] + z[i]) - Y[i,j] <- rpois(1, lambda = temp_mean) - } - } + Y <- radEmu:::simulate_data(n = n, + J = J, + X = X, + b0 = rnorm(J), + b1 = seq(1,10,length.out = J), + distn = "Poisson", + mean_z = 5) + ml_fit <- emuFit_micro(X, Y, constraint_fn = function(x) mean(x), maxit = 200, tolerance = 1e-6, verbose= FALSE) + ml_fit_direct <- emuFit_micro(X, Y, constraint_fn = function(x) mean(x), @@ -67,42 +59,136 @@ test_that("With or without 'working_constraint' we get same results", { }) + +test_that("PL fit with categorical predictor matches analytical form of MPLE in this case, + and does NOT match MLE when group sizes are equal", { + set.seed(90333) + X <- cbind(1,rep(c(0,1),each = 20)) + z <- rnorm(40) + J <- 10 + p <- 2 + n <- 40 + b0 <- rnorm(J) + b1 <- seq(1,10,length.out = J)/5 + b <- rbind(b0,b1) + Y <- matrix(NA,ncol = J, nrow = 40) + + for(i in 1:40){ + for(j in 1:J){ + temp_mean <- exp(X[i,,drop = FALSE]%*%b[,j,drop = FALSE] + z[i]) + Y[i,j] <- rpois(1, lambda = temp_mean) + } + } + pl_fit <- emuFit_micro_penalized(X, + Y, + B = matrix(rnorm(20),nrow = 2), + constraint_fn = function(x) mean(x), + maxit = 200, + tolerance = 1e-8, + verbose= FALSE) + + ml_fit <- emuFit_micro(X, + Y, + B = matrix(rnorm(20),nrow = 2), + constraint_fn = function(x) mean(x), + maxit = 200, + tolerance = 1e-8, + verbose= FALSE) + + cs_grp1 <- colSums(Y[1:20,]) + cs_grp2 <- colSums(Y[21:40,]) + + bhat1 <- log(cs_grp1 + 0.5) - log(cs_grp1[1] + 0.5) + bhat2 <- log(cs_grp2 + 0.5) - log(cs_grp2[1] +0.5) + bhat2 <- bhat2 - bhat1 + bhat1 <- bhat1 - mean(bhat1) + bhat2 <- bhat2 - mean(bhat2) + analytical_B <- rbind(bhat1,bhat2) + + expect_true(max(abs(pl_fit$B - analytical_B))< 1e-7) + expect_true(max(abs(ml_fit - analytical_B))>0.01) + }) + + +test_that("PL fit with categorical predictor matches analytical form of MPLE in this case, + and does NOT match MLE when group sizes are unequal", { + set.seed(90333) + X <- cbind(1,rep(c(0,0,0,1),each = 10)) + z <- rnorm(40) + J <- 10 + p <- 2 + n <- 40 + b0 <- rnorm(J) + b1 <- seq(1,10,length.out = J)/5 + b <- rbind(b0,b1) + Y <- matrix(NA,ncol = J, nrow = 40) + + for(i in 1:40){ + for(j in 1:J){ + temp_mean <- exp(X[i,,drop = FALSE]%*%b[,j,drop = FALSE] + z[i]) + Y[i,j] <- rpois(1, lambda = temp_mean) + } + } + pl_fit <- emuFit_micro_penalized(X, + Y, + B = matrix(rnorm(20),nrow = 2), + constraint_fn = function(x) mean(x), + maxit = 200, + tolerance = 1e-8, + verbose= FALSE) + + ml_fit <- emuFit_micro(X, + Y, + B = matrix(rnorm(20),nrow = 2), + constraint_fn = function(x) mean(x), + maxit = 200, + tolerance = 1e-8, + verbose= FALSE) + + cs_grp1 <- colSums(Y[1:30,]) + cs_grp2 <- colSums(Y[31:40,]) + + bhat1 <- log(cs_grp1 + 0.5) - log(cs_grp1[1] + 0.5) + bhat2 <- log(cs_grp2 + 0.5) - log(cs_grp2[1] +0.5) + bhat2 <- bhat2 - bhat1 + bhat1 <- bhat1 - mean(bhat1) + bhat2 <- bhat2 - mean(bhat2) + analytical_B <- rbind(bhat1,bhat2) + + expect_true(max(abs(pl_fit$B - analytical_B))< 1e-7) + expect_true(max(abs(ml_fit - analytical_B))>0.01) + }) + test_that("We get same results with and without warm start", { set.seed(4323) X <- cbind(1,rep(c(0,1),each = 20)) - z <- rnorm(40) +5 J <- 10 - p <- 2 n <- 40 - b0 <- rnorm(J) - b1 <- seq(1,10,length.out = J) - b <- rbind(b0,b1) - Y <- matrix(NA,ncol = J, nrow = 40) - - for(i in 1:40){ - for(j in 1:J){ - temp_mean <- exp(X[i,,drop = FALSE]%*%b[,j,drop = FALSE] + z[i]) - Y[i,j] <- rpois(1, lambda = temp_mean) - } - } - - + Y <- radEmu:::simulate_data(n = n, + J = J, + X = X, + b0 = rnorm(J), + b1 = 1:J, + distn = "Poisson", + mean_z = 2) + + # may need to have large number of iterations and small tolerance ml_fit <- emuFit_micro(X, Y, constraint_fn = function(x) mean(x), - maxit = 200, - tolerance = 1e-6, + maxit = 1e3, + tolerance = 1e-14, verbose = FALSE) + ml_fit_direct <- emuFit_micro(X, Y, constraint_fn = function(x) mean(x), - maxit = 200, + maxit = 1e3, warm_start = FALSE, - tolerance = 1e-6, + tolerance = 1e-14, verbose = FALSE) - - expect_equal(ml_fit,ml_fit_direct, tolerance = 1e-6) + expect_equal(ml_fit, ml_fit_direct, tolerance = 1e-6) }) @@ -110,21 +196,16 @@ test_that("We get same results with and without warm start", { test_that("We get a fit if we don't specify constraint", { set.seed(4323) X <- cbind(1,rep(c(0,1),each = 20)) - z <- rnorm(40) +5 J <- 10 - p <- 2 n <- 40 - b0 <- rnorm(J) - b1 <- seq(1,10,length.out = J) - b <- rbind(b0,b1) - Y <- matrix(NA,ncol = J, nrow = 40) - - for(i in 1:40){ - for(j in 1:J){ - temp_mean <- exp(X[i,,drop = FALSE]%*%b[,j,drop = FALSE] + z[i]) - Y[i,j] <- rpois(1, lambda = temp_mean) - } - } + Y <- radEmu:::simulate_data(n = n, + J = J, + X = X, + b0 = rnorm(J), + b1 = 1:J, + distn = "Poisson", + mean_z = 2) + ml_fit <- emuFit_micro(X, Y, maxit = 200, @@ -139,20 +220,18 @@ test_that("ML fit to simple example give reasonable output with J >> n", { set.seed(4323) n <- 10 X <- cbind(1,rep(c(0,1),each = n/2)) - z <- rnorm(n) +8 J <- 1000 b0 <- rnorm(J) b1 <- seq(-5,5,length.out = J) b <- rbind(b0,b1) - Y <- matrix(NA,ncol = J, nrow = n) - - for(i in 1:n){ - for(j in 1:J){ - temp_mean <- exp(X[i,,drop = FALSE]%*%b[,j,drop = FALSE] + z[i]) - Y[i,j] <- rpois(1, lambda = temp_mean) - } - } - + + Y <- radEmu:::simulate_data(n = n, + J = J, + X = X, + b0 = b0, + b1 = b1, + distn = "Poisson", + mean_z = 10) ml_fit <- emuFit_micro(X, Y, @@ -160,12 +239,8 @@ test_that("ML fit to simple example give reasonable output with J >> n", { maxit = 500, tolerance = 1e-3, verbose = FALSE) - - + expect_true(max(abs(ml_fit[2,] - b1))<.5) - - - }) # # diff --git a/tests/testthat/test-emuFit_micro_penalized.R b/tests/testthat/test-emuFit_micro_penalized.R index e97a752..2cd0f60 100644 --- a/tests/testthat/test-emuFit_micro_penalized.R +++ b/tests/testthat/test-emuFit_micro_penalized.R @@ -48,27 +48,24 @@ test_that("Penalized estimation reduces to Haldane correction in saturated case test_that("PL fit to simple example returns reasonable values", { set.seed(4323) X <- cbind(1,rep(c(0,1),each = 20)) - z <- rnorm(40) +10 J <- 10 - p <- 2 n <- 40 - b0 <- rnorm(J) - b1 <- seq(1,10,length.out = J) - b <- rbind(b0,b1) - Y <- matrix(NA,ncol = J, nrow = 40) - - for(i in 1:40){ - for(j in 1:J){ - temp_mean <- exp(X[i,,drop = FALSE]%*%b[,j,drop = FALSE] + z[i]) - Y[i,j] <- rpois(1, lambda = temp_mean) - } - } + b1 <- seq(1,5,length.out = J) - + mean(seq(1,5,length.out = J)) + Y <- radEmu:::simulate_data(n = n, + J = J, + X = X, + b0 = rnorm(J), + b1 = b1, + distn = "Poisson", + mean_z = 50) + pl_fit <- emuFit_micro_penalized(X, Y, B = matrix(rnorm(20),nrow = 2), constraint_fn = function(x) mean(x), maxit = 200, - tolerance = 1e-5, + tolerance = 1e-10, verbose= FALSE) # plot(b1-mean(b1),ml_fit[2,]) @@ -82,31 +79,27 @@ test_that("PL fit to simple example returns numerically identical results regardless of whether we use computationally efficient augmentation or older less efficient implementation (and that both substantially differ from MLE", { set.seed(4323) - X <- cbind(1,rep(c(0,1),each = 5)) - z <- rnorm(10) +5 J <- 10 - p <- 2 n <- 10 - b0 <- rnorm(J) - b1 <- seq(1,5,length.out = J) - b1 <- b1 - mean(b1) - b <- rbind(b0,b1) - Y <- matrix(NA,ncol = J, nrow = n) - - for(i in 1:n){ - for(j in 1:J){ - temp_mean <- exp(X[i,,drop = FALSE]%*%b[,j,drop = FALSE] + z[i]) - Y[i,j] <- rnbinom(1, mu= temp_mean,size = 2)*rbinom(1,1,0.4) - } - } - - ml_fit <- emuFit_micro(X, - Y, + X <- cbind(1,rep(c(0,1),each = n/2)) + Y <- radEmu:::simulate_data(n = n, + J = J, + X = X, + b0 = rnorm(J), + b1 = seq(1,5,length.out = J), + distn = "ZINB", + zinb_size = 2, + zinb_zero_prop = 0.7, + mean_z = 10) + + ml_fit <- emuFit_micro(X = as.matrix(X), + Y = as.matrix(Y), B = NULL, constraint_fn = function(x) mean(x), maxit = 1000, tolerance = 1e-3, verbose= FALSE) + pl_fit_new <- emuFit_micro_penalized(X, Y, B = NULL, @@ -121,7 +114,7 @@ less efficient implementation (and that both substantially differ from MLE", { Y, B = NULL, constraint_fn = function(x) mean(x), - maxit = 1000, + maxit = 10000, tolerance = 1e-3, verbose= TRUE, use_legacy_augmentation = TRUE))) @@ -138,22 +131,17 @@ regardless of whether we use computationally efficient augmentation or older less efficient implementation (and that both substantially differ from MLE", { set.seed(4323) X <- cbind(1,rnorm(10)) - z <- rnorm(10) +5 J <- 10 - p <- 2 n <- 10 - b0 <- rnorm(J) - b1 <- seq(1,5,length.out = J) - b1 <- b1 - mean(b1) - b <- rbind(b0,b1) - Y <- matrix(NA,ncol = J, nrow = n) - - for(i in 1:n){ - for(j in 1:J){ - temp_mean <- exp(X[i,,drop = FALSE]%*%b[,j,drop = FALSE] + z[i]) - Y[i,j] <- rnbinom(1, mu= temp_mean,size = 2)*rbinom(1,1,0.2) - } - } + Y <- radEmu:::simulate_data(n = n, + J = J, + X = X, + b0 = rnorm(J), + b1 = seq(1,5,length.out = J), + distn = "ZINB", + zinb_size = 2, + zinb_zero_prop = 0.7, + mean_z = 10) ml_fit <- emuFit_micro(X, Y, @@ -162,6 +150,7 @@ less efficient implementation (and that both substantially differ from MLE", { maxit = 10000, tolerance = 0.01, verbose= FALSE) + pl_fit_new <- emuFit_micro_penalized(X, Y, B = NULL, diff --git a/tests/testthat/test-emuFit_phyloseq.R b/tests/testthat/test-emuFit_phyloseq.R index a676694..cf37957 100644 --- a/tests/testthat/test-emuFit_phyloseq.R +++ b/tests/testthat/test-emuFit_phyloseq.R @@ -13,9 +13,14 @@ test_that("emuFit works with a phyloseq object", { wirbel_smaller <- phyloseq::prune_samples(wirbel_sample$Country == "FRA" & wirbel_sample$Gender == "F", wirbel_small) - fit <- emuFit(wirbel_small, formula = ~ Group, run_score_tests = FALSE) + fit <- emuFit(wirbel_small, formula = ~ Group, run_score_tests = FALSE, tolerance = 0.01) expect_true(is.matrix(fit$B)) + wirbel_transpose <- wirbel_small + phyloseq::otu_table(wirbel_transpose) <- t(phyloseq::otu_table(wirbel_transpose)) + fit_transpose <- emuFit(wirbel_transpose, formula = ~ Group, run_score_tests = FALSE, tolerance = 0.01) + expect_true(all.equal(fit$coef, fit_transpose$coef)) + } else { expect_error(stop("You don't have phyloseq installed")) } diff --git a/tests/testthat/test-f_info.R b/tests/testthat/test-f_info.R index 88d20ae..7dd5197 100644 --- a/tests/testthat/test-f_info.R +++ b/tests/testthat/test-f_info.R @@ -33,17 +33,21 @@ test_that("computing information matrix gives same result regardless of method", test_that("Computed information is equal to numerical derivative with categorical predictor", { set.seed(59542234) n <- 2 - p <- 2 - X <- cbind(1,rep(c(0,1),each = n/2)) J <- 2 - z <- rnorm(n) +3 + X <- cbind(1,rep(c(0,1),each = n/2)) b0 <- rnorm(J) b1 <- seq(1,10,length.out = J) b1 <- b1 - mean(b1) b0 <- b0 - mean(b0) - b <- rbind(b0,b1) - Y <- matrix(NA,ncol = J, nrow = n) - + Y <- radEmu:::simulate_data(n = n, + J = J, + X = X, + b0 = b0, + b1 = b1, + distn = "Poisson", + mean_z = 3) + + k_constr <- 2 j_constr <- 1 p <- 2 @@ -66,17 +70,6 @@ test_that("Computed information is equal to numerical derivative with categorica maxit = 1000 inner_maxit = 25 - Y[] <- 0 - for(i in 1:n){ - while(sum(Y[i,])==0){ - for(j in 1:J){ - temp_mean <- exp(X[i,,drop = FALSE]%*%b[,j,drop = FALSE] + z[i]) - Y[i,j] <- rpois(1, lambda = temp_mean) - # Y[i,j] <- rnbinom(1,mu = temp_mean, size = 3)*rbinom(1,1,0.6) - } - } - } - full_fit <- #suppressMessages( emuFit_micro_penalized(X = X, Y = Y, diff --git a/tests/testthat/test-fit_null.R b/tests/testthat/test-fit_null.R index fce8ec2..209b005 100644 --- a/tests/testthat/test-fit_null.R +++ b/tests/testthat/test-fit_null.R @@ -1,16 +1,19 @@ test_that("we get same null fit with different j_ref", { set.seed(59542234) n <- 10 - p <- 2 - X <- cbind(1,rep(c(0,1),each = n/2)) J <- 5 - z <- rnorm(n) +8 + X <- cbind(1,rep(c(0,1),each = n/2)) b0 <- rnorm(J) - b1 <- seq(1,5,length.out = J) + b1 <- seq(1,10,length.out = J) b1 <- b1 - mean(b1) b0 <- b0 - mean(b0) - b <- rbind(b0,b1) - Y <- matrix(NA,ncol = J, nrow = n) + Y <- radEmu:::simulate_data(n = n, + J = J, + X = X, + b0 = b0, + b1 = b1, + distn = "Poisson", + mean_z = 8) k_constr <- 2 j_constr <- 1 @@ -34,17 +37,6 @@ test_that("we get same null fit with different j_ref", { maxit = 1000 inner_maxit = 25 - Y[] <- 0 - for(i in 1:n){ - while(sum(Y[i,])==0){ - for(j in 1:J){ - temp_mean <- exp(X[i,,drop = FALSE]%*%b[,j,drop = FALSE] + z[i]) - Y[i,j] <- rpois(1, lambda = temp_mean) - # Y[i,j] <- rnbinom(1,mu = temp_mean, size = 3)*rbinom(1,1,0.6) - } - } - } - full_fit <- #suppressMessages( emuFit_micro_penalized(X = X, Y = Y, diff --git a/tests/testthat/test-fit_null_repar.R b/tests/testthat/test-fit_null_repar.R new file mode 100644 index 0000000..5d78ddb --- /dev/null +++ b/tests/testthat/test-fit_null_repar.R @@ -0,0 +1,238 @@ +test_that("we get same null fit using fit_null and fit_null_repar", { + set.seed(59542234) + n <- 10 + p <- 2 + X <- cbind(1,rep(c(0,1),each = n/2)) + J <- 5 + z <- rnorm(n) +8 + b0 <- rnorm(J) + b1 <- seq(1,5,length.out = J) + b1 <- b1 - mean(b1) + b0 <- b0 - mean(b0) + b <- rbind(b0,b1) + Y <- matrix(NA,ncol = J, nrow = n) + + k_constr <- 2 + j_constr <- 1 + p <- 2 + + # constraint_fn <- function(x){ pseudohuber_center(x,0.1)} + constraint_fn <- function(x){pseudohuber_center(x,0.1)} + + ##### Arguments to fix: + + constraint_grad_fn <- function(x){dpseudohuber_center_dx(x,0.1)} + + rho_init = 1 + tau = 1.2 + kappa = 0.8 + obj_tol = 100 + B_tol <- 1e-3 + constraint_tol = 1e-5 + init_tol = 1e6 + c1 = 1e-4 + maxit = 1000 + inner_maxit = 25 + + Y[] <- 0 + for(i in 1:n){ + while(sum(Y[i,])==0){ + for(j in 1:J){ + temp_mean <- exp(X[i,,drop = FALSE]%*%b[,j,drop = FALSE] + z[i]) + Y[i,j] <- rpois(1, lambda = temp_mean) + # Y[i,j] <- rnbinom(1,mu = temp_mean, size = 3)*rbinom(1,1,0.6) + } + } + } + + full_fit <- #suppressMessages( + emuFit_micro_penalized(X = X, + Y = Y, + B = NULL, + constraint_fn = mean, + tolerance = 1e-7, + verbose = FALSE, + j_ref = 5) + + B <- full_fit$B + Y_aug <- full_fit$Y_augmented + + X_cup <- X_cup_from_X(X,J) + + null_fit <- fit_null(B = B, + Y = Y_aug, + X = X , + X_cup = X_cup, + k_constr = k_constr, + j_constr = j_constr, + j_ref = 5, + constraint_fn = constraint_fn, + constraint_tol = 1e-4, + B_tol = 1e-4, + constraint_grad_fn = constraint_grad_fn, + verbose = FALSE, + trackB = FALSE) + + null_fit_repar <- + fit_null_repar(B = B, + Y = Y_aug, + X = X, + k_constr = k_constr, + j_constr = j_constr, + j_ref =5, + constraint_fn = constraint_fn, + constraint_grad_fn = constraint_grad_fn, + maxit = 1000, + verbose = FALSE, + trackB = TRUE, + tolerance = 1e-2, + method = "fisher", + starting_stepsize = 0.5) + + expect_equal(null_fit_repar$B,null_fit$B,tolerance = 0.1) +}) + + + + +test_that("Vast majority of deviation from zero in null fits computed via +fit_null and fit_null_repar (with continuous covariate) are explained by +derivative of constraint function (as they should be at a constrained +approximate optimum).", { + set.seed(59542234) + n <- 10 + p <- 2 + X <- cbind(1,rnorm(10)) + J <- 25 + z <- rnorm(n) +8 + b0 <- rnorm(J) + b1 <- seq(1,5,length.out = J) + b1 <- b1 - mean(b1) + # b1[1] <- pseudohuber_center(b1[-1],0.1) + b0 <- b0 - mean(b0) + b <- rbind(b0,b1) + Y <- matrix(NA,ncol = J, nrow = n) + + k_constr <- 2 + j_constr <- 1 + p <- 2 + + constraint_fn <- function(x){ pseudohuber_center(x,0.1)} + # constraint_fn <- function(x){pseudohuber_center(x,0.1)} + # constraint_fn <- function(x){mean(x)} + + ##### Arguments to fix: + + constraint_grad_fn <- function(x){dpseudohuber_center_dx(x,0.1)} + # constraint_grad_fn <- function(x){rep(1/length(x),length(x))} + + rho_init = 1 + tau = 1.2 + kappa = 0.8 + obj_tol = 100 + B_tol <- 1e-3 + constraint_tol = 1e-5 + init_tol = 1e6 + c1 = 1e-4 + maxit = 1000 + inner_maxit = 25 + + Y[] <- 0 + for(i in 1:n){ + while(sum(Y[i,])==0){ + for(j in 1:J){ + temp_mean <- exp(X[i,,drop = FALSE]%*%b[,j,drop = FALSE] + z[i]) + Y[i,j] <- rpois(1, lambda = temp_mean) + # Y[i,j] <- rnbinom(1,mu = temp_mean, size = 3)*rbinom(1,1,0.6) + } + } + } + + full_fit <- #suppressMessages( + emuFit_micro_penalized(X = X, + Y = Y, + B = NULL, + constraint_fn = mean, + tolerance = 1e-6, + verbose = FALSE, + j_ref = 25) + + B <- full_fit$B + Y_aug <- full_fit$Y_augmented + + X_cup <- X_cup_from_X(X,J) + + null_fit <- fit_null(B = B, + Y = Y_aug, + X = X , + X_cup = X_cup, + k_constr = k_constr, + j_constr = j_constr, + j_ref = 25, + constraint_fn = constraint_fn, + constraint_tol = 1e-3, + B_tol = 1e-2, + rho_init = 100, + kappa = kappa, + tau = 1.2, + inner_tol = 0.01, + inner_maxit = 50, + maxit = 1e4, + constraint_grad_fn = constraint_grad_fn, + verbose = FALSE, + trackB = TRUE) + + null_fit_repar <- + fit_null_repar(B = B, + Y = Y_aug, + X = X, + k_constr = k_constr, + j_constr = j_constr, + j_ref =25, + constraint_fn = constraint_fn, + constraint_grad_fn = constraint_grad_fn, + maxit = 10000, + verbose = FALSE, + trackB = TRUE, + tolerance = 1e-4, + method = "fisher", + starting_stepsize = 1, + do_shift = TRUE) + + + + null_fit_grad <- numDeriv::grad(function(x){ + profile_ll(X = X, Y = Y_aug, B = cbind(rbind(x[1:(J - 1)], + x[(J -1) + 1:(J - 1)]), + 0)) + }, + c(null_fit$B[1,1:(J - 1)],null_fit$B[2,1:(J - 1)]) + ) + + null_fit_repar_grad <- numDeriv::grad(function(x){ + profile_ll(X = X, Y = Y_aug, B = cbind(rbind(x[1:(J - 1)], + x[(J -1) + 1:(J - 1)]), + 0)) + }, + c(null_fit_repar$B[1,1:(J - 1)],null_fit_repar$B[2,1:(J - 1)]) + ) + + cg_term <- constraint_grad_fn(null_fit$B[2,]) + cg_term[1] <- cg_term[1] - 1 + cg_term <- cg_term[1:24] + nf_b2 <- null_fit_grad[25:48] + cg_lm <- lm( nf_b2~cg_term - 1) + + + + cg_repar_term <- constraint_grad_fn(null_fit_repar$B[2,]) + cg_repar_term[1] <- cg_repar_term[1] - 1 + cg_repar_term <- cg_repar_term[1:24] + + nfr_b2 <- null_fit_repar_grad[25:48] + cg_repar_lm <- lm( nfr_b2~cg_repar_term - 1) + + + expect_equal(summary(cg_lm)$r.squared,1,tolerance = 0.0001) + expect_equal(summary(cg_repar_lm)$r.squared,1,tolerance = 0.0001) +}) diff --git a/tests/testthat/test-get_score_stat.R b/tests/testthat/test-get_score_stat.R index b41498e..be272bd 100644 --- a/tests/testthat/test-get_score_stat.R +++ b/tests/testthat/test-get_score_stat.R @@ -1,17 +1,19 @@ test_that("Robust score statistic is invariant to reference taxon", { set.seed(59542234) n <- 10 - p <- 2 - X <- cbind(1,rep(c(0,1),each = n/2)) J <- 5 - z <- rnorm(n) +8 + X <- cbind(1,rep(c(0,1),each = n/2)) b0 <- rnorm(J) b1 <- seq(1,10,length.out = J) b1 <- b1 - mean(b1) b0 <- b0 - mean(b0) - b <- rbind(b0,b1) - Y <- matrix(NA,ncol = J, nrow = n) - + Y <- radEmu:::simulate_data(n = n, + J = J, + X = X, + b0 = b0, + b1 = b1, + distn = "Poisson", + mean_z = 8) k_constr <- 2 j_constr <- 1 p <- 2 @@ -34,17 +36,6 @@ test_that("Robust score statistic is invariant to reference taxon", { maxit = 1000 inner_maxit = 25 - Y[] <- 0 - for(i in 1:n){ - while(sum(Y[i,])==0){ - for(j in 1:J){ - temp_mean <- exp(X[i,,drop = FALSE]%*%b[,j,drop = FALSE] + z[i]) - Y[i,j] <- rpois(1, lambda = temp_mean) - # Y[i,j] <- rnbinom(1,mu = temp_mean, size = 3)*rbinom(1,1,0.6) - } - } - } - full_fit <- #suppressMessages( emuFit_micro_penalized(X = X, Y = Y, @@ -96,7 +87,7 @@ test_that("Robust score statistic is invariant to reference taxon", { n = n, p = p, I = NULL, - Dy = NULL) + Dy = NULL)$score_stat ) expect_true(sd(score_stats)<1e-5) @@ -112,17 +103,19 @@ under null when Poisson assumption is met", { skip("Skipping test that requires 1000 simulations be run") set.seed(595434) n <- 10 - p <- 2 - X <- cbind(1,rep(c(0,1),each = n/2)) J <- 25 - z <- rnorm(n) +8 + X <- cbind(1,rep(c(0,1),each = n/2)) b0 <- rnorm(J) b1 <- seq(1,10,length.out = J) - # b1 <- c(-3,-2,0,5,0) b1 <- b1 - mean(b1) - # b0 <- b0 - mean(b0) - b <- rbind(b0,b1) - Y <- matrix(NA,ncol = J, nrow = n) + b0 <- b0 + Y <- radEmu:::simulate_data(n = n, + J = J, + X = X, + b0 = b0, + b1 = b1, + distn = "Poisson", + mean_z = 8) k_constr <- 2 j_constr <- 13 @@ -149,16 +142,15 @@ under null when Poisson assumption is met", { nsim <- 1000 score_stats <- numeric(nsim) for(sim in 1:nsim){ - Y[] <- 0 - for(i in 1:n){ - while(sum(Y[i,])==0){ - for(j in 1:J){ - temp_mean <- exp(X[i,,drop = FALSE]%*%b[,j,drop = FALSE] + z[i]) - # Y[i,j] <- rpois(1, lambda = temp_mean) - Y[i,j] <- rnbinom(1,mu = temp_mean, size = 3)*rbinom(1,1,0.6) - } - } - } + Y <- radEmu:::simulate_data(n = n, + J = J, + X = X, + b0 = b0, + b1 = b1, + distn = "ZINB", + zinb_size = 3, + zinb_zero_prop = 0.6, + mean_z = 8) full_fit <- #suppressMessages( emuFit_micro_penalized(X = X, @@ -193,7 +185,7 @@ under null when Poisson assumption is met", { c1 = c1, maxit = maxit, inner_maxit = inner_maxit, - verbose = FALSE) + verbose = FALSE)$score_stat j_ref <- 5 @@ -212,7 +204,7 @@ under null when Poisson assumption is met", { check_influence = FALSE, I = NULL, Dy = NULL, - model_based = FALSE) + model_based = FALSE)$score_stat x <- seq(-10,5,0.01) # hist(log(score_stats[1:sim]),breaks = 20,freq = FALSE) @@ -233,17 +225,19 @@ under null when Poisson assumption is met", { test_that("model-based score statistic is invariant to reference taxon", { set.seed(59542234) n <- 10 - p <- 2 - X <- cbind(1,rep(c(0,1),each = n/2)) J <- 5 - z <- rnorm(n) +8 + X <- cbind(1,rep(c(0,1),each = n/2)) b0 <- rnorm(J) - b1 <- seq(1,3,length.out = J) - # b1 <- c(-3,-2,0,5,0) + b1 <- seq(1,10,length.out = J) b1 <- b1 - mean(b1) - # b0 <- b0 - mean(b0) - b <- rbind(b0,b1) - Y <- matrix(NA,ncol = J, nrow = n) + b0 <- b0 + Y <- radEmu:::simulate_data(n = n, + J = J, + X = X, + b0 = b0, + b1 = b1, + distn = "Poisson", + mean_z = 8) k_constr <- 2 j_constr <- 3 @@ -267,18 +261,6 @@ test_that("model-based score statistic is invariant to reference taxon", { maxit = 1000 inner_maxit = 25 - - Y[] <- 0 - for(i in 1:n){ - while(sum(Y[i,])==0){ - for(j in 1:J){ - temp_mean <- exp(X[i,,drop = FALSE]%*%b[,j,drop = FALSE] + z[i]) - Y[i,j] <- rpois(1, lambda = temp_mean) - # Y[i,j] <- rnbinom(1,mu = temp_mean, size = 3)*rbinom(1,1,0.6) - } - } - } - full_fit <- #suppressMessages( emuFit_micro_penalized(X = X, Y = Y, @@ -329,7 +311,7 @@ test_that("model-based score statistic is invariant to reference taxon", { n = n, p = p, I = NULL, - Dy = NULL) + Dy = NULL)$score_stat }) expect_true(sd(score_stats)<1e-8) diff --git a/tests/testthat/test-macro_fisher_null.R b/tests/testthat/test-macro_fisher_null.R index f67a2ea..077a6d6 100644 --- a/tests/testthat/test-macro_fisher_null.R +++ b/tests/testthat/test-macro_fisher_null.R @@ -107,7 +107,7 @@ test_that("We take same step as we'd take using numerical derivatives when gap, min_ratio <- min(update$update/n_update,na.rm = TRUE) - expect_equal(max_ratio,min_ratio,tolerance = 1e-4) + expect_equal(max_ratio,min_ratio,tolerance = 1e-3) }) diff --git a/tests/testthat/test-micro_wald.R b/tests/testthat/test-micro_wald.R index fd4e90b..320b8f1 100644 --- a/tests/testthat/test-micro_wald.R +++ b/tests/testthat/test-micro_wald.R @@ -1,17 +1,76 @@ -test_that("wald test gives semi-reasonable output", { +# Here, we manually insert data used in test before the simulate_data() function +# was changed. This is a temporary solution until we want to do away with the +# specific p-value confirmations +Y_old1 <- + matrix( + c( 1748, 0, 0, 1287, 0, 0, 0, 0, 1875, 0, + 8286, 124134, 0, 1716, 0, 0, 4529, 3564, 21975, 0, + 4096, 43569, 0, 0, 0, 0, 2885, 2250, 1685, 217, + 4289, 122134, 416, 0, 5936, 0, 0, 0, 0, 0, + 1122, 15785, 920, 1354, 0, 0, 0, 0, 1654, 0, + 30007, 99540, 0, 8783, 0, 0, 13382, 20486, 28670, 0, + 5087, 0, 0, 3040, 0, 3954, 11802, 0, 9331, 0, + 3841, 41104, 1931, 6271, 33557, 1338, 0, 3442, 0, 0, + 3059, 0, 1279, 2274, 11459, 886, 2144, 5133, 4765, 0, + 3105, 0, 0, 0, 0, 426, 2622, 5103, 0, 0, + 80, 0, 0, 26431, 0, 0, 92214, 299927, 0, 0, + 32, 0, 0, 5186, 4065, 0, 18326, 34273, 0, 76950, + 0, 0, 2, 4147, 0, 0, 6183, 17407, 76158, 0, + 20, 572, 21, 0, 5391, 496, 11737, 0, 0, 19911, + 13, 0, 49, 0, 6721, 0, 0, 17896, 85495, 33042, + 0, 1497, 41, 6450, 8997, 709, 0, 149402, 247396, 43690, + 0, 0, 0, 0, 9225, 508, 12808, 42592, 51659, 68014, + 30, 0, 89, 0, 13951, 840, 7604, 0, 93587, 43885, + 41, 314, 0, 1483, 4061, 680, 4684, 0, 0, 6086, + 54, 0, 85, 0, 3871, 0, 9348, 25395, 84277, 0), + nrow = 20, ncol = 10, byrow = TRUE) + +Y_old2 <- + matrix( + c(6051354, 628305, 417132, 0, 0, 20886, 0, 237, 0, 54, + 190080, 0, 0, 327, 234, 5265, 106, 177, 0, 45, + 0, 0, 37588, 578, 933, 16147, 0, 691, 906, 1277, + 0, 0, 0, 0, 0, 1788, 1122, 0, 7657, 34478, + 0, 0, 16470, 9412, 0, 23456, 0, 17289, 0, 0, + 0, 1459, 0, 1479, 3321, 14736, 1319, 0, 0, 0, + 13, 13, 0, 324, 1832, 0, 0, 0, 0, 0, + 0, 22292, 21209, 0, 0, 2387, 139, 109, 6, 0, + 0, 17457, 0, 715, 1298, 0, 245, 0, 0, 127, + 0, 2562, 0, 0, 0, 1752, 522, 0, 0, 0, + 153, 113, 226, 0, 566, 2176, 861, 0, 0, 0, + 0, 6512, 110014, 10616, 6482, 86103, 2774, 18807, 7569, 8734, + 457, 0, 1240, 351, 401, 0, 3035, 10084, 4080, 14265, + 2631304, 255618, 0, 2257, 0, 0, 1351, 798, 0, 0, + 0, 0, 0, 0, 0, 0, 145, 0, 5, 0, + 110612, 8278, 0, 1309, 0, 5002, 0, 1018, 0, 0, + 0, 0, 0, 0, 5779, 230361, 0, 431527, 450787, 0, + 0, 0, 22245, 139, 0, 0, 151, 0, 0, 0, + 180, 169, 2659, 1115, 0, 6693, 12595, 125758, 107518, 425831, + 0, 0, 20670, 1216, 2086, 13617, 3777, 0, 1562, 4535), + nrow = 20, ncol = 10, byrow = TRUE) + + +test_that("wald test gives semi-reasonable output with categorical covariate", { set.seed(343234) n <- 20 X <- cbind(1,rep(c(0,1),each = n/2)) J <- 10 - z <- rnorm(n) +8 b0 <- rnorm(10) b1 <- 1:10 b1 <- b1 - mean(b1) b1[5] <- pseudohuber_center(b1[-5],0.1) b0 <- b0 - mean(b0) - b <- rbind(b0,b1) - Y <- matrix(NA,ncol = 10, nrow = n) + b <- rbind(b0, b1) + Y <- radEmu:::simulate_data(n = n, + J = J, + X = X, + b0 = b0, + b1 = b1, + distn = "ZINB", + zinb_size = 3, + zinb_zero_prop = 0.6, + mean_z = 8) k_constr <- 2 j_constr <- 5 @@ -28,18 +87,8 @@ test_that("wald test gives semi-reasonable output", { X_cup <- X_cup_from_X(X,J) - Y[] <- 0 - for(i in 1:n){ - while(sum(Y[i,])==0){ - for(j in 1:10){ - temp_mean <- exp(X[i,,drop = FALSE]%*%b[,j,drop = FALSE] + z[i]) - Y[i,j] <- rnbinom(1,mu = temp_mean, size = 3)*rbinom(1,1,0.6) - } - } - } - full_fit <- emuFit_micro_penalized(X = X, - Y = Y, + Y = Y_old1, B = NULL, constraint_fn = constraint_fn, tolerance = 1e-3, @@ -61,6 +110,30 @@ test_that("wald test gives semi-reasonable output", { expect_true(ncol(wald_result$I) ==20) expect_equal(wald_result$coefficients$pval, 0.61, tolerance = 0.02) + + + full_fit_new <- emuFit_micro_penalized(X = X, + Y = Y, + B = NULL, + constraint_fn = constraint_fn, + tolerance = 1e-3, + verbose = FALSE) + + wald_result_new <- + micro_wald(Y = full_fit_new$Y, + X, + X_cup = X_cup, + B = full_fit_new$B, + test_kj = data.frame(k = 2, j = 4), + constraint_fn = constraint_fn, + constraint_grad_fn = constraint_grad_fn, + nominal_coverage = 0.95) + + expect_true(is.data.frame(wald_result_new$coefficients)) + expect_true(wald_result_new$coefficients$pval>0.1) + expect_true(is.list(wald_result_new)) + expect_true(ncol(wald_result_new$I) ==20) + }) @@ -70,14 +143,21 @@ test_that("wald test gives semi-reasonable output with continuous covariate", { n <- 20 X <- cbind(1,rnorm(n)) J <- 10 - z <- rnorm(n) +8 b0 <- rnorm(10) b1 <- 1:10 b1 <- b1 - mean(b1) b1[5] <- pseudohuber_center(b1[-5],0.1) b0 <- b0 - mean(b0) - b <- rbind(b0,b1) - Y <- matrix(NA,ncol = 10, nrow = n) + b <- rbind(b0, b1) + Y <- radEmu:::simulate_data(n = n, + J = J, + X = X, + b0 = b0, + b1 = b1, + distn = "ZINB", + zinb_size = 3, + zinb_zero_prop = 0.6, + mean_z = 8) k_constr <- 2 j_constr <- 5 @@ -94,18 +174,8 @@ test_that("wald test gives semi-reasonable output with continuous covariate", { X_cup <- X_cup_from_X(X,J) - Y[] <- 0 - for(i in 1:n){ - while(sum(Y[i,])==0){ - for(j in 1:10){ - temp_mean <- exp(X[i,,drop = FALSE]%*%b[,j,drop = FALSE] + z[i]) - Y[i,j] <- rnbinom(1,mu = temp_mean, size = 3)*rbinom(1,1,0.6) - } - } - } - full_fit <- emuFit_micro_penalized(X = X, - Y = Y, + Y = Y_old2, B = NULL, constraint_fn = constraint_fn, tolerance = 1e-3, @@ -134,7 +204,39 @@ test_that("wald test gives semi-reasonable output with continuous covariate", { expect_true(is.list(wald_result)) expect_true(ncol(wald_result$I) ==20) expect_equal(wald_result$coefficients$pval, 0.11, tolerance = 0.03) - expect_true(wald_result_for_an_alternative$coefficients$pval < 0.01) + + + full_fit_new <- emuFit_micro_penalized(X = X, + Y = Y, + B = NULL, + constraint_fn = constraint_fn, + tolerance = 1e-3, + verbose = FALSE) + + wald_result_new <- micro_wald(Y = full_fit_new$Y, + X, + X_cup = X_cup, + B = full_fit_new$B, + test_kj = data.frame(k = 2, j = 4), + constraint_fn = constraint_fn, + constraint_grad_fn = constraint_grad_fn, + nominal_coverage = 0.95) + + wald_result_for_an_alternative_new <- micro_wald(Y = full_fit_new$Y, + X, + X_cup = X_cup, + B = full_fit_new$B, + test_kj = data.frame(k = 2, j = 10), + constraint_fn = constraint_fn, + constraint_grad_fn = constraint_grad_fn, + nominal_coverage = 0.95) + + + expect_true(is.data.frame(wald_result_new$coefficients)) + expect_true(is.list(wald_result_new)) + expect_true(ncol(wald_result_new$I) ==20) + expect_true(wald_result_for_an_alternative_new$coefficients$pval < 0.01) + }) diff --git a/tests/testthat/test-null_repar_one.R b/tests/testthat/test-null_repar_one.R new file mode 100644 index 0000000..865ec13 --- /dev/null +++ b/tests/testthat/test-null_repar_one.R @@ -0,0 +1,92 @@ +# test_that("we get same null fit using fit_null and fit_null_repar", { +# set.seed(59542234) +# n <- 10 +# p <- 2 +# X <- cbind(1,rep(c(0,1),each = n/2)) +# J <- 5 +# z <- rnorm(n) +8 +# b0 <- rnorm(J) +# b1 <- seq(1,5,length.out = J) +# b1 <- b1 - mean(b1) +# b0 <- b0 - mean(b0) +# b <- rbind(b0,b1) +# Y <- matrix(NA,ncol = J, nrow = n) +# +# k_constr <- 2 +# j_constr <- 1 +# p <- 2 +# +# # constraint_fn <- function(x){ pseudohuber_center(x,0.1)} +# constraint_fn <- function(x){pseudohuber_center(x,0.1)} +# +# ##### Arguments to fix: +# +# constraint_grad_fn <- function(x){dpseudohuber_center_dx(x,0.1)} +# # constraint_grad_fn <- function(x){ rep(1/length(x), length(x))} +# rho_init = 1 +# tau = 1.2 +# kappa = 0.8 +# obj_tol = 100 +# B_tol <- 1e-3 +# constraint_tol = 1e-5 +# init_tol = 1e6 +# c1 = 1e-4 +# maxit = 1000 +# inner_maxit = 25 +# +# Y[] <- 0 +# for(i in 1:n){ +# while(sum(Y[i,])==0){ +# for(j in 1:J){ +# temp_mean <- exp(X[i,,drop = FALSE]%*%b[,j,drop = FALSE] + z[i]) +# Y[i,j] <- rpois(1, lambda = temp_mean) +# # Y[i,j] <- rnbinom(1,mu = temp_mean, size = 3)*rbinom(1,1,0.6) +# } +# } +# } +# +# full_fit <- #suppressMessages( +# emuFit_micro_penalized(X = X, +# Y = Y, +# B = NULL, +# constraint_fn = mean, +# tolerance = 1e-7, +# verbose = FALSE, +# j_ref = 5) +# +# B <- full_fit$B +# Y_aug <- full_fit$Y_augmented +# +# X_cup <- X_cup_from_X(X,J) +# +# null_fit_repar <- +# fit_null_repar(B = B, +# Y = Y_aug, +# X = X, +# k_constr = k_constr, +# j_constr = j_constr, +# j_ref =5, +# constraint_fn = constraint_fn, +# constraint_grad_fn = constraint_grad_fn, +# maxit = 1000, +# verbose = TRUE, +# trackB = TRUE, +# tolerance = 1e-7, +# method = "exact", +# starting_stepsize = 0.5) +# +# a_dl_augs <- lapply(2:4, +# function(j){ +# d_l_aug_dB(Bj, +# Bj_constr_no_k_constr, +# z, +# Y, +# X, +# Bk_constr_no_j_j_constr, +# k_constr, +# j, +# j_constr, +# constraint_fn, +# constraint_grad_fn) +# +# }) diff --git a/tests/testthat/test-plot_emuFit.R b/tests/testthat/test-plot_emuFit.R new file mode 100644 index 0000000..a79dbb4 --- /dev/null +++ b/tests/testthat/test-plot_emuFit.R @@ -0,0 +1,53 @@ +set.seed(11) +J <- 6 +n <- 12 +X <- cbind(1,rnorm(n)) +Y <- radEmu:::simulate_data(n = n, + J = J, + X = X, + b0 = rnorm(10), + b1 = 1:10 - mean(1:10), + distn = "ZINB", + zinb_size = 3, + zinb_zero_prop = 0.6, + mean_z = 8) + +fitted_model <- emuFit(Y = Y, + X = X, + formula = ~group, + data = covariates, + verbose = FALSE, + B_null_tol = 1e-2, + tolerance = 0.01, + tau = 2, + return_wald_p = FALSE, + compute_cis = TRUE, + run_score_tests = FALSE, + use_fullmodel_info = FALSE, + use_fullmodel_cov = FALSE, + return_both_score_pvals = FALSE) + +test_that("plot() returns data frame and plot", { + plot_out <- plot(x = fitted_model) + expect_true(is.data.frame(plot_out$data)) + expect_true(all(sapply(plot_out$plots, ggplot2::is.ggplot))) +}) + + +test_that("plot() returns error when plot_key does not match coefficient table", { + expect_error({ + plot(x = fitted_model, + plot_key = list(c("First Covariate" = "covariate1"))) + }) +}) + +test_that("plot() returns error when coefficient is included multiple times in plot_key", { + expect_error({ + plot(x = fitted_model, + plot_key = list(c("First Covariate" = "covariate_1"), + c("Second Covariate" = "covariate_1"))) + }) +}) + + + diff --git a/tests/testthat/test-profile_ll.R b/tests/testthat/test-profile_ll.R new file mode 100644 index 0000000..b61c7bd --- /dev/null +++ b/tests/testthat/test-profile_ll.R @@ -0,0 +1,9 @@ +test_that("profile_ll returns a number", { + + set.seed(88442) + Y <- matrix(rexp(20),ncol = 5) + X <- cbind(1,rep(0:1,2)) + B <- matrix(rnorm(10),nrow = 2) + + expect_true(is.numeric(profile_ll(X = X, Y = Y, B = B))) +}) diff --git a/tests/testthat/test-row_name_matching.R b/tests/testthat/test-row_name_matching.R new file mode 100644 index 0000000..fbddd47 --- /dev/null +++ b/tests/testthat/test-row_name_matching.R @@ -0,0 +1,75 @@ +dat <- data.frame(cov1 = rep(c("A", "B", "C"), each = 6), + cov2 = rep(c("D", "E"), each = 9), + cov3 = rnorm(18), + cov4 = rnorm(18), + cov5 <- rep(c("G", "H", "I"), 6)) + +form <- ~ cov1 + cov2 + cov3 + cov4 + cov5 +X.based <- model.matrix(form, dat) + +Y <- matrix(rpois(18*6, 3), nrow = 18, ncol = 6) +colnames(Y) <- paste0("category_", 1:ncol(Y)) +rownames(Y) <- paste0("sample_", 1:nrow(Y)) + +test_that("emuFit handles missing row names", { + X1 <- matrix(X.based, nrow = 18) # No row names + + expect_message(emuFit(Y = Y, X = X1), + "Row names are missing from the covariate matrix X. We will assume the rows are in the same order as in the response matrix Y. You are responsible for ensuring the order of your observations is the same in both matrices.") +}) + +test_that("emuFit throws error on duplicate row names", { + X2 <- X.based + rownames(X2) <- rownames(Y) + rownames(X2)[5] <- "sample_4" #Repeating one of the sample labels + + expect_error(emuFit(Y = Y, X = X2), + "Covariate matrix X has duplicated row names. Please ensure all row names are unique before refitting the model.") +}) + +test_that("emuFit subsets to common row names with warning", { + X3 <- X.based + rownames(X3) <- rownames(Y) + X3 <- X3[c(1:4,7:14,16:18), , drop = FALSE] + Y3 <- Y[c(1:2,5:18), , drop = FALSE] + + expect_warning(emuFit(Y = Y3, X = X3), + regexp = "According to the rownames, there are observations that are missing either in the covariate matrix \\(X\\) and/or the response matrix \\(Y\\)\\. We will subset to common rows only, resulting in [0-9]+ samples\\.") +}) + +test_that("emuFit reorders rows of X when", { + X4 <- X.based + rownames(X4) <- rownames(Y) + X4.p <- X4[c(1,(nrow(Y):2)), , drop = FALSE] + + expect_message(model.p <- emuFit(Y = Y, X = X4.p), + "There is a different row ordering between the covariate matrix \\(X\\) and the response matrix \\(Y\\)\\. Covariate data will be reordered to match response data\\.") + model.o <- emuFit(Y = Y, X = X4) + expect_equal(model.o$coef$estimate, model.p$coef$estimate) +}) + +test_that("emuFit does not reorder rows of X when match_row_names is FALSE", { + X5 <- X.based + X5.p <- X.based + rownames(X5) <- rownames(Y) + rownames(X5.p) <- rownames(Y)[c(1,nrow(Y):2)] + + model.o <- emuFit(Y = Y, X = X5) + model.p <- emuFit(Y = Y, X = X5.p, match_row_names = FALSE) + + expect_silent(model.o <- emuFit(Y = Y, X = X5)) + expect_silent(model.p <- emuFit(Y = Y, X = X5.p, match_row_names = FALSE)) + + expect_equal(model.o$coef$estimate, model.p$coef$estimate) +}) + +test_that("emuFit stops when match_row_names is FALSE, but nrow does not coincide",{ + + X6 <- X.based + rownames(X6) <- rownames(Y) + X6 <- X6[c(1,(nrow(Y):2)),] + X6 <- X6[(1:16), , drop = FALSE] + + expect_error(emuFit(Y = Y, X = X6, match_row_names = FALSE), + "The number of rows does not match between the covariate matrix \\(X\\) and the response matrix \\(Y\\), and subsetting/matching by row name has been disabled\\. Please resolve this issue before refitting the model\\.") +}) diff --git a/tests/testthat/test-score_test.R b/tests/testthat/test-score_test.R index 27f88ef..92188b2 100644 --- a/tests/testthat/test-score_test.R +++ b/tests/testthat/test-score_test.R @@ -2,18 +2,24 @@ test_that("We get same score test results regardless of whether we provide inver we do *not* get same results if we use incorrect info", { set.seed(343234) - X <- cbind(1,rep(c(0,1),each = 20)) J <- 10 n <- 40 - p <- 2 - z <- rnorm(40) +8 - b0 <- rnorm(10) - b1 <- 1:10 + X <- cbind(1,rep(c(0,1),each = n/2)) + b0 <- rnorm(J) + b1 <- 1:J b1 <- b1 - mean(b1) b1[5] <- pseudohuber_center(b1[-5],0.1) b0 <- b0 - mean(b0) - b <- rbind(b0,b1) - Y <- matrix(NA,ncol = 10, nrow = 40) + b <- rbind(b0, b1) + Y <- radEmu:::simulate_data(n = n, + J = J, + X = X, + b0 = b0, + b1 = b1, + distn = "ZINB", + zinb_size = 3, + zinb_zero_prop = 0.6, + mean_z = 4) k_constr <- 2 j_constr <- 5 @@ -25,15 +31,6 @@ we do *not* get same results if we use incorrect info", { ##### Arguments to fix: - - for(i in 1:40){ - for(j in 1:10){ - temp_mean <- exp(X[i,,drop = FALSE]%*%b[,j,drop = FALSE] + z[i]) - # Y[i,j] <- rpois(1, lambda = temp_mean) - Y[i,j] <- rnbinom(1,mu = temp_mean, size = 3)*rbinom(1,1,0.6) - } - } - full_fit <- emuFit(X = X, Y = Y, B = NULL, @@ -87,7 +84,7 @@ we do *not* get same results if we use incorrect info", { J = J, n = n, p = p, - I_inv = I_inv) + I_inv = I_inv)$score_stat score_stat_with_other_mat <- get_score_stat(Y = Y_aug, @@ -102,7 +99,7 @@ we do *not* get same results if we use incorrect info", { J = J, n = n, p = p, - I_inv = diag(rep(1,18))) + I_inv = diag(rep(1,18)))$score_stat expect_equal(score_stat_with_I_inv,score_test_as_is$score_stat) expect_true(score_stat_with_other_mat !=score_test_as_is$score_stat) diff --git a/tests/testthat/test-simulate_data.R b/tests/testthat/test-simulate_data.R new file mode 100644 index 0000000..986dd62 --- /dev/null +++ b/tests/testthat/test-simulate_data.R @@ -0,0 +1,33 @@ +test_that("simulating data works", { + J <- 10; n <- 60 + # cluster membership + cluster <- rep(1:4, each = n/4) + # intercepts for each category + b0 <- rnorm(J) + # coefficients for X1 for each category + b1 <- seq(1, 5, length.out = J) + + dat1 <- radEmu:::simulate_data(n = n, J = J, b0 = b0, b1 = b1, distn = "Poisson", + mean_z = 8) + dat2 <- radEmu:::simulate_data(n = n, J = J, b0 = b0, b1 = b1, distn = "Poisson", + mean_z = 8, cluster = cluster) + dat3 <- radEmu:::simulate_data(n = n, J = J, b0 = b0, b1 = b1, distn = "ZINB", + mean_z = 8, zinb_size = 10, zinb_zero_prop = 0.5) + dat4 <- radEmu:::simulate_data(n = n, J = J, b0 = b0, b1 = b1, distn = "ZINB", cluster = cluster, + mean_z = 8, zinb_size = 10, zinb_zero_prop = 0.5) + + # make sure data are the correct dimensions + expect_true(nrow(dat1) == n & nrow(dat2) == n & nrow(dat3) == n & nrow(dat4) == n & + ncol(dat1) == J & ncol(dat2) == J & ncol(dat3) == J & ncol(dat4) == J) + # check that ZINB data is more sparse than Poisson data + expect_true(mean(dat3 == 0) > mean(dat1 == 0)) + expect_true(mean(dat4 == 0) > mean(dat2 == 0)) + +}) + +test_that("simulating b's works", { + b_res <- get_sim_bs(10) + expect_true(length(b_res[[1]]) == 10) + expect_true(length(b_res[[2]]) == 10) +}) + \ No newline at end of file diff --git a/tests/testthat/test-zero_comparison_check.R b/tests/testthat/test-zero_comparison_check.R new file mode 100644 index 0000000..62de43f --- /dev/null +++ b/tests/testthat/test-zero_comparison_check.R @@ -0,0 +1,109 @@ +dat <- data.frame(cov1 = rep(c("A", "B", "C"), each = 6), + cov2 = rep(c("D", "E"), each = 9), + cov3 = rnorm(18), + cov4 = rnorm(18), + cov5 <- rep(c("G", "H", "I"), 6)) +form1 <- ~ cov1 +form2 <- ~ cov2 +form3 <- ~ cov1 + cov3 +form4 <- ~ cov1 + cov3 + cov4 +form5 <- ~ cov1 + cov2 + cov3 + cov4 + cov5 +X1 <- model.matrix(form1, dat) +X1_base <- matrix(as.vector(X1), nrow = nrow(X1)) +colnames(X1_base) <- colnames(X1) +X2 <- model.matrix(form2, dat) +X2_base <- matrix(as.vector(X2), nrow = nrow(X2)) +colnames(X2_base) <- colnames(X2) +X3 <- model.matrix(form3, dat) +X3_base <- matrix(as.vector(X3), nrow = nrow(X3)) +colnames(X3_base) <- colnames(X3) +X4 <- model.matrix(form4, dat) +X4_base <- matrix(as.vector(X4), nrow = nrow(X4)) +colnames(X4_base) <- colnames(X4) +X5 <- model.matrix(form5, dat) +Y <- matrix(rpois(18*6, 3), nrow = 18, ncol = 6) +colnames(Y) <- paste0("category_", 1:ncol(Y)) +Y0 <- Y +Y0[, c(1, 4)] <- 0 +Y0[15, 1] <- 2 +Y0[10, 4] <- 3 + +### check functionality when X does not have "assign" attribute from `model.matrix` +test_that("zero_comparison_check function works", { + + zero_comparison_res <- zero_comparison_check(X = X1_base, Y = Y0) + expect_true(zero_comparison_res$zero_comparison[1]) + expect_true(zero_comparison_res$zero_comparison[10]) + expect_false(zero_comparison_res$zero_comparison[2]) + +}) + +test_that("zero_comparison column is added to coef when it should be", { + + # column doesn't exist with a 2 level covariate + emuRes1 <- emuFit(Y = Y0, X = X2_base, run_score_tests = FALSE, compute_cis = FALSE) + expect_false("zero_comparison" %in% names(emuRes1$coef)) + + # column doesn't exist when no zero-comparison parameters + emuRes2 <- emuFit(Y = Y, X = X1_base, run_score_tests = FALSE, compute_cis = FALSE) + expect_false("zero_comparison" %in% names(emuRes2$coef)) + + # column does exist when there are zero-comparison parameters + emuRes3 <- emuFit(Y = Y0, X = X1_base, run_score_tests = FALSE, compute_cis = FALSE) + expect_true("zero_comparison" %in% names(emuRes3$coef)) + + # column does exist in the presence of continuous covariate + emuRes4 <- emuFit(Y = Y0, X = X4_base, run_score_tests = FALSE, compute_cis = FALSE) + expect_true("zero_comparison" %in% names(emuRes4$coef)) + +}) + +test_that("remove_zero_comparison_pvals argument works in different ways", { + + expect_error(emuFit(Y = Y0, X = X1_base, remove_zero_comparison_pvals = 15)) + + emuRes1 <- emuFit(Y = Y0, X = X1_base, remove_zero_comparison_pvals = TRUE, + test_kj = data.frame(k = 2, j = 1)) + expect_true(is.na(emuRes1$coef$pval[1])) + emuRes2 <- emuFit(Y = Y0, X = X1_base, remove_zero_comparison_pvals = TRUE, + test_kj = data.frame(k = 2, j = 1), use_fullmodel_info = TRUE, + return_both_score_pvals = TRUE) + expect_true(is.na(emuRes2$coef$score_pval_full_info[1]) & + is.na(emuRes2$coef$score_pval_null_info[1])) + emuRes3 <- emuFit(Y = Y0, X = X1_base, remove_zero_comparison_pvals = 0.5, + test_kj = data.frame(k = 2, j = 1:2)) + expect_true(is.na(emuRes3$coef$pval[1]) || emuRes3$coef$pval[1] > 0.5) + expect_false(is.na(emuRes3$coef$pval[2])) + emuRes4 <- emuFit(Y = Y0, X = X1_base, remove_zero_comparison_pvals = FALSE, + test_kj = data.frame(k = 2, j = 1)) + expect_true(emuRes4$coef$zero_comparison[1] & !is.na(emuRes4$coef$pval[1])) + +}) + +### check functionality when X does have "assign" attribute from `model.matrix` +test_that("zero_comparison column is added to coef when it should be, model.matrix X", { + + # column doesn't exist with a 2 level covariate + emuRes1 <- emuFit(Y = Y0, X = X2, run_score_tests = FALSE, compute_cis = FALSE) + expect_false("zero_comparison" %in% names(emuRes1$coef)) + + # column doesn't exist when no zero-comparison parameters + emuRes2 <- emuFit(Y = Y, X = X1, run_score_tests = FALSE, compute_cis = FALSE) + expect_false("zero_comparison" %in% names(emuRes2$coef)) + + # column does exist when there are zero-comparison parameters + emuRes3 <- emuFit(Y = Y0, X = X1, run_score_tests = FALSE, compute_cis = FALSE) + expect_true("zero_comparison" %in% names(emuRes3$coef)) + + # column does exist in the presence of continuous covariate + emuRes4 <- emuFit(Y = Y0, X = X4, run_score_tests = FALSE, compute_cis = FALSE) + expect_true("zero_comparison" %in% names(emuRes4$coef)) + + # column does exist when there are multiple category covariates + emuRes5 <- emuFit(Y = Y0, X = X5, run_score_tests = FALSE, compute_cis = FALSE) + expect_true("zero_comparison" %in% names(emuRes5$coef)) + expect_true(emuRes5$coef$zero_comparison[1]) + expect_false(emuRes5$coef$zero_comparison[13]) + expect_true(emuRes5$coef$zero_comparison[31]) + expect_false(emuRes5$coef$zero_comparison[32]) +}) diff --git a/vignettes/.gitignore b/vignettes/.gitignore new file mode 100644 index 0000000..097b241 --- /dev/null +++ b/vignettes/.gitignore @@ -0,0 +1,2 @@ +*.html +*.R diff --git a/vignettes/intro_radEmu.Rmd b/vignettes/intro_radEmu.Rmd index bd04222..61812af 100644 --- a/vignettes/intro_radEmu.Rmd +++ b/vignettes/intro_radEmu.Rmd @@ -122,6 +122,18 @@ Again, while we would generally fit a model using all of our samples, for this t ch_study_obs <- which(wirbel_sample$Country %in% c("CHI")) ``` +Next, we want to confirm that all samples have at least one non-zero count across the categories we've chosen and that all categories have at least one non-zero count across the samples we've chosen. + +```{r} +small_Y <- wirbel_otu[ch_study_obs, restricted_mOTU_names] +sum(rowSums(small_Y) == 0) # no samples have a count sum of 0 +sum(colSums(small_Y) == 0) # one category has a count sum of 0 + +category_to_rm <- which(colSums(small_Y) == 0) +small_Y <- small_Y[, -category_to_rm] +sum(colSums(small_Y) == 0) +``` + The function that we use to fit our model is called `emuFit`. It can accept your data in various forms, and here we will show how to use it with data frames as input. Check out the `phyloseq` vignette if you'd like to know how `radEmu` plays with `phyloseq` objects! One version of inputs to `emuFit` are - `formula`: This is a formula telling radEmu what predictors to use in its model. We are using Group, which is an indicator for case (CRC) vs control (CTR). @@ -136,7 +148,7 @@ and some optional arguments include ```{r} ch_fit <- emuFit(formula = ~ Group, data = wirbel_sample[ch_study_obs, ], - Y = wirbel_otu[ch_study_obs, restricted_mOTU_names], + Y = small_Y, run_score_tests = FALSE) ``` Let's check out what this object looks like! @@ -146,30 +158,29 @@ ch_fit ``` -The way to access estimated coefficients and confidence intervals from the model is with `ch_fit$coef`. Let's plot our results: -```{r, fig.width= 7, fig.height = 4} -ch_df <- ch_fit$coef %>% - mutate(Genus = mOTU_name_df %>% - filter(genus_name %in% chosen_genera) %>% - pull(genus_name)) %>% - # add genus name to output from emuFit +Now, we can easily visualize our results using the `plot` function! +So that we only produce the plot, and not the dataframe used to produce the plots, +we will explicitly extract the plots using the `$` operator to get the `plots` component. + +```{r, fig.height = 6, fig.width = 6} +plot(ch_fit)$plots +``` + +In the plot above, it is a bit difficult to read the taxa names on the y-axis. We can create a data frame that maps the full taxa names in our data to simplified labels and plot using those instead, as shown below. +```{r, fig.height = 6, fig.width = 6} +#create data frame that has simplified names for taxa +taxa_names <- data.frame("category" = ch_fit$coef$category) %>% mutate(cat_small = stringr::str_remove(paste0("mOTU_", stringr::str_split(category, 'mOTU_v2_', simplify = TRUE)[, 2]), - "\\]")) %>% - mutate(cat_small = factor(cat_small, levels = cat_small[order(Genus)])) - # reorder mOTU categories by genus - -ggplot(ch_df) + - geom_point(aes(x = cat_small, y = estimate,color = Genus), size = .5) + - geom_errorbar(aes(x = cat_small, ymin = lower, ymax = upper, color = Genus), width = .25) + - theme_bw() + - theme(axis.text.x = element_text(angle = 45, hjust = 1)) + - labs(x = "Category", - y = "Estimate") + - coord_cartesian(ylim = c(-5,10)) + "\\]")) + + +#produce plot with cleaner taxa labels +plot(ch_fit, taxon_names = taxa_names)$plots ``` + Interestingly, we estimate a meta-mOTU "unknown Eubacterium [meta_mOTU_v2_7116]" assigned to Eubacteria to have a much higher ratio of abundance (comparing CRC group to control) than is typical across the mOTUs we included in this analysis. The confidence interval for this effect does not include zero -- but (!!!) the kind of confidence interval that is returned by default by emuFit is not extremely reliable when counts are very skewed or sample size is small-to-moderate. @@ -192,7 +203,7 @@ two_robust_score_tests <- emuFit(formula = ~ Group, B = ch_fit, test_kj = data.frame(k = covariate_to_test, j = taxa_to_test), - Y = as.matrix(wirbel_otu[ch_study_obs, restricted_mOTU_names])) + Y = small_Y) ``` Let's take a look at the test output. @@ -232,7 +243,7 @@ We could run robust score tests for every taxon in this analysis, but it will ta test_all <- emuFit(formula = ~ Group, data = wirbel_sample[ch_study_obs, ], B = ch_fit, - Y = wirbel_otu[ch_study_obs, restricted_mOTU_names], + Y = small_Y, run_score_tests = TRUE) ``` diff --git a/vignettes/intro_radEmu_with_phyloseq.Rmd b/vignettes/intro_radEmu_with_phyloseq.Rmd index 4690842..6f5e282 100644 --- a/vignettes/intro_radEmu_with_phyloseq.Rmd +++ b/vignettes/intro_radEmu_with_phyloseq.Rmd @@ -106,8 +106,7 @@ head(phyloseq::tax_table(wirbel_phylo)) `radEmu` is a package that can be used to estimate fold-differences in the abundance of microbial taxa between levels of a covariate. In this analysis, the covariate that we are primarily interested in is whether a sample is from a case of colorectal cancer or a control. We will make control ("CTR") the reference category: ```{r, eval = phy} -phyloseq::sample_data(wirbel_phylo)$Group <- factor(phyloseq::sample_data(wirbel_phylo)$Group, - levels = c("CTR","CRC")) +phyloseq::sample_data(wirbel_phylo)$Group <- factor(phyloseq::sample_data(wirbel_phylo)$Group, levels = c("CTR","CRC")) ``` While in general we would fit a model to all mOTUs, we are going to subset to some specific genera for the purposes of this tutorial. Let's look at *Eubacterium*, *Porphyromonas*, *Faecalibacteria*, and *Fusobacterium* for now. @@ -123,6 +122,16 @@ Again, while we would generally fit a model using all of our samples, for this t wirbel_china <- phyloseq::subset_samples(wirbel_restrict, Country == "CHI") ``` +Next, we want to confirm that all samples have at least one non-zero count across the categories we've chosen and that all categories have at least one non-zero count across the samples we've chosen. + +```{r} +sum(rowSums(phyloseq::otu_table(wirbel_china)) == 0) # no samples have a count sum of 0 +sum(colSums(phyloseq::otu_table(wirbel_china)) == 0) # one category has a count sum of 0 +category_to_rm <- names(which(colSums(phyloseq::otu_table(wirbel_china)) == 0)) +wirbel_china <- phyloseq::subset_taxa(wirbel_china, species != category_to_rm) +sum(colSums(phyloseq::otu_table(wirbel_china)) == 0) # now no categories have a count sum of 0 +``` + The function that we use to fit our model is called `emuFit`. It can accept your data in various forms, and here we will show how to use it with a `phyloseq` object as input. @@ -132,26 +141,12 @@ ch_fit <- emuFit(formula = ~ Group, run_score_tests = FALSE) ``` -The way to access estimated coefficients and confidence intervals from the model is with `ch_fit$coef`. Let's plot our results: - -```{r, fig.width= 7, fig.height = 4, eval = phy} -ch_df <- ch_fit$coef %>% - mutate(Genus = as.vector(phyloseq::tax_table(wirbel_china)[, 6])) %>% - # add genus name to output from emuFit - mutate(cat_small = stringr::str_remove(paste0("mOTU_", - stringr::str_split(category, 'mOTU_v2_', simplify = TRUE)[, 2]), - "\\]")) %>% - mutate(cat_small = factor(cat_small, levels = cat_small[order(Genus)])) - # reorder mOTU categories by genus - -ggplot(ch_df) + - geom_point(aes(x = cat_small, y = estimate, color = Genus), size = .5) + - geom_errorbar(aes(x = cat_small, ymin = lower, ymax = upper, color = Genus), width = .25) + - theme_bw() + - theme(axis.text.x = element_text(angle = 45, hjust = 1)) + - labs(x = "Category", - y = "Estimate") + - coord_cartesian(ylim = c(-5,10)) +The way to access estimated coefficients and confidence intervals from the model is with `ch_fit$coef`. + +Now, we can easily visualize our results using the `plot.emuFit` function! + +```{r, fig.height = 6, fig.width = 6} +plot(ch_fit)$plots ``` If you'd like to see more explanations of the `radEmu` software and additional analyses of this data, check out the vignette "intro_radEmu.Rmd". diff --git a/vignettes/parallel_radEmu.Rmd b/vignettes/parallel_radEmu.Rmd new file mode 100644 index 0000000..6430d72 --- /dev/null +++ b/vignettes/parallel_radEmu.Rmd @@ -0,0 +1,186 @@ +--- +title: "Parallelizing computation for score tests with radEmu" +author: "Sarah Teichman" +date: "`r Sys.Date()`" +output: rmarkdown::html_vignette +vignette: > + %\VignetteIndexEntry{Parallelizing computation for score tests with radEmu} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + +```{r, include = FALSE} +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>" +) +``` + +First, we will install `radEmu`, if we haven't already. + +```{r, eval = FALSE} +# if (!require("remotes", quietly = TRUE)) +# install.packages("remotes") +# +# remotes::install_github("statdivlab/radEmu") +``` + +Next, we can load `radEmu` as well as the `tidyverse` package suite. + +```{r setup, message = FALSE} +library(magrittr) +library(dplyr) +library(ggplot2) +library(stringr) +library(radEmu) +``` + +Finally, we will load the package `parallel`. + +```{r, message = FALSE} +library(parallel) +``` + +## Introduction + +In this vignette we will introduce parallel computing in order to do more efficient computation for score tests. Our recommendation for testing statistical hypotheses with small to moderate sample sizes with `radEmu` is to run a robust score test. While this test performs well (we like that it controls Type I error rate even in small samples!), it takes some time to run, because we need to fit the model under each null hypothesis. For differential abundance analysis, we often want to run a hypothesis test for each category (taxon, gene, etc) that we care about, so this adds up quickly. +In order to improve computational efficiency, we can run these score tests in parallel using the `parallel` R package. This will let us take advantage of having additional cores on personal computers or computing clusters. Note that we will be using the `mclapply` function from the `parallel` package, which works on Mac and Linux machines, but does not on Windows. If you are using a Windows machine and would like a vignette about parallel computing on Windows, please let us know by opening an issue. + +We recommend that before working through this vignette, you start the an introduction to `radEmu` in "intro_radEmu.Rmd." + +## Setting up our radEmu model and running a single score test + +We'll use the same Wirbel et al. data as in the introduction vignette. Recall that the [dataset published by Wirbel et al. (2019)](https://www.nature.com/articles/s41591-019-0406-6) is from a meta-analysis of case-controls comparing participants with and without colorectal cancer. + +```{r} +# load in sample data +data("wirbel_sample") +# set group to be a factor with levels CTR for control and CRC for cancer +wirbel_sample$Group <- factor(wirbel_sample$Group, levels = c("CTR","CRC")) +# load in abundance data +data("wirbel_otu") +# save mOTU names +mOTU_names <- colnames(wirbel_otu) +# consider taxa in the following genera +chosen_genera <- c("Eubacterium", "Faecalibacterium", "Fusobacterium", "Porphyromonas") +# get taxonomy information from mOTU names +mOTU_name_df <- data.frame(name = mOTU_names) %>% + mutate(base_name = stringr::str_remove(mOTU_names, "unknown ") %>% + stringr::str_remove("uncultured ")) %>% + mutate(genus_name = stringr::word(base_name, 1)) +# restrict to names in chosen genera +restricted_mOTU_names <- mOTU_name_df %>% + filter(genus_name %in% chosen_genera) %>% + pull(name) +# pull out observations from a chinese study within the meta-analysis +ch_study_obs <- which(wirbel_sample$Country %in% c("CHI")) +# make count matrix for chosen samples and genera +small_Y <- wirbel_otu[ch_study_obs, restricted_mOTU_names] +# check for samples with only zero counts +sum(rowSums(small_Y) == 0) # no samples have a count sum of 0 +# check for genera with only zero counts +sum(colSums(small_Y) == 0) # one category has a count sum of 0 +# remove the one genus with only zero counts +category_to_rm <- which(colSums(small_Y) == 0) +small_Y <- small_Y[, -category_to_rm] +``` + +Now that we've processed our data, we can fit the `radEmu` model. Here we just want to get estimates for our parameters and their standard errors, but we will avoid running score tests by setting `run_score_tests = FALSE`. + +```{r} +ch_fit <- emuFit(formula = ~ Group, + data = wirbel_sample[ch_study_obs, ], + Y = small_Y, + run_score_tests = FALSE) +``` + +In the introduction vignette we found that a meta-mOTU "unknown Eubacterium [meta_mOTU_v2_7116]" assigned to Eubacteria has a much higher ratio of abundance (comparing CRC group to control) than is typical across the mOTUs we included in this analysis, based on the parameter estimates in `ch_fit`. We can run a robust score test to test whether the differential abundance of this mOTU between cases and controls is significantly different from the differential abundance of a typical mOTU in our analysis. + +In order to run a single robust score test, we will set `run_score_tests = TRUE` and include the argument `test_kj`. Instead of re-estimating parameters in our model, we will provide `ch_fit` to the argument `fitted_model` and set `refit = FALSE`. + +```{r} +mOTU_to_test <- which(str_detect(restricted_mOTU_names, "7116")) +ch_fit$B %>% rownames +covariate_to_test <- which("GroupCRC" == ch_fit$B %>% rownames) +robust_score <- emuFit(formula = ~ Group, + data = wirbel_sample[ch_study_obs, ], + fitted_model = ch_fit, + refit = FALSE, + test_kj = data.frame(k = covariate_to_test, + j = mOTU_to_test), + Y = small_Y) +robust_score$coef$pval[mOTU_to_test] +``` + +Now, we can see that it took a little while to run our robust score test. If we investigate the coefficient table in our `robust_score` output, we can see that we have a p-value of `R robust_score$coef$pval[mOTU_to_test]` from our test. + +## Running robust score tests in parallel + +Now, let's run some tests in parallel. We will be parallelizing our code over $j$ in the argument `test_kj`. We will assume that you have one covariate that you want to test, corresponding with a specific column of your design matrix $k$. However, if you want to run tests for multiple columns of your design matrix, then you can parallelize over pairs of $k$ and $j$ in the argument $test_kj$. + +Let's say that I want to run score tests for the first five mOTUs in our dataset. First, I need to check the cores on my computer to see a reasonable amount of cores to parallelize over. I tend to use one fewer core than the number of cores that I have available. + +```{r} +ncores <- parallel::detectCores() - 1 +ncores +``` + +Next, I will write a function that will be called in parallel. This function will fit the model under the null and calculate the robust score test statistics. Note that the output of this function is an `emuFit` object. + +```{r} +emuTest <- function(category) { + score_res <- emuFit(formula = ~ Group, + data = wirbel_sample[ch_study_obs, ], + fitted_model = ch_fit, + refit = FALSE, + test_kj = data.frame(k = covariate_to_test, + j = category), + Y = small_Y) + return(score_res) +} +``` + +Now, we can run our score tests in parallel. We'll just do the first five. It may take a minute or so, depending on your machine. + +```{r} +if (.Platform$OS.type != "windows") { + # run if we are on a Mac or Linux machine + score_res <- mclapply(1:5, + emuTest, + mc.cores = ncores) +} else { + # don't run if we are on a Windows machine + score_res <- NULL +} +``` + +Now, we can see that this barely took more time than running a single score test, because we were able to parallelize over cores (on my laptop, I'm using more than five cores, so I can run all five tests at the same time). Each p-value can be pulled out of the list as follows: + + +```{r} +if (!is.null(score_res)) { + c(score_res[[1]]$coef$pval[1], ## robust score test p-value for the first taxon + score_res[[2]]$coef$pval[2]) ## robust score test p-value for the second taxon +} +``` + +To help organise this information, we can make a coefficient matrix that combines the information from each component in our list: + +```{r} +if (!is.null(score_res)) { + full_score <- sapply(1:length(score_res), + function(x) score_res[[x]]$coef$score_stat[x]) + full_pval <- sapply(1:length(score_res), + function(x) score_res[[x]]$coef$pval[x]) + full_coef <- ch_fit$coef %>% + dplyr::select(-score_stat, -pval) %>% + filter(category_num %in% 1:5) %>% + mutate(score_stat = full_score, + pval = full_pval) + full_coef +} +``` + +The column containing our p-values is called `pval`. + +Happy testing! \ No newline at end of file diff --git a/vignettes/radEmu_clustered_data.Rmd b/vignettes/radEmu_clustered_data.Rmd new file mode 100644 index 0000000..70ffa71 --- /dev/null +++ b/vignettes/radEmu_clustered_data.Rmd @@ -0,0 +1,132 @@ +--- +title: "Using radEmu on clustered data" +author: "Sarah Teichman and Amy Willis" +date: "`r Sys.Date()`" +output: rmarkdown::html_vignette +vignette: > + %\VignetteIndexEntry{Using radEmu on clustered data} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + +```{r, include = FALSE} +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>" +) +``` + +First, we will install `radEmu`, if we haven't already. + +```{r, eval = FALSE} +# if (!require("remotes", quietly = TRUE)) +# install.packages("remotes") +# +# remotes::install_github("statdivlab/radEmu") +``` + +Next, we can load `radEmu` as well as the `tidyverse` package suite. + +```{r setup, message = FALSE} +library(magrittr) +library(dplyr) +library(ggplot2) +library(stringr) +library(radEmu) +``` + +## Introduction + +In this vignette we will explore the use of `radEmu` with clustered data. Data with cluster dependence is a common phenomenon in microbiome studies. This type of dependence can come from experimental factors such as shared cages or tanks for study animals. It can also arise from repeated measurements in longitudinal studies. + +Sometimes people deal with clustered data via random effects. You might have seen the syntax `+ (1 | cluster)`. We (the `radEmu` developers) don't want to make strong assumptions (such as the random effects being normally distributed), so we handle cluster dependence using a GEE framework. We like this approach because it's robust to many forms of model misspecification, and for this reason, we find it superior to random effects. + +When dependence is not accounted for, statistical inference in `radEmu` (and most statistical methods!) assume that all samples are independent. If you have cluster dependence but assume independence, you'll have anti-conservative inference (i.e., p-values that are smaller than they should be). **Therefore, we strongly recommend adjusting for cluster dependence if it arose in your sample collection!** + +Note that cluster dependence won't change your estimates ($\hat{\beta}_j$'s), but it will (most likely) change your p-values. + +Luckily, we have tools to account for cluster dependence implemented in `radEmu`! This argument is only implemented from `radEmu` v1.2.0.0 forward, so if you are having trouble using the `cluster` argument, check that you reinstalled `radEmu` recently. + +TLDR; the basic syntax is as follows. `my_clusters` is a vector of length $n$ (your number of samples), with observations from the same cluster having the same value in `my_clusters` (e.g. `my_clusters = c(1, 1, 2, 2, 3, 4)`). + +``` +emuFit(formula = ~ covariate, + data = my_data_frame, + Y = my_microbial_abundances, + cluster = my_clusters) +``` + +Fun fact! We implemented this functionality because of user requests. Therefore, if there's something that you'd like to see that you don't see, [let us know](https://github.com/statdivlab/radEmu/issues) and we'll see what we can do! + +## Generating data with cluster dependence + +To start, let's generate a toy example (10 categories, 60 samples) in which there is cluster dependence within our data. The way we simulate data isn't important; it's just an illustration. **radEmu can handle more taxa and samples,** we just did this so that the vignette builds quickly. + +```{r} +J <- 10; n <- 60 +# generate design matrix +set.seed(10) +X <- cbind(1, rnorm(n)) +cov_dat <- data.frame(cov = X[, 2]) +# cluster membership +cluster <- rep(1:4, each = n/4) +cluster_named <- paste("cage", cluster, sep = "") +cov_dat$cluster <- cluster +# intercepts for each category +b0 <- rnorm(J) +# coefficients for X1 for each category +b1 <- seq(1, 5, length.out = J) +# mean center the coefficients +b1 <- b1 - mean(b1) +# set the coefficient for the 3rd category to 4 (why not!?) +# Note that because of the constraint, we're only able to estimate +# b1 - mean(b1), which is ~3.9. +b1[3] <- 4 +# generate B coefficient matrix +b <- rbind(b0, b1) + +# simulate data according to a zero-inflated negative binomial distribution +# the mean model used to simulate this data takes into account the cluster membership +Y <- radEmu:::simulate_data(n = n, + J = J, + b0 = b0, + b1 = b1, + distn = "ZINB", + zinb_size = 10, + zinb_zero_prop = 0.3, + mean_z = 5, + X = X, + cluster = cluster) +``` + +Let's just pause to look at the elements of `cluster`: + +```{r} +table(cluster_named) +``` + +So all of the observations from the first cage (or person, tank, whatever...) have `cage1` in their corresponding `cluster_named` variable. + +Let's fit a model to this data! We know that the log-fold difference between in the abundance of category 3 when comparing samples that differ by 1 unit in $X$ is `r round((b1 - mean(b1))[3], 1)`, so if we have good power, we will reject the null that $\beta_{X_1, 3} = 0$. We fit a model including cluster dependence as follows: + +```{r} +ef_cluster <- emuFit(formula = ~ cov, + data = cov_dat, + Y = Y, + cluster = cluster_named, + test_kj = data.frame(k = 2, j = 3)) +``` + +You can check out the full object, but our estimate is `r round(ef_cluster$coef$estimate[3], 1)` and a p-value for testing that this parameter equals zero is `r round(ef_cluster$coef$pval[3], 3)`. Not too shabby, especially considering that about half of our observations are zero, and we have a lot of noise in our data (arising from a negative binomial simulation scheme). + +Let's also compare that to a situation where we mistakenly ignore clustering. In this case, we expect to have a smaller p-value, because we are saying that we have more independent observations. + +```{r} +ef_no_cluster <- emuFit(formula = ~ cov, + data = cov_dat, + Y = Y, + test_kj = data.frame(k = 2, j = 3)) +``` + +When we ignore clustering, we get an estimate of the log fold difference in category 3 across values of our covariate of `r round(ef_no_cluster$coef$estimate[3], 1)` and a p-value of `r round(ef_no_cluster$coef$pval[3], 3)`. So we can see that our estimates are the same whether or not we account for cluster, but our p-values are different (because we are pretending that we have more evidence in the absence of clustering). +