R Under development (unstable) (2024-11-15 r87338 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2024 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> 
> ##
> ## A Monte Carlo study to test measurement error models
> ##
> 
> library(geostan)
This is geostan version 0.8.0
> 
> ## no. iterations
> M = 15
> 
> ## a regular grid
> ncol = 20
> nrow = 20
> sars <- prep_sar_data2(ncol, nrow, quiet = TRUE)
> cars <- prep_car_data2(ncol, nrow, quiet = TRUE)
> W <- sars$W
> N <- nrow(W)
> 
> # iterate:
> b = 0.5
> g = -0.5
> rho.x = 0 ## checking ME, not checking spatial models
> rho.y = 0.6
> sigma.y = 0.3
> sigma.me = 0.3
> pars <- c(const = 0,
+           gamma = g,
+           beta = b,
+           rho = rho.y,
+           scale = sigma.y)
> 
> res <- sapply(1:M, FUN = function(i) {
+ 
+     x <- sim_sar(w=W, rho=rho.x)
+     Wx <- (W %*% x)[,1]    
+     mu <- b * x + g * Wx    
+     y <- sim_sar(w=W, rho=rho.y, mu = mu, sigma = sigma.y, type = "SEM")
+     x = x + rnorm(N, sd = sigma.me)        
+     dat <- data.frame(y, x)    
+ 
+     ME <- prep_me_data(se = data.frame(x = rep(sigma.me, N)))
+     if (rho.x > 0) ME <- prep_me_data(se = data.frame(x = rep(sigma.me, N)), car_parts = cars)
+ 
+     fit_sem_me <- geostan::stan_sar(y ~ x,
+                                  data = dat,
+                                  type = "SDEM",
+                                  sar_parts = sars,
+                                  ME = ME,
+                                  chains = 1,
+                                  iter = 600,
+                                  quiet = TRUE,
+                                  slim = TRUE
+                                  ) |>
+         suppressWarnings()
+ 
+     fit_sem <- geostan::stan_sar(y ~ x,
+                                   data = dat,
+                                   type = "SDEM",
+                                   sar_parts = sars,
+                                   ## ME = ME, ##
+                                   chains = 1,
+                                   iter = 600,
+                                   quiet = TRUE,
+                                   slim = TRUE
+                                   ) |>
+         suppressWarnings()
+     ## fit_glm0 <- geostan::stan_glm(y ~ x,
+         ##                          slx = ~ x,
+         ##                          data = dat,
+         ##                          ## ME = ME, ##
+         ##                          C = cars$C,
+         ##                          chains = 1,
+         ##                          iter = 600,
+         ##                          quiet = TRUE,
+         ##                          slim = TRUE
+         ##                          ) |>
+         ## suppressWarnings()
+     ## fit_car <- geostan::stan_car(y ~ x,
+     ##                              slx = ~ x,
+     ##                              data = dat,
+     ##                              car_parts = cars,
+     ##                              ME = ME,
+     ##                              chains = 1,
+     ##                              iter = 900,
+     ##                              quiet = TRUE,
+     ##                              slim = TRUE
+     ##                              ) |>
+     ##     suppressWarnings()
+ 
+     ## fit_glm <- geostan::stan_glm(y ~ x,
+     ##                              slx = ~ x,
+     ##                              data = dat,
+     ##                              ME = ME,
+     ##                              C = cars$C,
+     ##                              chains = 1,
+     ##                              iter = 900,
+     ##                              quiet = TRUE,
+     ##                              slim = TRUE
+     ##                              ) |>
+     ##     suppressWarnings()
+     
+     x <- c(
+         # GLM0 = fit_glm0$summary[c('intercept', 'w.x', 'x', 'rho', 'sigma'), 'mean'],
+        # GLM = fit_glm$summary[c('intercept', 'w.x', 'x', 'rho', 'sigma'), 'mean'],        
+         SEM_ME = fit_sem_me$summary[c('intercept', 'w.x', 'x', 'sar_rho', 'sar_scale'), 'mean'],
+         SEM = fit_sem$summary[c('intercept', 'w.x', 'x', 'sar_rho', 'sar_scale'), 'mean'] #,        
+        # CAR = fit_car$summary[c('intercept', 'w.x', 'x', 'car_rho', 'car_scale'), 'mean']
+         )
+     
+     return (x)
+     
+     })
> 
> 
> RMSE <- function(est, true) sqrt(mean(est - true)^2)
> 
> g_res <- res |>
+     subset(row.names(res) %in% c('SEM_ME2', 'SEM2'))
>    ## subset(row.names(res) %in% c('GLM02', 'GLM2', 'SEM2', 'CAR2'))
> b_res <- res |>
+     subset(row.names(res) %in% c('SEM_ME3', 'SEM3'))    
>      ##  subset(row.names(res) %in% c('GLM03', 'GLM2', 'SEM3', 'CAR3'))
> 
> cat("\n**\nRMSE of ME models: \n**\n")

**
RMSE of ME models: 
**
> 
> cat("\n**\nGamma\n**\n")

**
Gamma
**
> apply(g_res, 1, RMSE, g) |>
+     print(digits = 2)
SEM_ME2    SEM2 
 0.0057  0.0374 
> 
> cat("\n**\nBeta\n**\n")

**
Beta
**
> apply(b_res, 1, RMSE, b) |>
+     print(digits = 2)
SEM_ME3    SEM3 
 0.0052  0.0373 
> 
> 
> cat("\n**\nME models\nSEM Estimates (rounded) should be close to DGP parameters \n",M, "iterations\n**\n")

**
ME models
SEM Estimates (rounded) should be close to DGP parameters 
 15 iterations
**
> 
> est <- apply(res, 1, mean) 
> 
> est <- est[c('SEM2', 'SEM3', 'SEM_ME2', 'SEM_ME3')]
> 
> out <- data.frame(
+    # Mod = attributes(est)$names,
+     Par = rep(c('Gamma', 'Beta'), 2),
+     DGP = c(g, b, g, b),
+     Est = est
+ ) |>
+     transform(
+         Error = round(Est - DGP, 2)
+     ) 
> 
> out |>
+     print()
          Par  DGP        Est Error
SEM2    Gamma -0.5 -0.4625611  0.04
SEM3     Beta  0.5  0.4627315 -0.04
SEM_ME2 Gamma -0.5 -0.5056960 -0.01
SEM_ME3  Beta  0.5  0.5051826  0.01
> 
> write.csv(out, "check-ME-estimates-monte-carlo-output.csv")
> 
> proc.time()
   user  system elapsed 
  69.42    0.50   69.90