- translated remainer of .cpp file to R

- replace dependencies on TMB w RTMB
nwfsc-cb · Mar 28, 2024 · 70ae330 · 70ae330
1 parent ff0b21c
commit 70ae330
Show file tree

Hide file tree

Showing 11 changed files with 288 additions and 584 deletions.
diff --git a/.DS_Store b/.DS_Store
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -27,7 +27,7 @@ Imports:
     gnorm,
     methods,
     stats,
-    TMB (>= 1.7.20),
+    RTMB,
     nlme
 Suggests:
     testthat,

diff --git a/NAMESPACE b/NAMESPACE
@@ -25,7 +25,9 @@ export(qdgnorm)
 export(qdnorm)
 export(qdt)
 export(qthill)
+import(RTMB)
 import(ggplot2)
+import(gnorm)
 importFrom(TMB,MakeADFun)
 importFrom(TMB,sdreport)
 importFrom(dplyr,left_join)
@@ -44,4 +46,3 @@ importFrom(stats,qnorm)
 importFrom(stats,qt)
 importFrom(stats,rnorm)
 importFrom(stats,runif)
-useDynLib(phenomix, .registration = TRUE)
diff --git a/R/distributions.R b/R/distributions.R
@@ -18,7 +18,7 @@
 #' @importFrom stats qt
 #' @export
 qthill <- function(p, mean, sd, df) {
-  qt(p, df) * std_dev + mean
+  qt(p, df) * sd + mean
 }
 
 #' Calculate the density of a double normal distribution

diff --git a/R/fit.R b/R/fit.R
@@ -1,6 +1,3 @@
-#' @useDynLib phenomix, .registration = TRUE
-NULL
-
 #' Fitting function to be called by user
 #'
 #' This function creates a list of parameters, sets up TMB object and attempts to
@@ -15,34 +12,36 @@ NULL
 #' @param fit_model Whether to fit the model. If not, returns a list including the data, parameters, and initial values. Defaults to TRUE
 #' @importFrom stats runif rnorm
 #' @importFrom methods is
+#' @import RTMB
+#' @import gnorm
 #' @export
 #' @examples
 #' data(fishdist)
 #'
 #' # example of fitting fixed effects, no trends, no random effects
-#' set.seed(1)
-#' datalist <- create_data(fishdist[which(fishdist$year > 1970), ],
-#'   asymmetric_model = FALSE,
-#'   est_mu_re = FALSE, est_sigma_re = FALSE
-#' )
-#' fit <- fit(datalist)
+# set.seed(2)
+# datalist <- create_data(fishdist[which(fishdist$year > 1970), ],
+#   asymmetric_model = FALSE,
+#   est_mu_re = FALSE, est_sigma_re = FALSE
+# )
+# fit <- fit(datalist, silent=TRUE)
 #' #
 #' # # example of model with random effects in means only, and symmetric distribution
-#' # set.seed(1)
-#' # datalist <- create_data(fishdist[which(fishdist$year > 1970),], asymmetric_model = FALSE,
-#' #                      est_sigma_re = FALSE)
-#' # fit <- fit(datalist)
+# set.seed(1)
+# datalist <- create_data(fishdist[which(fishdist$year > 1970),], asymmetric_model = FALSE,
+#                      est_sigma_re = FALSE)
+# fit <- fit(datalist)
 #' # # example of model with random effects in variances
-#' # set.seed(1)
-#' # datalist <- create_data(fishdist[which(fishdist$year > 1970),], asymmetric_model = TRUE,
-#' #                          est_mu_re = TRUE)
-#' # fit <- fit(datalist)
+# set.seed(1)
+# datalist <- create_data(fishdist[which(fishdist$year > 1970),], asymmetric_model = TRUE,
+#                          est_mu_re = TRUE)
+# fit <- fit(datalist)
 #' #
 #' # # example of model with poisson response
-#' # set.seed(1)
-#' # datalist <- create_data(fishdist[which(fishdist$year > 1970),], asymmetric_model = FALSE,
-#' #                         est_sigma_trend=FALSE, est_mu_trend=FALSE, est_mu_re = TRUE,
-#' #                         family="poisson")
+# set.seed(1)
+# datalist <- create_data(fishdist[which(fishdist$year > 1970),], asymmetric_model = FALSE,
+#                         est_sigma_trend=FALSE, est_mu_trend=FALSE, est_mu_re = TRUE,
+#                         family="poisson")
 #' # fit <- fit(datalist)
 fit <- function(data_list,
                 silent = FALSE,
@@ -195,11 +194,271 @@ fit <- function(data_list,
   random <- random[-1]
   if (length(random) == 0) random <- NULL
 
-  obj <- TMB::MakeADFun(
-    data = data_list,
+  # perhaps it'd be better to split this function out?
+  f <- function(parms) {
+
+    RTMB::getAll(data_list, parms, warn=FALSE)
+    nll<-0
+    y <- RTMB::OBS(y) # optional, https://cran.r-project.org/web/packages/RTMB/vignettes/RTMB-introduction.html
+    # derived parameters
+    obs_sigma <- RTMB:::getValues(exp(log_obs_sigma))
+    tdf_1 <- exp(log_tdf_1) + 2
+    tdf_2 <- exp(log_tdf_2) + 2
+    beta_1 <- exp(log_beta_1)
+    beta_2 <- exp(log_beta_2)
+    if(share_shape==1) {
+      tdf_2 <- tdf_1
+      beta_2 <- beta_1
+    }
+
+    if(use_t_prior & tail_model == 1) {
+      #https:#statmodeling.stat.columbia.edu/2015/05/17/do-we-have-any-recommendations-for-priors-for-student_ts-degrees-of-freedom-parameter/
+        nll <- nll + dgamma(RTMB:::getValues(tdf_1), nu_prior[1], nu_prior[2], log=TRUE)
+        if(asymmetric == 1) nll <- nll + dgamma(RTMB:::getValues(tdf_2), nu_prior[1], nu_prior[2], log=TRUE)
+    }
+    if(use_beta_prior & tail_model == 2) {
+      # value ~ 2-3 is normal
+      nll <- nll + dgamma(RTMB:::getValues(beta_1), beta_prior[1], beta_prior[2], log=TRUE)
+      if(asymmetric == 1) nll <- nll + dgamma(RTMB:::getValues(beta_2), beta_prior[1], beta_prior[2], log=TRUE)
+    }
+
+    sigma1 <- sigma2 <- mu <- alpha1 <- alpha2 <- lower25 <- upper75 <- range <- rep(NA, nLevels)
+    n <- length(y)
+
+    # calculations for beta for gnorm dist if implemented
+    beta_ratio <- rep(NA,2)
+    if(tail_model == 2) {
+      beta_ratio[1] <- sqrt(exp(lgamma(1.0/RTMB:::getValues(beta_1))) / exp(lgamma(3.0/RTMB:::getValues(beta_1))))
+      if(asymmetric == 1) {
+        beta_ratio[2] <- sqrt(exp(lgamma(1.0/RTMB:::getValues(beta_2))) / exp(lgamma(3.0/RTMB:::getValues(beta_2))))
+      }
+    }
+
+    # random effects components
+    for(i in 1:nLevels) {
+      # random effects contributions of mean and sigma1
+      if(est_mu_re==1) {
+        nll <- nll + dnorm(RTMB:::getValues(mu_devs[i]), 0, exp(RTMB:::getValues(log_sigma_mu_devs)), log=TRUE)
+      }
+      if(est_sigma_re==1) {
+        nll <- dnorm(RTMB:::getValues(sigma1_devs[i]), 0, exp(RTMB:::getValues(log_sigma1_sd)), log=TRUE)
+        if(asymmetric == 1) {
+          nll <- dnorm(RTMB:::getValues(sigma2_devs[i]), 0, exp(RTMB:::getValues(log_sigma2_sd)), log=TRUE)
+        }
+      }
+    }
+
+    # fixed effects for mu and both sigmas
+    mu <- mu_mat * RTMB:::getValues(b_mu)
+    sigma1 <- sig_mat * RTMB:::getValues(b_sig1)
+    if(asymmetric == 1) sigma2 <- sig_mat * RTMB:::getValues(b_sig2)
+
+    for(i in 1:nLevels) {
+
+      if(est_sigma_re == 1) {
+        sigma1[i] <- sigma1[i] + RTMB:::getValues(sigma1_devs[i])
+        if(asymmetric == 1) sigma2[i] <- sigma2[i] + RTMB:::getValues(sigma2_devs[i])
+      }
+
+      # calculate alphas if the gnorm model is used
+      if(tail_model == 2) {
+        alpha1[i] <- sigma1[years[i]]*beta_ratio[1]
+        if(asymmetric==1) {
+          alpha2[i] <- sigma2[years[i]]*beta_ratio[2]
+        }
+      }
+
+      # trend in in normal space, e.g. not log-linear
+      if(est_mu_re == 1) {
+        mu[i] <- mu[i] + RTMB:::getValues(mu_devs[i])
+      }
+
+      # this is all for calculating quantiles of normal distributions
+      if(tail_model==0) {
+        if(asymmetric == 0) {
+          lower25[i] <- qnorm(0.25, mu[i], sigma1[i])
+          upper75[i] <- qnorm(0.75, mu[i], sigma1[i])
+        } else {
+          lower25[i] <- qdnorm(0.25, mu[i], sigma1[i], sigma2[i])
+          upper75[i] <- qdnorm(0.75, mu[i], sigma1[i], sigma2[i])
+        }
+      }
+      # this is all for calculating quantiles of Student-t distributions
+      if(tail_model==1) {
+        if(asymmetric == 0) {
+          lower25[i] <- qthill(0.25, RTMB:::getValues(tdf_1), mu[i], sigma1[i])
+          upper75[i] <- qthill(0.75, RTMB:::getValues(tdf_1), mu[i], sigma1[i])
+        } else {
+          lower25[i] <- qdt(0.25, mu[i], sigma1[i], sigma2[i], RTMB:::getValues(tdf_1), RTMB:::getValues(tdf_2))
+          upper75[i] <- qdt(0.75, mu[i], sigma1[i], sigma2[i], RTMB:::getValues(tdf_1), RTMB:::getValues(tdf_2))
+        }
+      }
+
+      # this is all for calculating quantiles of Student-t distributions
+      if(tail_model==2) {
+        if(asymmetric == 0) {
+          lower25[i] <- qgnorm(0.25, mu[i], sigma1[i]*beta_ratio[1], RTMB:::getValues(beta_1))
+          upper75[i] <- qgnorm(0.75, mu[i], sigma1[i]*beta_ratio[1], RTMB:::getValues(beta_1))
+        } else {
+          lower25[i] <- qdgnorm(0.25, mu[i], sigma1[i], sigma2[i], beta_ratio[1], beta_ratio[2], RTMB:::getValues(beta_1), RTMB:::getValues(beta_2))
+          upper75[i] <- qdgnorm(0.75, mu[i], sigma1[i], sigma2[i], beta_ratio[1], beta_ratio[2], RTMB:::getValues(beta_1), RTMB:::getValues(beta_2))
+        }
+      }
+      range[i] = upper75[i] - lower25[i]
+    }
+
+    # this is for the predictions
+    log_dens <- pred <- rep(NA, n)
+    for(i in 1:n) {
+      if(asymmetric == 1) {
+        # model is asymmetric, left side smaller / right side bigger
+        if(tail_model==0) {
+          log_dens[i] <- ddnorm(x[i], mu[years[i]], sigma1[years[i]], sigma2[years[i]])
+        }
+        if(tail_model==1) {
+          log_dens[i] <- ddt(x[i], mu[years[i]], sigma1[years[i]], sigma2[years[i]], RTMB:::getValues(tdf_1), RTMB:::getValues(tdf_2))
+        }
+        if(tail_model==2) {
+          log_dens[i] <- ddgnorm(x[i], mu[years[i]], alpha1[years[i]], alpha2[years[i]], RTMB:::getValues(beta_1), RTMB:::getValues(beta_2), sigma1[years[i]], sigma2[years[i]])
+        }
+        pred[i] <- log_dens[i] + RTMB:::getValues(theta[years[i]])
+
+      } else {
+        if(tail_model==0) {
+          # model is symmetric around mu, gaussian tails
+          log_dens[i] <- dnorm(x[i], mu[years[i]], sigma1[years[i]], log=TRUE)
+        } else {
+          if(tail_model==1) {
+            # model is symmetric around mu, student-t tails
+            log_dens[i] <- dt((x[i] - mu[years[i]]) / sigma1[years[i]], RTMB:::getValues(tdf_1), log=TRUE) - log(sigma1[years[i]])
+          } else {
+            # gnorm, copied from maryclare/gnorm
+            # alpha = sqrt( var * gamma(1/beta) / gamma(3/beta) ), alpha = sigma(1)*beta_ratio(1)
+            log_dens[i] <- dgnorm(x[i], mu[years[i]], alpha1[years[i]], RTMB:::getValues(beta_1))
+          }
+        }
+        pred[i] <- log_dens[i] + RTMB:::getValues(theta[years[i]])
+      }
+    }
+
+    # this is for the cumulative annual predictions
+    year_log_tot <- rep(0, nLevels)
+    year_tot <- rep(0, nLevels)
+    for(i in 1:nLevels) {
+
+      for(t in 1:366) {
+
+        if(asymmetric == 1) {
+          # model is asymmetric, left side smaller / right side bigger
+          if(tail_model==0) {
+            dens <- ddnorm(t, mu[i], sigma1[i], sigma2[i])
+          }
+          if(tail_model==1) {
+            dens <- ddt(t, mu[i], sigma1[i], sigma2[i], RTMB:::getValues(tdf_1), RTMB:::getValues(tdf_2))
+          }
+          if(tail_model==2) {
+            dens <- ddgnorm(t, mu[i], alpha1[i], alpha2[i], RTMB:::getValues(beta_1), RTMB:::getValues(beta_2), sigma1[i], sigma2[i])
+          }
+          year_log_tot[i] <- year_log_tot[i] + dens + RTMB:::getValues(theta[i])
+          year_tot[i] <- year_tot[i] + exp(dens + RTMB:::getValues(theta[i]))
+
+        } else {
+          if(tail_model==0) {
+            # model is symmetric around mu, gaussian tails
+            dens <- dnorm(t, mu[i], sigma1[i], log=TRUE)
+          } else {
+            if(tail_model==1) {
+              # model is symmetric around mu, student-t tails
+              dens <- dt((t - mu[i]) / sigma1[i], RTMB:::getValues(tdf_1), log=TRUE) - log(sigma1[i])
+            } else {
+              # gnorm, copied from maryclare/gnorm
+              # alpha = sqrt( var * gamma(1/beta) / gamma(3/beta) ), alpha = sigma(1)*beta_ratio(1)
+              dens <- dgnorm(t, mu[i], alpha1[i], RTMB:::getValues(beta_1))
+            }
+          }
+          year_log_tot[i] <- year_log_tot[i] + dens + RTMB:::getValues(theta[i])
+          year_tot[i] <- year_tot[i] + exp(dens + RTMB:::getValues(theta[i]))
+        }
+      }
+    }
+
+    # this is the likelihood
+    s1 <- 0
+    s2 <- 0
+
+    if(family==1) {
+      # gaussian, both data and predictions in log space
+      nll <- nll + sum(dnorm(y, pred, obs_sigma, log=TRUE))
+    }
+    if(family==2) {
+      for(i in 1:n) {
+        if(pred[i] > 20) pred[i] <- 20 # not needed
+        nll <- dpois(y[i], exp(pred[i]), log=TRUE)
+      }
+    }
+    if(family==3) {
+      for(i in 1:n) {
+        s1 = pred[i]
+        #s2 = s1 + pow(s1, Type(2))*obs_sigma
+        s2 <- 2 * s1 - RTMB:::getValues(log_obs_sigma) # log(var - mu)
+        nll <- nll + RTMB::dnbinom_robust(y[i], s1, s2, log=TRUE)
+      }
+    }
+    if(family==4) {
+      for(i in 1:n) {
+        nll <- nll + RTMB::dbinom_robust(y[i], 1, pred[i], log=TRUE)
+      }
+    }
+    if(family==5) {
+      # lognormal, both data and predictions in log space
+      nll <- sum(dnorm(log(y), pred, obs_sigma, log=TRUE))
+    }
+    # ADREPORT section
+    RTMB::ADREPORT(theta) # nuisance parameter
+    RTMB::ADREPORT(sigma1) # sigma, LHS
+    RTMB::ADREPORT(mu) # mean of curves by year
+    RTMB::ADREPORT(b_mu) # sigma, LHS
+    RTMB::ADREPORT(b_sig1) # mean of curves by year
+    RTMB::ADREPORT(year_tot) # mean of curves by year
+    RTMB::ADREPORT(year_log_tot) # mean of curves by year
+
+    if(family %in% c(2,4) == FALSE) {
+      RTMB::ADREPORT(obs_sigma) # obs sd (or phi, NB)
+    }
+    RTMB::ADREPORT(pred) # predictions in link space (log)
+
+    RTMB::ADREPORT(lower25) # lower quartile
+    RTMB::ADREPORT(upper75) # upper quartile
+    RTMB::ADREPORT(range) # diff between upper and lower quartiles
+    if(tail_model==1) {
+      RTMB::ADREPORT(tdf_1) # tdf for LHS
+    }
+    if(tail_model==2) {
+      RTMB::ADREPORT(beta_1) # tdf for LHS
+    }
+    if(asymmetric == 1) {
+      # these are only reported for asymmetric model
+      RTMB::ADREPORT(b_sig2) # mean of curves by year
+      if(est_sigma_re==1) {
+        RTMB::ADREPORT(sigma2) # same as above, but RHS optionally
+      }
+      if(tail_model==1) {
+        RTMB::ADREPORT(tdf_2)
+      }
+      if(tail_model==2) {
+        RTMB::ADREPORT(beta_2)
+      }
+    }
+    return (-nll)
+
+  }
+
+  obj <- RTMB::MakeADFun(
+    func = f,
+    #data = data_list,
     parameters = parameters,
     map = tmb_map,
-    DLL = "phenomix",
+    #DLL = "phenomix",
     random = random,
     silent = silent
   )
@@ -244,7 +503,7 @@ fit <- function(data_list,
       )
     }
 
-    sdreport <- TMB::sdreport(obj)
+    sdreport <- RTMB::sdreport(obj)
     mod_list <- c(mod_list, list(pars = pars, sdreport = sdreport, tmb_map = tmb_map, tmb_random = random, tmb_parameters = parameters))
   }
 

diff --git a/R/phenomix-package.R b/R/phenomix-package.R