From 76160b5172ebcce05dc6b1990f50418a8b813a67 Mon Sep 17 00:00:00 2001 From: erik-ringen Date: Fri, 11 Oct 2024 17:20:58 +0200 Subject: [PATCH 1/2] allow non-intervention values in calculate_theta() --- R/coev_calculate_theta.R | 208 ++++++++++++--------- man/coev_calculate_theta.Rd | 31 ++- tests/testthat/test-coev_calculate_theta.R | 11 ++ 3 files changed, 151 insertions(+), 99 deletions(-) diff --git a/R/coev_calculate_theta.R b/R/coev_calculate_theta.R index d9c9c03..0e79013 100644 --- a/R/coev_calculate_theta.R +++ b/R/coev_calculate_theta.R @@ -5,11 +5,10 @@ #' object. #' #' @param object An object of class \code{coevfit} -#' @param intervention_values A named list of variables and associated -#' intervention values for calculating equilibrium trait values. All -#' coevolving variables must be declared separately in the named list without +#' @param intervention_values Either \code{NULL} (the default) or a named list of variables and associated +#' intervention values for calculating equilibrium trait values. If \code{NULL}, calculates the equiilbirium states when all parameters are free to vary. Otherwise, all coevolving variables must be declared separately in the named list without #' repetition. If the intervention value for a particular variable is set to -#' NA, this variable is treated as a free variable. Otherwise, if the +#' NA, this variable is treated as a free variable. If the #' intervention value for a particular variable is specified, the variable is #' held constant at this trait value in the calculation. At least one variable #' must be declared as a free variable and at least one variable must be held @@ -20,11 +19,22 @@ #' @author Scott Claessens \email{scott.claessens@@gmail.com}, Erik Ringen #' \email{erikjacob.ringen@@uzh.ch} #' -#' @details The equilibrium trait value for a free trait \eqn{\eta_i} is +#' @details The equilibrium trait values for freely evolving traits \eqn{\mathbf{\eta}} are #' calculated using the following formula: -#' \deqn{\theta_{\eta_i} = \frac{-(\sum_{j \neq i}\textbf{A}[i,j]\eta_j) + -#' \textbf{b}_i}{\textbf{A}[i,i]}} -#' +#' \deqn{\mathbf{\theta} = -\mathbf{A}^{-1}\mathbf{b}} +#' +#' If we hold some variables constant at some value(s) (denoted \eqn{\eta_{h}}) and let others evolve freely (denoted \eqn{\eta_{f}}), we can partition the parameters as follows: +#' +#' \deqn{\mathbf{A}_{ff}: \text{selection coefficients to/from the free variables}} +#' \deqn{\mathbf{A}_{fh}: \text{selection coefficients to the free variables from the held variables}} +#' \deqn{\mathbf{b}_{f}: \text{continuous time intercepts for the free variables}} +#' +#' We can then calculate the equilibrium values for the free variables: +#' \deqn{\boldsymbol{\theta_f} = -\mathbf{A}_{ff}^{-1} \left( \mathbf{A}_{fh} \mathbf{\eta}_h + \mathbf{b}_f \right)} +#' +#' With the overall equilibrium vector being a mix of the free equilibria and the held values: +#' \deqn{\boldsymbol{\theta | \mathbf{\eta}_h} = \begin{bmatrix} \boldsymbol{\theta}_f \\\mathbf{\eta}_h\end{bmatrix}} +#' #' @references #' Ringen, E., Martin, J. S., & Jaeggi, A. (2021). Novel phylogenetic methods #' reveal that resource-use intensification drives the evolution of "complex" @@ -54,6 +64,11 @@ #' parallel_chains = 4, #' seed = 1 #' ) +#' +#' # calculate theta with no interventions +#' coev_calculate_theta( +#' object = fit +#' ) #' #' # calculate theta given intervention values #' coev_calculate_theta( @@ -66,7 +81,7 @@ #' } #' #' @export -coev_calculate_theta <- function(object, intervention_values) { +coev_calculate_theta <- function(object, intervention_values = NULL) { # stop if object is not of class coevfit if (!methods::is(object, "coevfit")) { stop2( @@ -76,95 +91,106 @@ coev_calculate_theta <- function(object, intervention_values) { ) ) } - # stop if intervention_values argument is not a named list - if (!is.list(intervention_values) | is.null(names(intervention_values))) { - stop2("Argument 'intervention_values' is not a named list.") - } - # stop if intervention_list contains variables not included in the model - if (!all(names(intervention_values) %in% names(object$variables))) { - stop2( - paste0( - "At least one variable in 'intervention_values' is not included in ", - "the fitted model." - ) - ) - } - # stop if any coevolving variables are not listed in intervention_values - if (any(!(names(object$variables) %in% names(intervention_values)))) { - stop2( - paste0( - "All coevolving variables must be included in ", - "argument 'intervention_values'." - ) - ) - } - # stop if repetition in intervention_values - if (any(duplicated(names(intervention_values)))) { - stop2("Argument 'intervention_values' contains duplicated variable names.") - } - # stop if any values in intervention_list are not of length one - if (any(unlist(lapply(intervention_values, length)) != 1)) { - stop2("Values in 'intervention_values' must each be of length one.") - } - # stop if any values in intervention_list are not NA or numeric - if (any(unlist(lapply(intervention_values, - function(x) !is.na(x) & !is.numeric(x))))) { - stop2("Values in 'intervention_values' must each be NA or numeric.") - } - # stop if all variables are held constant in intervention_list - if (mean(is.na(intervention_values)) == 0) { - stop2( - paste0( - "Argument 'intervention_values' must have at least one NA value ", - "declaring a free variable. If all variables are held constant, the ", - "system is already at equilibrium and there is nothing to compute." - ) - ) - } - # stop if all variables in intervention_values are free - if (mean(is.na(intervention_values)) == 1) { - stop2( - paste0( - "Argument 'intervention_values' must have at least one variable ", - "held constant (i.e., all values are NA)." - ) - ) - } - # construct intervention values x_hat - x_hat <- unlist(intervention_values[names(object$variables)]) # extract posterior draws post <- posterior::as_draws_rvars(object$fit$draws(variables = c("A", "b"))) A <- posterior::draws_of(post$A) b <- posterior::draws_of(post$b) - # partition the matrices - held_indices <- which(!is.na(x_hat)) - free_indices <- which(is.na(x_hat)) - A_free_free <- A[,free_indices, free_indices, drop = FALSE] - A_free_held <- A[,free_indices, held_indices, drop = FALSE] - b_free <- b[,free_indices, drop = FALSE] # construct theta matrix - theta <- matrix(NA, nrow = posterior::ndraws(post), ncol = length(x_hat)) - for (i in 1:posterior::ndraws(post)) { - # compute the equilibrium for the free variables - if (length(held_indices) == 1) { - equilibrium_free <- - -solve( A_free_free[i,,] ) %*% - (b_free[i,] + as.matrix(A_free_held[i,,]) %*% x_hat[!is.na(x_hat)]) - } else { - equilibrium_free <- - -solve( A_free_free[i,,] ) %*% - (b_free[i,] + A_free_held[i,,] %*% x_hat[!is.na(x_hat)]) + theta <- matrix(NA, nrow = posterior::ndraws(post), ncol = length(object$variables)) + + # Non-intervention + if (is.null(intervention_values)) { + for (i in 1:posterior::ndraws(post)) { + theta[i,] = -solve(A[i,,]) %*% b[i,] + } + } + # Intervention (at least one value held) + else{ + # stop if intervention_values argument is not a named list + if (!is.list(intervention_values) | is.null(names(intervention_values))) { + stop2("Argument 'intervention_values' is not a named list.") + } + # stop if intervention_list contains variables not included in the model + if (!all(names(intervention_values) %in% names(object$variables))) { + stop2( + paste0( + "At least one variable in 'intervention_values' is not included in ", + "the fitted model." + ) + ) + } + # stop if any coevolving variables are not listed in intervention_values + if (any(!(names(object$variables) %in% names(intervention_values)))) { + stop2( + paste0( + "All coevolving variables must be included in ", + "argument 'intervention_values'." + ) + ) + } + # stop if repetition in intervention_values + if (any(duplicated(names(intervention_values)))) { + stop2("Argument 'intervention_values' contains duplicated variable names.") + } + # stop if any values in intervention_list are not of length one + if (any(unlist(lapply(intervention_values, length)) != 1)) { + stop2("Values in 'intervention_values' must each be of length one.") + } + # stop if any values in intervention_list are not NA or numeric + if (any(unlist(lapply(intervention_values, + function(x) !is.na(x) & !is.numeric(x))))) { + stop2("Values in 'intervention_values' must each be NA or numeric.") + } + # stop if all variables are held constant in intervention_list + if (mean(is.na(intervention_values)) == 0) { + stop2( + paste0( + "Argument 'intervention_values' must have at least one NA value ", + "declaring a free variable. If all variables are held constant, the ", + "system is already at equilibrium and there is nothing to compute." + ) + ) + } + # stop if all variables in intervention_values are free + if (mean(is.na(intervention_values)) == 1) { + stop2( + paste0( + "Argument 'intervention_values' must have at least one variable ", + "held constant (i.e., all values are NA)." + ) + ) + } + # construct intervention values x_hat + x_hat <- unlist(intervention_values[names(object$variables)]) + # partition the matrices + held_indices <- which(!is.na(x_hat)) + free_indices <- which(is.na(x_hat)) + A_free_free <- A[,free_indices, free_indices, drop = FALSE] + A_free_held <- A[,free_indices, held_indices, drop = FALSE] + b_free <- b[,free_indices, drop = FALSE] + + for (i in 1:posterior::ndraws(post)) { + # compute the equilibrium for the free variables + if (length(held_indices) == 1) { + equilibrium_free <- + -solve( A_free_free[i,,] ) %*% + (b_free[i,] + as.matrix(A_free_held[i,,]) %*% x_hat[!is.na(x_hat)]) + } else { + equilibrium_free <- + -solve( A_free_free[i,,] ) %*% + (b_free[i,] + A_free_held[i,,] %*% x_hat[!is.na(x_hat)]) + } + # initialize the equilibrium vector + equilibrium <- rep(NA, length(x_hat)) + # fill in the computed equilibrium values for free variables + equilibrium[free_indices] <- equilibrium_free + # fill in the constant values for held variables + equilibrium[held_indices] <- x_hat[!is.na(x_hat)] + # add to theta matrix + theta[i,] = equilibrium } - # initialize the equilibrium vector - equilibrium <- rep(NA, length(x_hat)) - # fill in the computed equilibrium values for free variables - equilibrium[free_indices] <- equilibrium_free - # fill in the constant values for held variables - equilibrium[held_indices] <- x_hat[!is.na(x_hat)] - # add to theta matrix - theta[i,] = equilibrium } # add column names to theta matrix colnames(theta) <- names(object$variables) return(theta) -} + } diff --git a/man/coev_calculate_theta.Rd b/man/coev_calculate_theta.Rd index 7272eab..1d81edd 100644 --- a/man/coev_calculate_theta.Rd +++ b/man/coev_calculate_theta.Rd @@ -4,16 +4,15 @@ \alias{coev_calculate_theta} \title{Calculate equilibrium trait values (theta) for a fitted \code{coevfit} object} \usage{ -coev_calculate_theta(object, intervention_values) +coev_calculate_theta(object, intervention_values = NULL) } \arguments{ \item{object}{An object of class \code{coevfit}} -\item{intervention_values}{A named list of variables and associated -intervention values for calculating equilibrium trait values. All -coevolving variables must be declared separately in the named list without +\item{intervention_values}{Either \code{NULL} (the default) or a named list of variables and associated +intervention values for calculating equilibrium trait values. If \code{NULL}, calculates the equiilbirium states when all parameters are free to vary. Otherwise, all coevolving variables must be declared separately in the named list without repetition. If the intervention value for a particular variable is set to -NA, this variable is treated as a free variable. Otherwise, if the +NA, this variable is treated as a free variable. If the intervention value for a particular variable is specified, the variable is held constant at this trait value in the calculation. At least one variable must be declared as a free variable and at least one variable must be held @@ -28,10 +27,21 @@ a set of intervention values for other traits from a fitted \code{coevfit} object. } \details{ -The equilibrium trait value for a free trait \eqn{\eta_i} is +The equilibrium trait values for freely evolving traits \eqn{\mathbf{\eta}} are calculated using the following formula: -\deqn{\theta_{\eta_i} = \frac{-(\sum_{j \neq i}\textbf{A}[i,j]\eta_j) + - \textbf{b}_i}{\textbf{A}[i,i]}} +\deqn{\mathbf{\theta} = -\mathbf{A}^{-1}\mathbf{b}} + +If we hold some variables constant at some value(s) (denoted \eqn{\eta_{h}}) and let others evolve freely (denoted \eqn{\eta_{f}}), we can partition the parameters as follows: + +\deqn{\mathbf{A}_{ff}: \text{selection coefficients to/from the free variables}} +\deqn{\mathbf{A}_{fh}: \text{selection coefficients to the free variables from the held variables}} +\deqn{\mathbf{b}_{f}: \text{continuous time intercepts for the free variables}} + +We can then calculate the equilibrium values for the free variables: +\deqn{\boldsymbol{\theta_f} = -\mathbf{A}_{ff}^{-1} \left( \mathbf{A}_{fh} \mathbf{\eta}_h + \mathbf{b}_f \right)} + +With the overall equilibrium vector being a mix of the free equilibria and the held values: +\deqn{\boldsymbol{\theta | \mathbf{\eta}_h} = \begin{bmatrix} \boldsymbol{\theta}_f \\\mathbf{\eta}_h\end{bmatrix}} } \examples{ \dontrun{ @@ -50,6 +60,11 @@ fit <- coev_fit( seed = 1 ) +# calculate theta with no interventions +coev_calculate_theta( + object = fit + ) + # calculate theta given intervention values coev_calculate_theta( object = fit, diff --git a/tests/testthat/test-coev_calculate_theta.R b/tests/testthat/test-coev_calculate_theta.R index a53d41a..a51a73c 100644 --- a/tests/testthat/test-coev_calculate_theta.R +++ b/tests/testthat/test-coev_calculate_theta.R @@ -94,6 +94,8 @@ test_that("coev_calculate_theta() produces expected errors and output", { theta7 <- coev_calculate_theta(m7, list(w = NA, x = 0)) theta8 <- coev_calculate_theta(m8, list(x = NA, y = 0)) theta9 <- coev_calculate_theta(m9, list(x = NA, y = 0)) + theta1_null <- coev_calculate_theta(m1, NULL) + theta9_null <- coev_calculate_theta(m9, NULL) expect_no_error(theta1) expect_no_error(theta2) expect_no_error(theta3) @@ -103,6 +105,8 @@ test_that("coev_calculate_theta() produces expected errors and output", { expect_no_error(theta7) expect_no_error(theta8) expect_no_error(theta9) + expect_no_error(theta1_null) + expect_no_error(theta9_null) # output should be matrix of posterior draws expect_true(methods::is(theta1, "matrix")) expect_true(methods::is(theta2, "matrix")) @@ -113,6 +117,8 @@ test_that("coev_calculate_theta() produces expected errors and output", { expect_true(methods::is(theta7, "matrix")) expect_true(methods::is(theta8, "matrix")) expect_true(methods::is(theta9, "matrix")) + expect_true(methods::is(theta1_null, "matrix")) + expect_true(methods::is(theta9_null, "matrix")) # output column names should be equal to variable names expect_true(identical(colnames(theta1), names(m1$variables))) expect_true(identical(colnames(theta2), names(m2$variables))) @@ -123,6 +129,8 @@ test_that("coev_calculate_theta() produces expected errors and output", { expect_true(identical(colnames(theta7), names(m7$variables))) expect_true(identical(colnames(theta8), names(m8$variables))) expect_true(identical(colnames(theta9), names(m9$variables))) + expect_true(identical(colnames(theta9_null), names(m9$variables))) + expect_true(identical(colnames(theta1_null), names(m1$variables))) # variables should be correctly held constant in output expect_true(all(theta1[,c("v","w","x","y")] == 0)) expect_true(all(theta2[,"x"] == 0)) @@ -133,4 +141,7 @@ test_that("coev_calculate_theta() produces expected errors and output", { expect_true(all(theta7[,"x"] == 0)) expect_true(all(theta8[,"y"] == 0)) expect_true(all(theta9[,"y"] == 0)) + # free variables should not be constant in output + expect_false(all(theta1_null[,"y"] == 0)) + expect_false(all(theta9_null[,"y"] == 0)) }) From a29208e90b777416917f1ffe708d56e2772d0e99 Mon Sep 17 00:00:00 2001 From: Scott Claessens Date: Fri, 11 Oct 2024 16:51:35 +0100 Subject: [PATCH 2/2] Minor tidying --- R/coev_calculate_theta.R | 74 +++++++++++++--------- man/coev_calculate_theta.Rd | 51 +++++++++------ tests/testthat/test-coev_calculate_theta.R | 34 +++++++++- 3 files changed, 109 insertions(+), 50 deletions(-) diff --git a/R/coev_calculate_theta.R b/R/coev_calculate_theta.R index 0e79013..e356deb 100644 --- a/R/coev_calculate_theta.R +++ b/R/coev_calculate_theta.R @@ -5,36 +5,50 @@ #' object. #' #' @param object An object of class \code{coevfit} -#' @param intervention_values Either \code{NULL} (the default) or a named list of variables and associated -#' intervention values for calculating equilibrium trait values. If \code{NULL}, calculates the equiilbirium states when all parameters are free to vary. Otherwise, all coevolving variables must be declared separately in the named list without -#' repetition. If the intervention value for a particular variable is set to -#' NA, this variable is treated as a free variable. If the -#' intervention value for a particular variable is specified, the variable is -#' held constant at this trait value in the calculation. At least one variable -#' must be declared as a free variable and at least one variable must be held -#' constant (e.g., \code{list(var1 = NA, var2 = 0)}). +#' @param intervention_values Either \code{NULL} (the default) or a named list +#' of variables and associated intervention values for calculating equilibrium +#' trait values. If \code{NULL}, calculates the equilibrium states when all +#' parameters are free to vary. Otherwise, all coevolving variables must be +#' declared separately in the named list without repetition. If the +#' intervention value for a particular variable is set to NA, this variable is +#' treated as a free variable. If the intervention value for a particular +#' variable is specified, the variable is held constant at this trait value in +#' the calculation. At least one variable must be declared as a free variable +#' and at least one variable must be held constant (e.g., +#' \code{list(var1 = NA, var2 = 0)}). #' #' @return Posterior samples in matrix format #' #' @author Scott Claessens \email{scott.claessens@@gmail.com}, Erik Ringen #' \email{erikjacob.ringen@@uzh.ch} #' -#' @details The equilibrium trait values for freely evolving traits \eqn{\mathbf{\eta}} are -#' calculated using the following formula: +#' @details The equilibrium trait values for freely evolving traits +#' \eqn{\mathbf{\eta}} are calculated using the following formula: +#' #' \deqn{\mathbf{\theta} = -\mathbf{A}^{-1}\mathbf{b}} -#' -#' If we hold some variables constant at some value(s) (denoted \eqn{\eta_{h}}) and let others evolve freely (denoted \eqn{\eta_{f}}), we can partition the parameters as follows: -#' -#' \deqn{\mathbf{A}_{ff}: \text{selection coefficients to/from the free variables}} -#' \deqn{\mathbf{A}_{fh}: \text{selection coefficients to the free variables from the held variables}} -#' \deqn{\mathbf{b}_{f}: \text{continuous time intercepts for the free variables}} -#' -#' We can then calculate the equilibrium values for the free variables: -#' \deqn{\boldsymbol{\theta_f} = -\mathbf{A}_{ff}^{-1} \left( \mathbf{A}_{fh} \mathbf{\eta}_h + \mathbf{b}_f \right)} -#' -#' With the overall equilibrium vector being a mix of the free equilibria and the held values: -#' \deqn{\boldsymbol{\theta | \mathbf{\eta}_h} = \begin{bmatrix} \boldsymbol{\theta}_f \\\mathbf{\eta}_h\end{bmatrix}} -#' +#' +#' If we hold some variables constant at some value(s) (denoted +#' \eqn{\eta_{h}}) and let others evolve freely (denoted \eqn{\eta_{f}}), we +#' can partition the parameters as follows: +#' +#' - \eqn{\mathbf{A}_{ff}}: selection coefficients to/from the free +#' variables +#' - \eqn{\mathbf{A}_{fh}}: selection coefficients to the free variables +#' from the held variables +#' - \eqn{\mathbf{b}_{f}}: continuous time intercepts for the free +#' variables +#' +#' We can then calculate the equilibrium values for the free variables: +#' +#' \deqn{\boldsymbol{\theta_f} = -\mathbf{A}_{ff}^{-1} \left( \mathbf{A}_{fh} +#' \mathbf{\eta}_h + \mathbf{b}_f \right)} +#' +#' With the overall equilibrium vector being a mix of the free equilibria and +#' the held values: +#' +#' \deqn{\boldsymbol{\theta | \mathbf{\eta}_h} = \begin{bmatrix} +#' \boldsymbol{\theta}_f \\\mathbf{\eta}_h\end{bmatrix}} +#' #' @references #' Ringen, E., Martin, J. S., & Jaeggi, A. (2021). Novel phylogenetic methods #' reveal that resource-use intensification drives the evolution of "complex" @@ -64,7 +78,7 @@ #' parallel_chains = 4, #' seed = 1 #' ) -#' +#' #' # calculate theta with no interventions #' coev_calculate_theta( #' object = fit @@ -96,15 +110,15 @@ coev_calculate_theta <- function(object, intervention_values = NULL) { A <- posterior::draws_of(post$A) b <- posterior::draws_of(post$b) # construct theta matrix - theta <- matrix(NA, nrow = posterior::ndraws(post), ncol = length(object$variables)) - - # Non-intervention + theta <- matrix(NA, nrow = posterior::ndraws(post), + ncol = length(object$variables)) + # non-intervention if (is.null(intervention_values)) { for (i in 1:posterior::ndraws(post)) { theta[i,] = -solve(A[i,,]) %*% b[i,] } } - # Intervention (at least one value held) + # intervention (at least one value held) else{ # stop if intervention_values argument is not a named list if (!is.list(intervention_values) | is.null(names(intervention_values))) { @@ -130,7 +144,9 @@ coev_calculate_theta <- function(object, intervention_values = NULL) { } # stop if repetition in intervention_values if (any(duplicated(names(intervention_values)))) { - stop2("Argument 'intervention_values' contains duplicated variable names.") + stop2( + "Argument 'intervention_values' contains duplicated variable names." + ) } # stop if any values in intervention_list are not of length one if (any(unlist(lapply(intervention_values, length)) != 1)) { diff --git a/man/coev_calculate_theta.Rd b/man/coev_calculate_theta.Rd index 1d81edd..25389c9 100644 --- a/man/coev_calculate_theta.Rd +++ b/man/coev_calculate_theta.Rd @@ -9,14 +9,17 @@ coev_calculate_theta(object, intervention_values = NULL) \arguments{ \item{object}{An object of class \code{coevfit}} -\item{intervention_values}{Either \code{NULL} (the default) or a named list of variables and associated -intervention values for calculating equilibrium trait values. If \code{NULL}, calculates the equiilbirium states when all parameters are free to vary. Otherwise, all coevolving variables must be declared separately in the named list without -repetition. If the intervention value for a particular variable is set to -NA, this variable is treated as a free variable. If the -intervention value for a particular variable is specified, the variable is -held constant at this trait value in the calculation. At least one variable -must be declared as a free variable and at least one variable must be held -constant (e.g., \code{list(var1 = NA, var2 = 0)}).} +\item{intervention_values}{Either \code{NULL} (the default) or a named list +of variables and associated intervention values for calculating equilibrium +trait values. If \code{NULL}, calculates the equilibrium states when all +parameters are free to vary. Otherwise, all coevolving variables must be +declared separately in the named list without repetition. If the +intervention value for a particular variable is set to NA, this variable is +treated as a free variable. If the intervention value for a particular +variable is specified, the variable is held constant at this trait value in +the calculation. At least one variable must be declared as a free variable +and at least one variable must be held constant (e.g., +\code{list(var1 = NA, var2 = 0)}).} } \value{ Posterior samples in matrix format @@ -27,21 +30,33 @@ a set of intervention values for other traits from a fitted \code{coevfit} object. } \details{ -The equilibrium trait values for freely evolving traits \eqn{\mathbf{\eta}} are -calculated using the following formula: -\deqn{\mathbf{\theta} = -\mathbf{A}^{-1}\mathbf{b}} +The equilibrium trait values for freely evolving traits +\eqn{\mathbf{\eta}} are calculated using the following formula: -If we hold some variables constant at some value(s) (denoted \eqn{\eta_{h}}) and let others evolve freely (denoted \eqn{\eta_{f}}), we can partition the parameters as follows: +\deqn{\mathbf{\theta} = -\mathbf{A}^{-1}\mathbf{b}} -\deqn{\mathbf{A}_{ff}: \text{selection coefficients to/from the free variables}} -\deqn{\mathbf{A}_{fh}: \text{selection coefficients to the free variables from the held variables}} -\deqn{\mathbf{b}_{f}: \text{continuous time intercepts for the free variables}} +If we hold some variables constant at some value(s) (denoted +\eqn{\eta_{h}}) and let others evolve freely (denoted \eqn{\eta_{f}}), we +can partition the parameters as follows: +\itemize{ +\item \eqn{\mathbf{A}_{ff}}: selection coefficients to/from the free +variables +\item \eqn{\mathbf{A}_{fh}}: selection coefficients to the free variables +from the held variables +\item \eqn{\mathbf{b}_{f}}: continuous time intercepts for the free +variables +} We can then calculate the equilibrium values for the free variables: -\deqn{\boldsymbol{\theta_f} = -\mathbf{A}_{ff}^{-1} \left( \mathbf{A}_{fh} \mathbf{\eta}_h + \mathbf{b}_f \right)} -With the overall equilibrium vector being a mix of the free equilibria and the held values: -\deqn{\boldsymbol{\theta | \mathbf{\eta}_h} = \begin{bmatrix} \boldsymbol{\theta}_f \\\mathbf{\eta}_h\end{bmatrix}} +\deqn{\boldsymbol{\theta_f} = -\mathbf{A}_{ff}^{-1} \left( \mathbf{A}_{fh} + \mathbf{\eta}_h + \mathbf{b}_f \right)} + +With the overall equilibrium vector being a mix of the free equilibria and +the held values: + +\deqn{\boldsymbol{\theta | \mathbf{\eta}_h} = \begin{bmatrix} + \boldsymbol{\theta}_f \\\mathbf{\eta}_h\end{bmatrix}} } \examples{ \dontrun{ diff --git a/tests/testthat/test-coev_calculate_theta.R b/tests/testthat/test-coev_calculate_theta.R index a51a73c..3e8ea41 100644 --- a/tests/testthat/test-coev_calculate_theta.R +++ b/tests/testthat/test-coev_calculate_theta.R @@ -94,8 +94,15 @@ test_that("coev_calculate_theta() produces expected errors and output", { theta7 <- coev_calculate_theta(m7, list(w = NA, x = 0)) theta8 <- coev_calculate_theta(m8, list(x = NA, y = 0)) theta9 <- coev_calculate_theta(m9, list(x = NA, y = 0)) - theta1_null <- coev_calculate_theta(m1, NULL) - theta9_null <- coev_calculate_theta(m9, NULL) + theta1_null <- coev_calculate_theta(m1, intervention_values = NULL) + theta2_null <- coev_calculate_theta(m2, intervention_values = NULL) + theta3_null <- coev_calculate_theta(m3, intervention_values = NULL) + theta4_null <- coev_calculate_theta(m4, intervention_values = NULL) + theta5_null <- coev_calculate_theta(m5, intervention_values = NULL) + theta6_null <- coev_calculate_theta(m6, intervention_values = NULL) + theta7_null <- coev_calculate_theta(m7, intervention_values = NULL) + theta8_null <- coev_calculate_theta(m8, intervention_values = NULL) + theta9_null <- coev_calculate_theta(m9, intervention_values = NULL) expect_no_error(theta1) expect_no_error(theta2) expect_no_error(theta3) @@ -106,6 +113,13 @@ test_that("coev_calculate_theta() produces expected errors and output", { expect_no_error(theta8) expect_no_error(theta9) expect_no_error(theta1_null) + expect_no_error(theta2_null) + expect_no_error(theta3_null) + expect_no_error(theta4_null) + expect_no_error(theta5_null) + expect_no_error(theta6_null) + expect_no_error(theta7_null) + expect_no_error(theta8_null) expect_no_error(theta9_null) # output should be matrix of posterior draws expect_true(methods::is(theta1, "matrix")) @@ -118,6 +132,13 @@ test_that("coev_calculate_theta() produces expected errors and output", { expect_true(methods::is(theta8, "matrix")) expect_true(methods::is(theta9, "matrix")) expect_true(methods::is(theta1_null, "matrix")) + expect_true(methods::is(theta2_null, "matrix")) + expect_true(methods::is(theta3_null, "matrix")) + expect_true(methods::is(theta4_null, "matrix")) + expect_true(methods::is(theta5_null, "matrix")) + expect_true(methods::is(theta6_null, "matrix")) + expect_true(methods::is(theta7_null, "matrix")) + expect_true(methods::is(theta8_null, "matrix")) expect_true(methods::is(theta9_null, "matrix")) # output column names should be equal to variable names expect_true(identical(colnames(theta1), names(m1$variables))) @@ -129,8 +150,15 @@ test_that("coev_calculate_theta() produces expected errors and output", { expect_true(identical(colnames(theta7), names(m7$variables))) expect_true(identical(colnames(theta8), names(m8$variables))) expect_true(identical(colnames(theta9), names(m9$variables))) - expect_true(identical(colnames(theta9_null), names(m9$variables))) expect_true(identical(colnames(theta1_null), names(m1$variables))) + expect_true(identical(colnames(theta2_null), names(m2$variables))) + expect_true(identical(colnames(theta3_null), names(m3$variables))) + expect_true(identical(colnames(theta4_null), names(m4$variables))) + expect_true(identical(colnames(theta5_null), names(m5$variables))) + expect_true(identical(colnames(theta6_null), names(m6$variables))) + expect_true(identical(colnames(theta7_null), names(m7$variables))) + expect_true(identical(colnames(theta8_null), names(m8$variables))) + expect_true(identical(colnames(theta9_null), names(m9$variables))) # variables should be correctly held constant in output expect_true(all(theta1[,c("v","w","x","y")] == 0)) expect_true(all(theta2[,"x"] == 0))