R/linear_outcome_envir_interaction_mle_function.R

Defines functions calc.like.linear.log.envir.interaction linear.outcome.log.envir.interaction.sds linear.mles.log.envir.interaction p_vec_returner X_mat_returner

Documented in calc.like.linear.log.envir.interaction linear.mles.log.envir.interaction linear.outcome.log.envir.interaction.sds p_vec_returner X_mat_returner

# MLE Functions for Normal/Continuous Outcomes with logistic environment interaction
#' Function to output X matrices used in calculation of MLE's for linear outcome with logistic environment interaction
#'
#' Returns X matrices used in calculation of MLE's for linear outcome with logistic environment interaction
#'
#' @import MASS
#'
#' @param mod type of model
#' @param reduced logical, indicates whether the X matrix will be used for a reduced model
#'
#' @return A matrix to be used in MLE calculation for linear outcome with logistic environment interaction
#'
#' @examples
#' X_mat_returner(mod = "Dominant")
#'
#' @export
#'
X_mat_returner <- function(mod, reduced = F)
{

	X_mat_dom <- rbind(
		c(1, 0, 0, 0),
		c(1, 1, 0, 0),
		c(1, 1, 0, 0),
		c(1, 0, 1, 0),
		c(1, 1, 1, 1),
		c(1, 1, 1, 1))

	X_mat_rec <- rbind(
		c(1, 0, 0, 0),
		c(1, 0, 0, 0),
		c(1, 1, 0, 0),
		c(1, 0, 1, 0),
		c(1, 0, 1, 0),
		c(1, 1, 1, 1))

	X_mat_add <- rbind(
		c(1, 0, 0, 0),
		c(1, 1, 0, 0),
		c(1, 2, 0, 0),
		c(1, 0, 1, 0),
		c(1, 1, 1, 1),
		c(1, 2, 1, 2))

	X_mat_2df <- rbind(
		c(1, 0, 0, 0, 0, 0),
		c(1, 1, 0, 0, 0, 0),
		c(1, 0, 1, 0, 0, 0),
		c(1, 0, 0, 1, 0, 0),
		c(1, 1, 0, 1, 1, 0),
		c(1, 0, 1, 1, 0, 1))
	if(reduced) X_mat_2df <- X_mat_2df[,1:(ncol(X_mat_2df) - 1)]

	X_mat_null <- rbind(
		c(1, 0, 0, 0),
		c(1, 0, 0, 0),
		c(1, 0, 0, 0),
		c(1, 0, 0, 0),
		c(1, 0, 0, 0),
		c(1, 0, 0, 0))

	X_mat_list <- list("Dominant" = X_mat_dom, "Recessive" = X_mat_rec, "Additive" = X_mat_add, "2df" = X_mat_2df, "null" = X_mat_null)

	if(reduced) X_mat_list <- lapply(X_mat_list, function(amat) amat <- amat[,1:(ncol(amat) - 1)])

	return(X_mat_list[[mod]])
}

#' Function to output probability vector used in calculation of MLE's for linear outcome with logistic environment interaction
#'
#' Returns probability vector used in calculation of MLE's for linear outcome with logistic environment interaction
#'
#' @param MAF Minor allele frequency
#' @param P_e Population prevalence of logistic environmental factor
#'
#' @return A probability vector to be used in MLE calculation for linear outcome with logistic environment interaction
#'
#' @examples 
#' p_vec_returner(MAF = 0.1, P_e = 0.2)
#'
#' @export
#'
p_vec_returner <- function(MAF, P_e)
{
	p_AA0 <- (1 - MAF)^2 * (1 - P_e)
	p_AB0 <- 2*(1 - MAF)*MAF * (1 - P_e)
	p_BB0 <- MAF^2 *(1 - P_e)
	p_AA1 <- (1 - MAF)^2 * P_e
	p_AB1 <- 2*(1 - MAF)*MAF * P_e
	p_BB1 <- MAF^2 * P_e
	p_vec <- c(p_AA0, p_AB0, p_BB0, p_AA1, p_AB1, p_BB1)
	return(p_vec)
}

#' Function to calculate the standard deviation of y given x for linear models with logistic environment interaction
#'
#' Returns the standard deviation of y given x for linear models with logistic environment interaction
#'
#' @param MAF Minor allele Frequency
#' @param P_e Population prevalence of logistic environmental factor
#' @param ES_G Genetic Effect size
#' @param ES_E Environment Effect size
#' @param ES_GE Environment x Genetic interaction Effect size
#' @param Test.Model Test model
#' @param True.Model True model
#' @param reduced logical, indicates whether the X matrix will be used for a reduced model
#'
#' @return The standard deviation of y given x for linear models with logistic environment interaction
#'
#' @examples 
#' linear.mles.log.envir.interaction(MAF = 0.1, P_e = 0.2, 
#' 	ES_G = 1.2, ES_E = 1.3, ES_GE = 2, 
#' 	Test.Model = "Dominant", True.Model = "Additive")
#'
#' @export
#'
linear.mles.log.envir.interaction <- function(MAF, P_e, ES_G, ES_E, ES_GE, Test.Model, True.Model, reduced = F)
{

	# create elements we will use to calculate MLE
	p_vec <- p_vec_returner(MAF, P_e)
	W_mat <- diag(p_vec)
	X_matF <- X_mat_returner(Test.Model, reduced = reduced) # the MLE is calculated using the test model

	# the effect sizes are calculated from the true model
	ES_vec <- (X_mat_returner(True.Model) %*% c(0, ES_G, ES_E, ES_GE))[,1]

	beta_hat <- as.vector(ginv(crossprod(X_matF, W_mat) %*% X_matF) %*% crossprod(X_matF, W_mat) %*% ES_vec)
	return(beta_hat)
}

#' Function to calculate the standard deviation of y given x for linear models with logistic environment interaction
#'
#' Returns the standard deviation of y given x for linear models with logistic environment interaction
#'
#' @param MAF Minor allele Frequency
#' @param P_e Population prevalence of logistic environmental factor
#' @param ES_G Genetic Effect size
#' @param ES_E Environment Effect size
#' @param ES_GE Environment x Genetic interaction Effect size
#' @param mod Test model
#' @param True.Model True model
#' @param sd_y Standard deviation of y
#' @param reduced logical, indicates whether the X matrix will be used for a reduced model
#'
#' @return The standard deviation of y given x for linear models with logistic environment interaction
#'
#' @examples
#' linear.outcome.log.envir.interaction.sds(MAF = 0.1, P_e = 0.2, sd_y = 10,
#' 	ES_G = 1.2, ES_E = 1.3, ES_GE = 2, mod = "Dominant", True.Model = "Additive")
#'
#' @export
#'
linear.outcome.log.envir.interaction.sds <- function(MAF, P_e, ES_G, ES_E, ES_GE, mod, True.Model, sd_y, reduced = F)
{

	p_vec <- p_vec_returner(MAF, P_e)

	ES_test <- linear.mles.log.envir.interaction(MAF = MAF, P_e = P_e, ES_G = ES_G, 
					ES_E = ES_E, ES_GE = ES_GE, Test.Model = mod, True.Model = True.Model, reduced = reduced)


	ES_vec <- (X_mat_returner(mod, reduced = reduced) %*% ES_test)[,1]
	ES_vec_truemodel  <- (X_mat_returner(True.Model) %*% c(0, ES_G, ES_E, ES_GE))[,1]

	y_bar <- p_vec_returner(MAF, P_e) * ES_vec_truemodel
	y_bar_sum <- sum(y_bar)

	sd_y_x <- sqrt(sd_y^2  - sum(p_vec * (ES_vec - y_bar_sum)^2))

	return(sd_y_x)
}

#' Function to calculate the standard deviation of y given x for linear models with logistic environment interaction
#'
#' Returns the standard deviation of y given x for linear models with logistic environment interaction
#'
#' @param beta_hat Effect sizes from MLE
#' @param MAF Minor allele Frequency
#' @param P_e Population prevalence of logistic environmental factor
#' @param ES_G Genetic Effect size
#' @param ES_E Environment Effect size
#' @param ES_GE Environment x Genetic interaction Effect size
#' @param sd_y_x_truth Standard deviation of y for the true model
#' @param sd_y_x_model Standard deviation of y for the test model
#' @param Test.Model Test model
#' @param True.Model True model
#' @param reduced logical, indicates whether the X matrix will be used for a reduced model
#'
#' @return The standard deviation of y given x for linear models with logistic environment interaction
#'
#' @examples
#' beta_hat = linear.mles.log.envir.interaction(MAF = 0.1, P_e = 0.2, 
#' 	ES_G = 1.2, ES_E = 1.3, ES_GE = 2, 
#' 	Test.Model = "Dominant", True.Model = "Additive")
#' calc.like.linear.log.envir.interaction(beta_hat = beta_hat,
#' 	MAF = 0.1, P_e = 0.2, ES_G = 1.2, ES_E = 1.3,
#' 	ES_GE = 2, sd_y_x_truth = 9.947945, sd_y_x_model = 9.949468, 
#' 	True.Model = "Additive", Test.Model="Dominant")
#'
#' @export
#'
calc.like.linear.log.envir.interaction <- function(beta_hat, MAF, P_e, ES_G, ES_E, ES_GE, sd_y_x_truth, sd_y_x_model, Test.Model, True.Model, reduced  = F)
{
	p_vec <- p_vec_returner(MAF, P_e)

	betavec <- (X_mat_returner(Test.Model, reduced = reduced) %*% beta_hat)[,1]
	ES_vec <- (X_mat_returner(True.Model) %*% c(0, ES_G, ES_E, ES_GE))[,1]


	ll <- sum(p_vec * apply(cbind(ES_vec, betavec), 1, function(x) 
		expected.linear.ll(mean_truth = x[1], mean_model = x[2], sd_y_x_truth, sd_y_x_model)))

	return(ll)
}

Try the genpwr package in your browser

Any scripts or data that you put into this service are public.

genpwr documentation built on March 31, 2021, 1:06 a.m.