Skip to content

Commit

Permalink
paired data analysis ignored for the moment, next cycle of development
Browse files Browse the repository at this point in the history
  • Loading branch information
@snaketron
snaketron committed Mar 28, 2024
1 parent 309f546 commit e4b92fb
Show file tree
Hide file tree
Showing 12 changed files with 23 additions and 757 deletions.
9 changes: 3 additions & 6 deletions R/dgu.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,17 +6,15 @@ DGU <- function(ud,
mcmc_cores = 1,
hdi_lvl = 0.95,
adapt_delta = 0.95,
max_treedepth = 12,
paired = FALSE) {
max_treedepth = 12) {

# check inputs
check_dgu_input(ud = ud,
mcmc_chains = as.integer(x = mcmc_chains),
mcmc_cores = as.integer(x = mcmc_cores),
mcmc_steps = as.integer(x = mcmc_steps),
mcmc_warmup = as.integer(x = mcmc_warmup),
hdi_lvl = hdi_lvl,
paired = paired)
hdi_lvl = hdi_lvl)

udr <- ud
ud <- get_usage(u = udr)
Expand All @@ -28,8 +26,7 @@ DGU <- function(ud,
# get model
m <- get_model(has_conditions = ud$has_conditions,
has_replicates = ud$has_replicates,
has_balanced_replicates = ud$has_balanced_replicates,
paired = paired)
has_balanced_replicates = ud$has_balanced_replicates)

# fit model
glm <- sampling(object = m$model,
Expand Down
9 changes: 3 additions & 6 deletions R/loo.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,17 +8,15 @@ LOO <- function(ud,
mcmc_cores = 1,
hdi_lvl = 0.95,
adapt_delta = 0.95,
max_treedepth = 12,
paired = FALSE) {
max_treedepth = 12) {

# check inputs
check_dgu_input(ud = ud,
mcmc_chains = as.integer(x = mcmc_chains),
mcmc_cores = as.integer(x = mcmc_cores),
mcmc_steps = as.integer(x = mcmc_steps),
mcmc_warmup = as.integer(x = mcmc_warmup),
hdi_lvl = hdi_lvl,
paired = paired)
hdi_lvl = hdi_lvl)

# process data
udp <- get_usage(u = ud)
Expand Down Expand Up @@ -59,8 +57,7 @@ LOO <- function(ud,
mcmc_cores = mcmc_cores,
hdi_lvl = hdi_lvl,
adapt_delta = adapt_delta,
max_treedepth = max_treedepth,
paired = paired)
max_treedepth = max_treedepth)

if(is.data.frame(out$gu)==TRUE) {
out$gu$loo_id <- rs[r]
Expand Down
21 changes: 2 additions & 19 deletions R/utils_input.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,16 +8,14 @@ check_dgu_input <- function(ud,
mcmc_cores,
mcmc_steps,
mcmc_warmup,
hdi_lvl,
paired) {
hdi_lvl) {

if(missing(ud) || is.null(ud) ||
missing(mcmc_chains) || is.null(mcmc_chains) ||
missing(mcmc_steps) || is.null(mcmc_steps) ||
missing(mcmc_warmup) || is.null(mcmc_warmup) ||
missing(mcmc_cores) || is.null(mcmc_cores) ||
missing(hdi_lvl) || is.null(hdi_lvl) ||
missing(paired) || is.null(paired)) {
missing(hdi_lvl) || is.null(hdi_lvl)) {
stop("arguments must be specified")
}

Expand All @@ -27,7 +25,6 @@ check_dgu_input <- function(ud,
check_mcmc_chains(mcmc_chains = mcmc_chains)
check_mcmc_cores(mcmc_cores = mcmc_cores)
check_hdi(hdi_lvl = hdi_lvl)
check_paired(paired = paired)
}


Expand Down Expand Up @@ -171,17 +168,3 @@ check_ud <- function(ud) {
}
}
}



# Description:
# paired input check
check_paired <- function(paired) {
if(length(paired) != 1) {
stop('paired must be a logical (TRUE or FALSE)')
}

if(is.logical(paired) == FALSE) {
stop('paired must be a logical (TRUE or FALSE)')
}
}
81 changes: 7 additions & 74 deletions R/utils_usage.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,7 @@ get_usage <- function(u) {
# if replicate column is provided BUT only one replicate is available
# per individual -> do analysis without replicates
k <- u[duplicated(u[,c("individual_id","condition","replicate_id")])==FALSE,
c("individual_id")]
c("individual_id")]
if(all(table(k)==1)) {
has_replicates <- FALSE
}
Expand Down Expand Up @@ -131,34 +131,14 @@ get_usage <- function(u) {
}



# Description:
# checks for structural problems in the data, and mismatch with given inputs
check_ud_content <- function(ud, paired) {
has_conditions <- ud$has_conditions
has_replicates <- ud$has_replicates
ud <- M$ud$proc_ud

table(ud$sample_id)
diag(table(ud$sample_id, ud$individual_id))
identical(diag(table(ud$sample_id, ud$individual_id)))
}



# Description:
# get the appropriate model
get_model <- function(has_replicates,
has_conditions,
has_balanced_replicates,
paired) {
has_balanced_replicates) {

model_type <- ifelse(test = has_conditions, yes = "DGU", no = "GU")

if(paired == TRUE & has_balanced_replicates == FALSE) {
stop("For paired analysis with replicates, you need balanced replicates!")
}

if(model_type == "GU") {
if(has_replicates) {
model <- stanmodels$gu_rep
Expand Down Expand Up @@ -187,18 +167,6 @@ get_model <- function(has_replicates,
"Yhat_rep", "Yhat_rep_prop", "Yhat_condition_prop",
"log_lik", "dgu", "dgu_prob", "theta")
model_name <- "DGU_rep"
if(paired) {
model <- stanmodels$dgu_paired_rep
pars <- c("phi", "kappa", "alpha",
"sigma_condition", "sigma_individual", "sigma_alpha",
"sigma_alpha_rep", "sigma_beta_rep",
"alpha_sample", "beta_sample",
"alpha_individual", "beta_individual",
"beta_condition",
"Yhat_rep", "Yhat_rep_prop", "Yhat_condition_prop",
"log_lik", "dgu", "dgu_prob", "theta")
model_name <- "DGU_paired_rep"
}
}
else {
model <- stanmodels$dgu
Expand All @@ -208,15 +176,6 @@ get_model <- function(has_replicates,
"Yhat_rep", "Yhat_rep_prop", "Yhat_condition_prop",
"log_lik", "dgu", "dgu_prob", "theta")
model_name <- "DGU"
if(paired) {
model <- stanmodels$dgu_paired
pars <- c("phi", "kappa", "alpha",
"sigma_condition", "sigma_individual", "sigma_alpha",
"alpha_individual", "beta_individual", "beta_condition",
"Yhat_rep", "Yhat_rep_prop", "Yhat_condition_prop",
"log_lik", "dgu", "dgu_prob", "theta")
model_name <- "DGU_paired"
}
}
}

Expand All @@ -225,22 +184,18 @@ get_model <- function(has_replicates,
model_type = model_type,
pars = pars,
has_replicates = has_replicates,
has_conditions = has_conditions,
paired = paired))
has_conditions = has_conditions))
}



# Description:
# get the appropriate model
get_model_debug <- function(has_replicates, has_conditions,
has_balanced_replicates, paired) {
get_model_debug <- function(has_replicates,
has_conditions,
has_balanced_replicates) {
model_type <- ifelse(test = has_conditions, yes = "DGU", no = "GU")

if(paired == TRUE & has_balanced_replicates == FALSE) {
stop("For paired analysis with replicates, you need balanced replicates!")
}

if(model_type == "GU") {
if(has_replicates) {
model <- rstan::stan_model(file = "inst/stan/gu_rep.stan")
Expand Down Expand Up @@ -270,18 +225,6 @@ get_model_debug <- function(has_replicates, has_conditions,
"Yhat_rep", "Yhat_rep_prop", "Yhat_condition_prop",
"log_lik", "dgu", "dgu_prob", "theta")
model_name <- "DGU_rep"
if(paired) {
model <- rstan::stan_model(file = "inst/stan/dgu_paired_rep.stan")
pars <- c("phi", "kappa", "alpha",
"sigma_condition", "sigma_individual", "sigma_alpha",
"sigma_alpha_rep", "sigma_beta_rep",
"alpha_sample", "beta_sample",
"alpha_individual", "beta_individual",
"beta_condition",
"Yhat_rep", "Yhat_rep_prop", "Yhat_condition_prop",
"log_lik", "dgu", "dgu_prob", "theta")
model_name <- "DGU_paired_rep"
}
}
else {
model <- rstan::stan_model(file = "inst/stan/dgu.stan")
Expand All @@ -291,15 +234,6 @@ get_model_debug <- function(has_replicates, has_conditions,
"Yhat_rep", "Yhat_rep_prop", "Yhat_condition_prop",
"log_lik", "dgu", "dgu_prob", "theta")
model_name <- "DGU"
if(paired) {
model <- rstan::stan_model(file = "inst/stan/dgu_paired.stan")
pars <- c("phi", "kappa", "alpha",
"sigma_condition", "sigma_individual", "sigma_alpha",
"alpha_individual", "beta_individual", "beta_condition",
"Yhat_rep", "Yhat_rep_prop", "Yhat_condition_prop",
"log_lik", "dgu", "dgu_prob", "theta")
model_name <- "DGU_paired"
}
}
}

Expand All @@ -308,6 +242,5 @@ get_model_debug <- function(has_replicates, has_conditions,
model_type = model_type,
pars = pars,
has_replicates = has_replicates,
has_conditions = has_conditions,
paired = paired))
has_conditions = has_conditions))
}
Binary file modified data/d_zibb_5.RData
View file
Binary file not shown.
Loading

0 comments on commit e4b92fb

Please sign in to comment.