From 77e62bd30b9a54c07fbc9a6ebd5ca5fa6db90a24 Mon Sep 17 00:00:00 2001 From: Josh Sumner <51797700+joshqsumner@users.noreply.github.com> Date: Tue, 5 Aug 2025 11:27:32 -0500 Subject: [PATCH 1/2] consistent interactions in frem --- NEWS.md | 4 +++ R/frem.R | 76 ++++++++++++++++++-------------------------------------- 2 files changed, 28 insertions(+), 52 deletions(-) diff --git a/NEWS.md b/NEWS.md index 062ac1c7..31304421 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,7 @@ +# pcvr 1.3.0.4 + +Simplified logic in `frem` to always use an interaction term between design variables. + # pcvr 1.3.0.3 Adding option to return exact parameters per each subject simulated in `growthSim`. diff --git a/R/frem.R b/R/frem.R index dc02ef79..4fc3790c 100644 --- a/R/frem.R +++ b/R/frem.R @@ -90,13 +90,7 @@ frem <- function(df, des, phenotypes, timeCol = NULL, cor = TRUE, returnData = F df <- formatted$data timeCol <- formatted$timeCol #* `Make formulas` - ext <- FALSE - if (length(des) == 2) { - ind_fmla <- paste0("(1|", des[1], ")+(1|", des[2], ")+(1|", des[1], ":", des[2], ")") - } else { - ind_fmla <- paste(paste0("(1|", des, ")"), collapse = "+") - ext <- TRUE - } + ind_fmla <- paste(paste0("(1|", c(des, paste(des, collapse = ":")), ")"), collapse = "+") #* `Find time and subset data` if (is.null(time)) { dat <- na.omit(df[df[[timeCol]] == max(df[[timeCol]]), c(des, phenotypes, timeCol)]) @@ -111,33 +105,18 @@ frem <- function(df, des, phenotypes, timeCol = NULL, cor = TRUE, returnData = F LONGITUDINAL <- TRUE dat <- na.omit(df[, c(des, phenotypes, timeCol)]) } - #* `Partition Variance` + H2 <- .partitionVarianceFrem(dat, timeCol, phenotypes, ind_fmla, des, ...) + colnames(H2) <- c(des, "Interaction", "Unexplained", timeCol, "singular", "Phenotypes") - H2 <- .partitionVarianceFrem(dat, timeCol, phenotypes, ind_fmla, ext, des, ...) - - if (!ext) { - colnames(H2) <- c(des[1], des[2], "Interaction", "Unexplained", timeCol, "singular", "Phenotypes") - } else { - colnames(H2) <- c(des, "Unexplained", timeCol, "singular", "Phenotypes") - } ordering <- H2[H2[[timeCol]] == max(H2[[timeCol]]), ] H2$Phenotypes <- ordered(H2$Phenotypes, levels = ordering$Phenotypes[order(ordering$Unexplained)]) h2_melt <- data.frame(data.table::melt(as.data.table(H2), id = c("Phenotypes", "singular", timeCol))) - - if (!ext) { - h2_melt$variable <- ordered(h2_melt$variable, - levels = c("Unexplained", des[1], des[2], "Interaction") - ) - } else { - h2_melt$variable <- ordered(h2_melt$variable, - levels = c("Unexplained", des) - ) - } + h2_melt$variable <- ordered(h2_melt$variable, + levels = c("Unexplained", des, "Interaction") + ) anova_dat <- h2_melt - #* `Plot Variance` - plotHelperOutputs <- .fremPlotHelper( LONGITUDINAL, anova_dat, markSingular, dummyX, timeCol, dat, phenotypes, cor, combine @@ -277,7 +256,7 @@ frem <- function(df, des, phenotypes, timeCol = NULL, cor = TRUE, returnData = F #' @keywords internal #' @noRd -.partitionVarianceFrem <- function(dat, timeCol, phenotypes, ind_fmla, ext, des, ...) { +.partitionVarianceFrem <- function(dat, timeCol, phenotypes, ind_fmla, des, ...) { H2 <- data.frame(do.call(rbind, lapply(sort(unique(dat[[timeCol]])), function(tm) { sub <- dat[dat[[timeCol]] == tm, ] des_bools <- unlist(lapply(des, function(var) { @@ -299,32 +278,25 @@ frem <- function(df, des, phenotypes, timeCol = NULL, cor = TRUE, returnData = F re <- lme4::VarCorr(model) res <- attr(lme4::VarCorr(model), "sc")^2 - if (!ext) { - interaction.var <- as.numeric(attr(re[[which(grepl(":", names(re)))]], "stddev"))^2 - des1.var <- as.numeric(attr(re[[des[1]]], "stddev"))^2 - des2.var <- as.numeric(attr(re[[des[2]]], "stddev"))^2 - - tot.var <- sum(as.numeric(re), res) - unexp <- 1 - sum(as.numeric(re)) / sum(as.numeric(re), res) - - h2 <- c( - (des1.var / tot.var), - (des2.var / tot.var), - (interaction.var / tot.var), - unexp, - tm, - singular + interaction.var <- as.numeric(attr(re[[which(grepl(":", names(re)))]], "stddev"))^2 + des.var <- unlist( + lapply( + des, + function(ds) { + return(as.numeric(attr(re[[ds]], "stddev"))^2) + } ) - } else { - var <- lapply(des, function(i) { - return(as.numeric(attr(re[[i]], "stddev"))^2) - }) - - tot.var <- sum(as.numeric(re), res) - unexp <- 1 - sum(as.numeric(re)) / sum(as.numeric(re), res) + ) + tot.var <- sum(as.numeric(re), res) + unexp <- 1 - sum(as.numeric(re)) / sum(as.numeric(re), res) - h2 <- c(unlist(var) / tot.var, unexp, tm, singular) - } + h2 <- c( + (des.var / tot.var), + (interaction.var / tot.var), + unexp, + tm, + singular + ) return(h2) })) return(pheno_df) From 4d195d6106f2d42fca68165fb1b9bcfd834282fb Mon Sep 17 00:00:00 2001 From: Josh Sumner <51797700+joshqsumner@users.noreply.github.com> Date: Tue, 5 Aug 2025 13:15:45 -0500 Subject: [PATCH 2/2] adapting to work with 1L design variables --- R/frem.R | 21 +++++++++++++++------ 1 file changed, 15 insertions(+), 6 deletions(-) diff --git a/R/frem.R b/R/frem.R index 4fc3790c..d4d8b830 100644 --- a/R/frem.R +++ b/R/frem.R @@ -90,7 +90,11 @@ frem <- function(df, des, phenotypes, timeCol = NULL, cor = TRUE, returnData = F df <- formatted$data timeCol <- formatted$timeCol #* `Make formulas` - ind_fmla <- paste(paste0("(1|", c(des, paste(des, collapse = ":")), ")"), collapse = "+") + int <- paste(des, collapse = ":") + if (!grepl(":", int)) { + int <- NULL + } # no interaction for 1L des + ind_fmla <- paste(paste0("(1|", c(des, int), ")"), collapse = "+") #* `Find time and subset data` if (is.null(time)) { dat <- na.omit(df[df[[timeCol]] == max(df[[timeCol]]), c(des, phenotypes, timeCol)]) @@ -107,7 +111,6 @@ frem <- function(df, des, phenotypes, timeCol = NULL, cor = TRUE, returnData = F } #* `Partition Variance` H2 <- .partitionVarianceFrem(dat, timeCol, phenotypes, ind_fmla, des, ...) - colnames(H2) <- c(des, "Interaction", "Unexplained", timeCol, "singular", "Phenotypes") ordering <- H2[H2[[timeCol]] == max(H2[[timeCol]]), ] H2$Phenotypes <- ordered(H2$Phenotypes, levels = ordering$Phenotypes[order(ordering$Unexplained)]) @@ -277,8 +280,6 @@ frem <- function(df, des, phenotypes, timeCol = NULL, cor = TRUE, returnData = F } re <- lme4::VarCorr(model) res <- attr(lme4::VarCorr(model), "sc")^2 - - interaction.var <- as.numeric(attr(re[[which(grepl(":", names(re)))]], "stddev"))^2 des.var <- unlist( lapply( des, @@ -287,12 +288,15 @@ frem <- function(df, des, phenotypes, timeCol = NULL, cor = TRUE, returnData = F } ) ) + if (length(des) > 1) { + interaction.var <- as.numeric(attr(re[[which(grepl(":", names(re)))]], "stddev"))^2 + des.var <- c(des.var, interaction.var) + } tot.var <- sum(as.numeric(re), res) - unexp <- 1 - sum(as.numeric(re)) / sum(as.numeric(re), res) + unexp <- 1 - (sum(as.numeric(re)) / sum(as.numeric(re), res)) h2 <- c( (des.var / tot.var), - (interaction.var / tot.var), unexp, tm, singular @@ -302,5 +306,10 @@ frem <- function(df, des, phenotypes, timeCol = NULL, cor = TRUE, returnData = F return(pheno_df) }))) H2$Phenotypes <- rep(phenotypes, length.out = nrow(H2)) + des_colnames <- des + if (length(des) > 1) { + des_colnames <- c(des, "Interaction") + } + colnames(H2) <- c(des_colnames, "Unexplained", timeCol, "singular", "Phenotypes") return(H2) }