DImodels v1.3

rafamoral · Nov 21, 2023 · 047c027 · 047c027
1 parent 872e549
commit 047c027
Show file tree

Hide file tree

Showing 17 changed files with 1,355 additions and 410 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: DImodels
 Type: Package
 Title: Diversity-Interactions (DI) Models
-Version: 1.2.1
+Version: 1.3
 Date: 2023-05-20
 Authors@R: c(person("Rafael", "de Andrade Moral", role = c("aut", "cre"), email = "rafael.deandrademoral@mu.ie"), person("John", "Connolly", role = "aut"), person("Rishabh", "Vishwakarma", role = "ctb"), person("Caroline", "Brophy", role = "aut"))
 Author: Rafael de Andrade Moral [aut, cre], John Connolly [aut], Rishabh Vishwakarma [ctb], Caroline Brophy [aut]
@@ -20,3 +20,4 @@ Description: The 'DImodels' package is suitable for analysing data from biodiver
 License: GPL (>=2)
 NeedsCompilation: no
 Config/testthat/edition: 3
+RoxygenNote: 7.2.3
diff --git a/R/DI.R b/R/DI.R
diff --git a/R/DI_internal.R b/R/DI_internal.R
@@ -44,14 +44,14 @@ DI_check_and_fit <- function(fmla, y, block, density, treat, family, data, FG) {
   return(mod)
 }
 
-proflik_theta <- function(theta, obj, family, int_terms, DImodel, nSpecies, FGnames) {
+proflik_theta <- function(theta, obj, family, int_terms, prop, DImodel, nSpecies, FGnames) {
   mm <- model.matrix(obj)
 
   if(DImodel %in% c("E","AV")) {
     data_theta_E_AV <- obj$data
     data_theta <- data.frame(mm)
-    data_theta[1:nSpecies] <- data_theta[1:nSpecies]^theta
-    new_E_AV <- DI_data_E_AV_internal(prop = 1:nSpecies, data = data_theta)
+    #data_theta[1:nSpecies] <- data_theta[1:nSpecies]^theta
+    new_E_AV <- DI_data_E_AV(prop = prop, data = data_theta_E_AV, theta = theta)
     data_theta_E_AV$E <- new_E_AV$E
     data_theta_E_AV$AV <- new_E_AV$AV
     fitted_model_theta <- glm(formula(obj), family = family, data = data_theta_E_AV)
@@ -64,10 +64,14 @@ proflik_theta <- function(theta, obj, family, int_terms, DImodel, nSpecies, FGna
     sigma_hat <- sqrt(sum(fitted_model_theta$residuals^2)/(fitted_model_theta$df.residual) * (n - p)/n)
     llik <- sum(dnorm(fitted_model_theta$y, mu_hat, sigma_hat, log = TRUE))
   } else if(DImodel == "FG") {
-    data_theta_FG <- obj$data
-    data_theta <- data.frame(mm)
-    data_theta[1:nSpecies] <- data_theta[1:nSpecies]^theta
-    new_FG <- DI_data_FG_internal(prop = 1:nSpecies, FG = FGnames, data = data_theta)
+    #data_theta_FG <- obj$data
+    data_theta_FG <- data.frame(mm, check.names = FALSE)
+    colnames(data_theta_FG) <- gsub("`", "", colnames(data_theta_FG))
+    # Add original props to calculate interactions
+    data_theta_FG <- cbind(data_theta_FG, obj$data[, prop])
+
+    #data_theta[1:nSpecies] <- data_theta[1:nSpecies]^theta
+    new_FG <- DI_data_FG(prop = prop, FG = FGnames, data = data_theta_FG[, prop], theta = theta)
     FG_ <- new_FG$FG
     ## if we have column names starting with FG_ already
     ## in data_theta_FG, then substitute columns else it's all good
@@ -84,6 +88,10 @@ proflik_theta <- function(theta, obj, family, int_terms, DImodel, nSpecies, FGna
     }
     old_formula <- formula(obj)
     new_formula <- paste0(old_formula[2], " ~ ", old_formula[3])
+
+    data_theta_FG$y <- obj$y
+    names(data_theta_FG)[length(names(data_theta_FG))] <- paste(old_formula[2])
+
     fitted_model_theta <- glm(as.formula(new_formula), family = family, data = data_theta_FG)
     #mu_hat <- fitted(fitted_model_theta)
     #sigma_hat <- sqrt(sum(fitted_model_theta$residuals^2)/(fitted_model_theta$df.residual-1))
@@ -95,8 +103,10 @@ proflik_theta <- function(theta, obj, family, int_terms, DImodel, nSpecies, FGna
     llik <- sum(dnorm(fitted_model_theta$y, mu_hat, sigma_hat, log = TRUE))
   } else if(DImodel == "ADD") {
     #data_theta_ADD <- obj$data
+    # Drop existing _add columns to avoid conflicts
     data_theta <- data.frame(mm, check.names = FALSE)
-    new_ADD <- DI_data_ADD_theta(prop = 1:nSpecies, data = data_theta, theta = theta)
+    data_theta <- cbind(data_theta, obj$data[, prop])
+    new_ADD <- DI_data_ADD_theta(prop = prop, data = data_theta, theta = theta)
     ADD_cols_in_the_data <- int_terms #grep("_add", colnames(data_theta))
     if(length(ADD_cols_in_the_data) > 0) {
       j <- 1
@@ -172,7 +182,8 @@ DI_theta <- function(obj, DImodel, FGnames, prop, nSpecies, family) {
   #options(warn = -1)
   upper_boundary <- 1.5
   theta_info <- get_theta_info(upper_boundary = upper_boundary, DImodel = DImodel,
-                               obj = obj, family = family, int_terms = int_terms,
+                               obj = obj, prop = prop, 
+                               family = family, int_terms = int_terms, 
                                nSpecies = nSpecies, FGnames = FGnames)
   #options(warn = 0)
   theta_hat <- theta_info$theta_hat
@@ -184,16 +195,18 @@ DI_theta <- function(obj, DImodel, FGnames, prop, nSpecies, family) {
   if(DImodel %in% c("E","AV")) {
     data_theta_E_AV <- obj$data
     data_theta <- data.frame(mm)
-    data_theta[1:nSpecies] <- data_theta[1:nSpecies]^theta_hat
-    new_E_AV <- DI_data_E_AV_internal(prop = 1:nSpecies, data = data_theta)
+    #data_theta[1:nSpecies] <- data_theta[1:nSpecies]^theta_hat
+    new_E_AV <- DI_data_E_AV(prop = prop, data = data_theta_E_AV, theta = theta_hat)
     data_theta_E_AV$E <- new_E_AV$E
     data_theta_E_AV$AV <- new_E_AV$AV
     mod_theta <- glm(formula(obj), family = family, data = data_theta_E_AV)
   } else if(DImodel == "FG") {
-    data_theta_FG <- obj$data
-    data_theta <- data.frame(mm)
-    data_theta[1:nSpecies] <- data_theta[1:nSpecies]^theta_hat
-    new_FG <- DI_data_FG_internal(prop = 1:nSpecies, FG = FGnames, data = data_theta)
+    #data_theta_FG <- obj$data
+    data_theta_FG <- data.frame(mm, check.names = FALSE)
+    colnames(data_theta_FG) <- gsub("`", "", colnames(data_theta_FG))
+    data_theta_FG <- cbind(data_theta_FG, obj$data[, prop])
+    #data_theta[1:nSpecies] <- data_theta[1:nSpecies]^theta_hat
+    new_FG <- DI_data_FG(prop = prop, FG = FGnames, theta = theta_hat, data = data_theta_FG)
     FG_ <- new_FG$FG
     ## if we have column names starting with FG_ already
     ## in data_theta_FG, then substitute columns, else it's all good
@@ -207,11 +220,14 @@ DI_theta <- function(obj, DImodel, FGnames, prop, nSpecies, family) {
     }
     old_formula <- formula(obj)
     new_formula <- paste0(old_formula[2], " ~ ", old_formula[3])
+    data_theta_FG$y <- obj$y
+    names(data_theta_FG)[length(names(data_theta_FG))] <- paste(old_formula[2])
     mod_theta <- glm(as.formula(new_formula), family = family, data = data_theta_FG)
   } else if(DImodel == "ADD") {
     #data_theta_ADD <- obj$data
     data_theta <- data.frame(mm, check.names = FALSE)
-    new_ADD <- DI_data_ADD_theta(prop = 1:nSpecies, data = data_theta, theta = theta_hat)
+    data_theta <- cbind(data_theta, obj$data[, prop])
+    new_ADD <- DI_data_ADD_theta(prop = prop, data = data_theta, theta = theta_hat)
     ADD_cols_in_the_data <- int_terms #grep("_add", colnames(data_theta))
     if(length(ADD_cols_in_the_data) > 0) {
       j <- 1
@@ -250,15 +266,15 @@ DI_theta <- function(obj, DImodel, FGnames, prop, nSpecies, family) {
   return(mod_theta)
 }
 
-get_theta_info <- function(upper_boundary, DImodel, obj, family, int_terms,
+get_theta_info <- function(upper_boundary, DImodel, obj, prop, family, int_terms,
                            nSpecies, FGnames) {
   optimum <- optimize(proflik_theta, interval = c(0.00001, upper_boundary), 
-                      maximum = TRUE, DImodel = DImodel, obj = obj, family = family, 
+                      maximum = TRUE, DImodel = DImodel, obj = obj, prop = prop, family = family, 
                       int_terms = int_terms, nSpecies = nSpecies, FGnames = FGnames)
   theta_hat <- optimum$maximum
   theta_grid <- seq(0.00001, upper_boundary + 1, length = 100)
   proflik_theta_vec <- Vectorize(proflik_theta, "theta")
-  profile_loglik <- proflik_theta_vec(theta = theta_grid, obj = obj, 
+  profile_loglik <- proflik_theta_vec(theta = theta_grid, obj = obj, prop = prop,
                                       family = family, int_terms = int_terms, DImodel = DImodel, 
                                       nSpecies = nSpecies, FGnames = FGnames)
 

diff --git a/R/autoDI.R b/R/autoDI.R
@@ -7,7 +7,7 @@
 # selection = "Ftest" to perform F (or X2) tests to select best model,
 #             "AIC" to use select model with smallest AIC
 
-autoDI <- function(y, prop, data, block, density, treat, FG = NULL, 
+autoDI <- function(y, prop, data, block, density, treat, ID, FG = NULL, 
                    selection = c("Ftest","AIC","AICc","BIC","BICc"), 
                    step0 = FALSE, step4 = TRUE) {
   if(missing(y)) stop("You must supply a response variable name or column index through the argument 'y'.\n")
@@ -57,20 +57,29 @@ autoDI <- function(y, prop, data, block, density, treat, FG = NULL,
   if(missing(FG)){
     FG <- NULL
   }
+
+  if(missing(ID)){
+    if(is.numeric(prop)){
+      ID <- paste0(colnames(data)[prop], "_ID")
+    }
+    else {
+      ID <- paste0(prop, "_ID")
+    }
+  }
 
   # Step 0
   if(step0) {
     autoDI_step0(y = y, block = block, density = density, treat = treat, family = family, data = data)
   }
 
   # Step 1
-  step1_model <- autoDI_step1(y = y, block = block, density = density, prop = prop, treat = treat, family = family, data = data, selection = selection)
+  step1_model <- autoDI_step1(y = y, block = block, density = density, prop = prop, treat = treat, ID = ID, family = family, data = data, selection = selection)
 
   selected_model <- step1_model$model
   theta_flag <- step1_model$theta_flag
 
   # step 2
-  step2_model <- autoDI_step2(y = y, block = block, density = density, prop = prop, treat = treat, family = family, FG = FG, data = data, theta_flag = theta_flag, selection = selection, selected_model = selected_model)
+  step2_model <- autoDI_step2(y = y, block = block, density = density, prop = prop, treat = treat, ID = ID, family = family, FG = FG, data = data, theta_flag = theta_flag, selection = selection, selected_model = selected_model)
 
   # step 3
   step3_model <- autoDI_step3(selected_model = step2_model, selection = selection, family = family)

diff --git a/R/autoDI_internal.R b/R/autoDI_internal.R
@@ -331,30 +331,30 @@ autoDI_step0 <- function(y, block, density, treat, family, data) {
   #options(scipen = 0)
 }
 
-autoDI_step1 <- function(y, block, density, prop, treat, family, data, selection) {
+autoDI_step1 <- function(y, block, density, prop, treat, ID, family, data, selection) {
 
   model_tag <- ifelse(length(prop) == 2, "FULL", "AV")
 
   if(is.na(treat)) {
       model_name <- paste0(model_tag, "_model")
       fit_theta <- suppressWarnings(suppressMessages(DI(y = y, block = block, density = density, prop = prop,
-                                       DImodel = model_tag, data = data, estimate_theta = TRUE)))
+                                                        ID = ID, DImodel = model_tag, data = data, estimate_theta = TRUE)))
       fit_notheta <- suppressWarnings(suppressMessages(DI(y = y, block = block, density = density, prop = prop,
-                                         DImodel = model_tag, data = data, estimate_theta = FALSE)))
+                                                          ID = ID, DImodel = model_tag, data = data, estimate_theta = FALSE)))
       fit_theta$DIcall <- call("DI", y = y, block = block, density = density, prop = prop,
-                               DImodel = model_tag, data = data, estimate_theta = TRUE)
+                               ID = ID, DImodel = model_tag, data = data, estimate_theta = TRUE)
       fit_notheta$DIcall <- call("DI", y = y, block = block, density = density, prop = prop,
-                                 DImodel = model_tag, data = data, estimate_theta = FALSE)
+                                 ID = ID, DImodel = model_tag, data = data, estimate_theta = FALSE)
     } else {
       model_name <- paste0(model_tag, "_model_treat")
       fit_theta <- suppressWarnings(suppressMessages(DI(y = y, block = block, density = density, prop = prop,
-                                       treat = treat, DImodel = model_tag, data = data, estimate_theta = TRUE)))
+                                                        ID = ID, treat = treat, DImodel = model_tag, data = data, estimate_theta = TRUE)))
       fit_notheta <- suppressWarnings(suppressMessages(DI(y = y, block = block, density = density, prop = prop,
-                                         treat = treat, DImodel = model_tag, data = data, estimate_theta = FALSE)))
+                                                          ID = ID, treat = treat, DImodel = model_tag, data = data, estimate_theta = FALSE)))
       fit_theta$DIcall <- call("DI", y = y, block = block, density = density, prop = prop,
-                               treat = treat, DImodel = model_tag, data = data, estimate_theta = TRUE)
+                               ID = ID, treat = treat, DImodel = model_tag, data = data, estimate_theta = TRUE)
       fit_notheta$DIcall <- call("DI", y = y, block = block, density = density, prop = prop,
-                                 treat = treat, DImodel = model_tag, data = data, estimate_theta = FALSE)
+                                 ID = ID, treat = treat, DImodel = model_tag, data = data, estimate_theta = FALSE)
     }
 
   message("\n", strrep("-", getOption("width")))
@@ -392,7 +392,7 @@ autoDI_step1 <- function(y, block, density, prop, treat, family, data, selection
 }
 
 autoDI_step2 <- function(y, block, density, 
-                         prop, treat, FG, 
+                         prop, treat, ID, FG, 
                          family, data, theta_flag,
                          selection, selected_model){
 
@@ -420,10 +420,10 @@ autoDI_step2 <- function(y, block, density,
   }
 
   ID_model <- suppressWarnings(suppressMessages(DI(y = y, block = block, density = density, 
-                 prop = prop, treat = treat, 
+                 prop = prop, treat = treat, ID = ID,
                  data = data, DImodel = 'ID')))
   STR_model <- suppressWarnings(suppressMessages(DI(y = y, block = block, density = density, 
-                  prop = prop, treat = treat, 
+                  prop = prop, treat = treat, ID = ID,
                   data = data, DImodel = 'STR')))
   model_list <- list('STR_model' = STR_model, 'ID_model' = ID_model)
 
@@ -432,28 +432,28 @@ autoDI_step2 <- function(y, block, density,
     theta_est <- coefficients(selected_model)['theta']
     if(AV_flag){
       AV_model <- suppressWarnings(suppressMessages(DI(y = y, block = block, density = density, 
-                     prop = prop, treat = treat, theta = theta_est, 
+                     prop = prop, treat = treat, ID = ID, theta = theta_est, 
                      data = data, DImodel = 'AV')))
       model_list[[length(model_list)+1]] <- AV_model
       names(model_list)[length(model_list)] <- 'AV_model'
     }
     if(FG_flag){
       FG_model <- suppressWarnings(suppressMessages(DI(y = y, block = block, density = density, 
-                     prop = prop, treat = treat, theta = theta_est,  
+                     prop = prop, treat = treat, ID = ID, theta = theta_est,  
                      data = data, FG = FG, DImodel = 'FG')))
       model_list[[length(model_list)+1]] <- FG_model
       names(model_list)[length(model_list)] <- 'FG_model'
     }
     if(ADD_flag){
       ADD_model <- suppressWarnings(suppressMessages(DI(y = y, block = block, density = density, 
-                      prop = prop, treat = treat, theta = theta_est, 
+                      prop = prop, treat = treat, ID = ID, theta = theta_est, 
                       data = data, DImodel = 'ADD')))
       model_list[[length(model_list)+1]] <- ADD_model
       names(model_list)[length(model_list)] <- 'ADD_model'
     }
     if(FULL_flag){
       FULL_model <- suppressWarnings(suppressMessages(DI(y = y, block = block, density = density, 
-                       prop = prop, treat = treat, theta = theta_est, 
+                       prop = prop, treat = treat, ID = ID, theta = theta_est, 
                        data = data, DImodel = 'FULL')))
       model_list[[length(model_list)+1]] <- FULL_model
       names(model_list)[length(model_list)] <- 'FULL_model'
@@ -467,28 +467,28 @@ autoDI_step2 <- function(y, block, density,
   } else {
     if(AV_flag){
       AV_model <- suppressWarnings(suppressMessages(DI(y = y, block = block, density = density, 
-                     prop = prop, treat = treat, estimate_theta = F,
+                     prop = prop, treat = treat, ID = ID, estimate_theta = F,
                      data = data, DImodel = 'AV')))
       model_list[[length(model_list)+1]] <- AV_model
       names(model_list)[length(model_list)] <- 'AV_model'
     }
     if(FG_flag){
       FG_model <- suppressWarnings(suppressMessages(DI(y = y, block = block, density = density, 
-                     prop = prop, treat = treat, estimate_theta = F, 
+                     prop = prop, treat = treat, ID = ID, estimate_theta = F, 
                      data = data, FG = FG, DImodel = 'FG')))
       model_list[[length(model_list)+1]] <- FG_model
       names(model_list)[length(model_list)] <- 'FG_model'
     }
     if(ADD_flag){
       ADD_model <- suppressWarnings(suppressMessages(DI(y = y, block = block, density = density, 
-                      prop = prop, treat = treat, estimate_theta = F, 
+                      prop = prop, treat = treat, ID = ID, estimate_theta = F, 
                       data = data, DImodel = 'ADD')))
       model_list[[length(model_list)+1]] <- ADD_model
       names(model_list)[length(model_list)] <- 'ADD_model'
     } 
     if(FULL_flag){
       FULL_model <- suppressWarnings(suppressMessages(DI(y = y, block = block, density = density, 
-                       prop = prop, treat = treat, estimate_theta = F, 
+                       prop = prop, treat = treat, ID = ID, estimate_theta = F, 
                        data = data, DImodel = 'FULL')))
       model_list[[length(model_list)+1]] <- FULL_model
       names(model_list)[length(model_list)] <- 'FULL_model'
@@ -537,13 +537,14 @@ autoDI_step2 <- function(y, block, density,
   estimate_theta <- theta_flag
   DImodel <-  fit_final$DIcall$DImodel
   the_final_model <- suppressWarnings(suppressMessages(DI(y = y, prop = prop, block = block,
-                                         FG = FG, density = density, treat = treat,
+                                         FG = FG, ID = ID, density = density, treat = treat,
                                          estimate_theta = estimate_theta,
                                          DImodel = DImodel, data = data)))
   the_final_model$DIcall$prop <- prop
   the_final_model$DIcall$DImodel <- DImodel
   the_final_model$DIcall$estimate_theta <- estimate_theta
   the_final_model$DIcall$FG <- FG
+  the_final_model$DIcall$ID <- ID
   the_final_model$DIcall$treat <- treat
   the_final_model$DIcall$block <- block
   the_final_model$DIcall$density <- density