fix dsemRTMB to work without loading locally

NAMESPACE

@@ -24,7 +24,9 @@ export(make_dsem_ram)
 export(parse_path)
 export(stepwise_selection)
 importFrom(Matrix,Cholesky)
+importFrom(Matrix,Diagonal)
 importFrom(Matrix,drop0)
+importFrom(Matrix,kronecker)
 importFrom(Matrix,mat2triplet)
 importFrom(Matrix,solve)
 importFrom(Matrix,sparseMatrix)

R/dsemRTMB.R

@@ -11,6 +11,8 @@
 #'        specifies a normal prior probability for the first path coefficient
 #'        with mean of zero and sd of 0.1
 #'
+#' @importFrom Matrix solve Diagonal sparseMatrix drop0 kronecker
+#'
 #' @details
 #' See \code{\link{dsem}} for details
 #'

R/make_matrices.R

@@ -39,17 +39,17 @@ function( beta_p,
   I_kk = Diagonal(nrow(P_kk))
 
   # Assemble
-  IminusP_kk = I_kk - P_kk
+  IminusP_kk = AD(I_kk - P_kk)
   invV_kk = AD(G_kk)
   invV_kk@x = 1 / G_kk@x^2
-  Q_kk = t(IminusP_kk) %*% invV_kk %*% IminusP_kk
+  #Q_kk = t(IminusP_kk) %*% invV_kk %*% IminusP_kk
 
   out = list(
     "P_kk" = P_kk,
     "G_kk" = G_kk,
     "invV_kk" = invV_kk,
-    "IminusP_kk" = IminusP_kk,
-    "Q_kk" = Q_kk
+    #"Q_kk" = Q_kk,     # NOT USED
+    "IminusP_kk" = IminusP_kk
   )
   return(out)
 }

scratch/test_dsemRTMB.R

@@ -2,8 +2,8 @@
 #devtools::install_local( R'(C:\Users\James.Thorson\Desktop\Git\dsem)', force=TRUE, dep=FALSE )
 
 library(dsem)
-library(RTMB)
-library(Matrix)
+#library(RTMB)
+#library(Matrix)
 
 # Define model
 sem = "
@@ -74,14 +74,59 @@ fit$obj$simulate()$y_tj
 # dsemRTMB
 ###################
 
+#devtools::install_local( R'(C:\Users\James.Thorson\Desktop\Git\dsem)', force=TRUE, dep=FALSE )
+
+library(dsem)
+
 #devtools::install_github( "kaskr/RTMB/RTMB", force=TRUE, dep=FALSE )
 
 # Files
-source( file.path(R'(C:\Users\James.Thorson\Desktop\Git\dsem\R)', "make_matrices.R") )
-source( file.path(R'(C:\Users\James.Thorson\Desktop\Git\dsem\R)', "compute_nll.R") )
-source( file.path(R'(C:\Users\James.Thorson\Desktop\Git\dsem\R)', "read_model.R") )
-source( file.path(R'(C:\Users\James.Thorson\Desktop\Git\dsem\R)', "dsemRTMB.R") )
-source( file.path(R'(C:\Users\James.Thorson\Desktop\Git\dsem\R)', "rgmrf.R") )
+if( FALSE ){
+  source( file.path(R'(C:\Users\James.Thorson\Desktop\Git\dsem\R)', "make_matrices.R") )
+  source( file.path(R'(C:\Users\James.Thorson\Desktop\Git\dsem\R)', "compute_nll.R") )
+  source( file.path(R'(C:\Users\James.Thorson\Desktop\Git\dsem\R)', "read_model.R") )
+  source( file.path(R'(C:\Users\James.Thorson\Desktop\Git\dsem\R)', "dsemRTMB.R") )
+  source( file.path(R'(C:\Users\James.Thorson\Desktop\Git\dsem\R)', "rgmrf.R") )
+}
+
+# Define model
+sem = "
+  # Link, lag, param_name
+  cprofits -> consumption, 0, a1
+  cprofits -> consumption, 1, a2
+  pwage -> consumption, 0, a3
+  gwage -> consumption, 0, a3
+  cprofits -> invest, 0, b1
+  cprofits -> invest, 1, b2
+  capital -> invest, 0, b3
+  gnp -> pwage, 0, c2
+  gnp -> pwage, 1, c3
+  time -> pwage, 0, c1
+"
+
+# Load data
+data(KleinI, package="AER")
+TS = ts(data.frame(KleinI, "time"=time(KleinI) - 1931))
+tsdata = TS[,c("time","gnp","pwage","cprofits",'consumption',
+               "gwage","invest","capital")]
+
+# Fit model
+fit0 = dsemRTMB( sem=sem,
+            tsdata = tsdata,
+            estimate_delta0 = TRUE,
+            family = rep("normal",ncol(tsdata)),
+            control = dsem_control( quiet=TRUE,
+                                    run_model = FALSE,
+                                    use_REML = TRUE,
+                                    gmrf_parameterization = "projection" ) )
+#ParHat = fit$internal$parhat
+#Rep = fit$obj$report()
+
+#
+Map = fit0$tmb_inputs$map
+Map$lnsigma_j = factor( rep(NA,ncol(tsdata)) )
+Params = fit0$tmb_inputs$parameters
+Params$lnsigma_j[] = log(0.1)
 
 # Define prior
 log_prior = function(p) dnorm( p$beta_z[1], mean=0, sd=0.1, log=TRUE)
@@ -95,7 +140,7 @@ fitRTMB = dsemRTMB( sem = sem,
                                     run_model = TRUE,
                                     use_REML = TRUE,
                                     gmrf_parameterization = "projection",
-                                    #constant_variance = c("conditional", "marginal", "diagonal")[2], # marginal not working
+                                    constant_variance = c("conditional", "marginal", "diagonal")[2], # marginal not working
                                     map = Map,
                                     parameters = Params ) )
 #Rep = fitRTMB$obj$report()