Merge pull request #162 from danforthcenter/frem-interactions

consistent interactions in frem

Author: joshqsumner

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`.

R/frem.R

@@ -90,13 +90,11 @@ 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
-  }
+  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)])
@@ -111,33 +109,17 @@ 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, ...)
 
-  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 +259,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) {
@@ -298,37 +280,36 @@ 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
+      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)
-
-        h2 <- c(unlist(var) / tot.var, unexp, tm, singular)
+      )
+      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))
+
+      h2 <- c(
+        (des.var / tot.var),
+        unexp,
+        tm,
+        singular
+      )
       return(h2)
     }))
     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)
 }