# function to parallelize my_func <- function(ind_b) { B <- length(ind_b) column_0a <- matrix(0, 16, B) column_0b <- matrix(0, 16, B) column_1 <- matrix(0, 16, B) column_2 <- matrix(0, 16, B) column_3 <- matrix(0, 16, B) column_4 <- matrix(0, 16, B) column_5 <- matrix(0, 16, B) column_6 <- matrix(0, 16, B) column_7 <- matrix(0, 16, B) for(k in 1:B) { set.seed(ind_b[k]) # column 0 epsi <- MASS::mvrnorm(283, c(0, 0), Sigma = sigma_bar_d) x_beta_epsi <- x_simu %*% beta_true + epsi y_lm_1 <- lm(x_beta_epsi[, 1] ~ x_simu - 1) y_lm_2 <- lm(x_beta_epsi[, 2] ~ x_simu - 1) resi <- cbind(residuals(y_lm_1), residuals(y_lm_2)) column_0a[, k] <- c(as.vector(y_lm_1$coefficients), as.vector(y_lm_2$coefficients), 0, 0, 0, 0, as.vector(t(resi) %*% resi / 279)) # column 0 epsi <- MASS::mvrnorm(283, c(0, 0), Sigma = sigma_d) x_beta_epsi <- x_simu %*% beta_true + epsi y_lm_1 <- lm(x_beta_epsi[, 1] ~ x_simu - 1) y_lm_2 <- lm(x_beta_epsi[, 2] ~ x_simu - 1) column_0b[, k] <- c(as.vector(y_lm_1$coefficients), as.vector(y_lm_2$coefficients), 0, 0, 0, 0, sum(residuals(y_lm_1)^2)/ 279, 0, 0, sum(residuals(y_lm_2)^2)/ 279) # column 1 y_N_mod_1 <- simu_spatial_multi_y(X = x_simu, beta_true = beta_true, method_simulate = "N", Sigma = sigma_bar_d, RHO = RHO_bar_d, W = W_simu) res <- estimate_spatial_multi_gen_N(Y = y_N_mod_1, X = x_simu, W = W_simu, ind_beta = matrix(c(T, T, T, T, T, T, T, T, T, T, T, T), 4, 3), ind_RHO = matrix(c(T, T, T, T), 2, 2), ind_GAMMA = matrix(c(F, F, F, F), 2, 2), method = "s2sls") column_1[, k] <- c(as.vector(res$est_beta), as.vector(res$est_RHO), as.vector(res$est_SIGMA)) # column 2 y_N_mod_2 <- simu_spatial_multi_y(X = x_simu, beta_true = beta_true, method_simulate = "N", Sigma = sigma_d, RHO = RHO_bar_d, W = W_simu) res <- estimate_spatial_multi_gen_N(Y = y_N_mod_2, X = x_simu, W = W_simu, ind_beta = matrix(c(T, T, T, T, T, T, T, T, T, T, T, T), 4, 3), ind_RHO = matrix(c(T, T, T, T), 2, 2), ind_GAMMA = matrix(c(F, F, F, F), 2, 2), method = "s2sls") column_2[, k] <- c(as.vector(res$est_beta), as.vector(res$est_RHO), as.vector(res$est_SIGMA)) # column 3 y_N_mod_3 <- simu_spatial_multi_y(X = x_simu, beta_true = beta_true, method_simulate = "N", Sigma = sigma_bar_d, RHO = RHO_d, W = W_simu) res <- estimate_spatial_multi_gen_N(Y = y_N_mod_3, X = x_simu, W = W_simu, ind_beta = matrix(c(T, T, T, T, T, T, T, T, T, T, T, T), 4, 3), ind_RHO = matrix(c(T, F, F, T), 2, 2), ind_GAMMA = matrix(c(F, F, F, F), 2, 2), method = "s2sls") column_3[, k] <- c(as.vector(res$est_beta), as.vector(res$est_RHO), as.vector(res$est_SIGMA)) # column 4 res <- estimate_spatial_multi_gen_N(Y = y_N_mod_3, X = x_simu, W = W_simu, ind_beta = matrix(c(T, T, T, T, T, T, T, T, T, T, T, T), 4, 3), ind_RHO = matrix(c(T, F, F, T), 2, 2), ind_GAMMA = matrix(c(F, F, F, F), 2, 2), method = "s3sls") column_4[, k] <- c(as.vector(res$est_beta), as.vector(res$est_RHO), as.vector(res$est_SIGMA)) # column 5 y_N_mod_4 <- simu_spatial_multi_y(X = x_simu, beta_true = beta_true, method_simulate = "N", Sigma = sigma_d, RHO = RHO_d, W = W_simu) res <- estimate_spatial_multi_gen_N(Y = y_N_mod_4, X = x_simu, W = W_simu, ind_beta = matrix(c(T, T, T, T, T, T, T, T, T, T, T, T), 4, 3), ind_RHO = matrix(c(T, F, F, T), 2, 2), ind_GAMMA = matrix(c(F, F, F, F), 2, 2), method = "s2sls") column_5[, k] <- c(as.vector(res$est_beta), as.vector(res$est_RHO), as.vector(res$est_SIGMA)) # column 6 res <- estimate_spatial_multi_gen_N(Y = y_N_mod_4, X = x_simu, W = W_simu, ind_beta = matrix(c(T, T, T, T, T, T, T, T, T, T, T, T), 4, 3), ind_RHO = matrix(c(T, F, F, T), 2, 2), ind_GAMMA = matrix(c(F, F, F, F), 2, 2), method = "s3sls") column_6[, k] <- c(as.vector(res$est_beta), as.vector(res$est_RHO), as.vector(res$est_SIGMA)) # column 7 res_a <- spatialreg::lagsarlm(y_N_mod_4[, 1] ~ x_simu - 1, listw = W1.listw) res_b <- spatialreg::lagsarlm(y_N_mod_4[, 2] ~ x_simu - 1, listw = W1.listw) column_7[, k] <- c(as.vector(res_a$coefficients), as.vector(res_b$coefficients), res_a$rho, 0, 0, res_b$rho, res_a$s2, 0, 0, res_b$s2) } return(list(column_0a = column_0a, column_0b = column_0b, column_1 = column_1, column_2 = column_2, column_3 = column_3, column_4 = column_4, column_5 = column_5, column_6 = column_6, column_7 = column_7)) } ############################################################################## ############################ Covariance matrix presented in the article ############################################################################## B <- 1000 P <- 4 library(parallel) cl <- parallel::makeCluster(P) clusterExport(cl, c("x_simu", "beta_true", "sigma_bar_d", "sigma_d", "RHO_bar_d", "RHO_d", "W_simu", "simu_spatial_multi_y", "estimate_spatial_multi_gen_N", "W1.listw")) clusterEvalQ(cl, library("MASS")) cut_B <- split(1:1000, sort(rank(1:B) %% P)) system.time( res_par <- clusterApply(cl, cut_B, fun = my_func) ) stopCluster(cl) # merge the results column_0a <- res_par[[1]]$column_0a column_0b <- res_par[[1]]$column_0b column_1 <- res_par[[1]]$column_1 column_2 <- res_par[[1]]$column_2 column_3 <- res_par[[1]]$column_3 column_4 <- res_par[[1]]$column_4 column_5 <- res_par[[1]]$column_5 column_6 <- res_par[[1]]$column_6 column_7 <- res_par[[1]]$column_7 for (k in 2:P) { column_0a <- cbind(column_0a, res_par[[k]]$column_0a) column_0b <- cbind(column_0b, res_par[[k]]$column_0b) column_1 <- cbind(column_1, res_par[[k]]$column_1) column_2 <- cbind(column_2, res_par[[k]]$column_2) column_3 <- cbind(column_3, res_par[[k]]$column_3) column_4 <- cbind(column_4, res_par[[k]]$column_4) column_5 <- cbind(column_5, res_par[[k]]$column_5) column_6 <- cbind(column_6, res_par[[k]]$column_6) column_7 <- cbind(column_7, res_par[[k]]$column_7) } # compute the RRMMSE # The true parameters in the 4 models true_prms_0a <- c(as.vector(beta_true), c(0, 0, 0, 0), as.vector(sigma_bar_d)) true_prms_0b <- c(as.vector(beta_true), c(0, 0, 0, 0), as.vector(sigma_d)) true_prms_1 <- c(as.vector(beta_true), as.vector(RHO_bar_d), as.vector(sigma_bar_d)) true_prms_2 <- c(as.vector(beta_true), as.vector(RHO_bar_d), as.vector(sigma_d)) true_prms_3 <- c(as.vector(beta_true), as.vector(RHO_d), as.vector(sigma_bar_d)) true_prms_4 <- c(as.vector(beta_true), as.vector(RHO_d), as.vector(sigma_d)) # initialization # biais biais_column_0a <- 0 biais_column_0b <- 0 biais_column_1 <- 0 biais_column_2 <- 0 biais_column_3 <- 0 biais_column_4 <- 0 biais_column_5 <- 0 biais_column_6 <- 0 biais_column_7 <- 0 # rmse rmse_column_0a <- 0 rmse_column_0b <- 0 rmse_column_1 <- 0 rmse_column_2 <- 0 rmse_column_3 <- 0 rmse_column_4 <- 0 rmse_column_5 <- 0 rmse_column_6 <- 0 rmse_column_7 <- 0 # beta_matrix <- matrix(0, B, 8) # RHO_matrix <- matrix(0, B, 4) # sigma_matrix <- matrix(0, B, 4) for (k in 1:B) { # column 0 parms_col_0a <- column_0a[, k] biais_column_0a <- biais_column_0a + (parms_col_0a - true_prms_0a) rmse_column_0a <- rmse_column_0a + (parms_col_0a - true_prms_0a) ^ 2 parms_col_0b <- column_0b[, k] biais_column_0b <- biais_column_0b + (parms_col_0b - true_prms_0b) rmse_column_0b <- rmse_column_0b + (parms_col_0b - true_prms_0b) ^ 2 # column 1 # beta_matrix[k, ] <- as.vector(column_1[[k]]$est_beta) # RHO_matrix[k, ] <- as.vector(column_1[[k]]$est_RHO) # sigma_matrix[k, ] <- as.vector(column_1[[k]]$est_SIGMA) parms_col_1 <- column_1[, k] biais_column_1 <- biais_column_1 + (parms_col_1 - true_prms_1) rmse_column_1 <- rmse_column_1 + (parms_col_1 - true_prms_1) ^ 2 # column 2 parms_col_2 <- column_2[, k] biais_column_2 <- biais_column_2 + (parms_col_2 - true_prms_2) rmse_column_2 <- rmse_column_2 + (parms_col_2 - true_prms_2) ^ 2 # column 3 parms_col_3 <- column_3[, k] biais_column_3 <- biais_column_3 + (parms_col_3 - true_prms_3) rmse_column_3 <- rmse_column_3 + (parms_col_3 - true_prms_3) ^ 2 # column 4 parms_col_4 <- column_4[, k] biais_column_4 <- biais_column_4 + (parms_col_4 - true_prms_3) rmse_column_4 <- rmse_column_4 + (parms_col_4 - true_prms_3) ^ 2 # column 5 parms_col_5 <- column_5[, k] biais_column_5 <- biais_column_5 + (parms_col_5 - true_prms_4) rmse_column_5 <- rmse_column_5 + (parms_col_5 - true_prms_4) ^ 2 # column 6 parms_col_6 <- column_6[, k] biais_column_6 <- biais_column_6 + (parms_col_6 - true_prms_4) rmse_column_6 <- rmse_column_6 + (parms_col_6 - true_prms_4) ^ 2 # column 7 parms_col_7 <- column_7[, k] biais_column_7 <- biais_column_7 + (parms_col_7 - true_prms_4) rmse_column_7 <- rmse_column_7 + (parms_col_7 - true_prms_4) ^ 2 } # Bias (res_biais <- data.frame(col_0a = biais_column_0a / B, col_0b = biais_column_0b / B, col_1 = biais_column_1 / B, col_2 = biais_column_2 / B, col_3 = biais_column_3 / B, col_4 = biais_column_4 / B, col_5 = biais_column_5 / B, col_6 = biais_column_6 / B, col_7 = biais_column_7 / B )) # RRMSE (rrmse_column_0a <- sqrt(rmse_column_0a / B) / abs(true_prms_0a) * 100) (rrmse_column_0b <- sqrt(rmse_column_0b / B) / abs(true_prms_0b) * 100) (rrmse_column_1 <- sqrt(rmse_column_1 / B) / abs(true_prms_1) * 100) (rrmse_column_2 <- sqrt(rmse_column_2 / B) / abs(true_prms_2) * 100) (rrmse_column_3 <- sqrt(rmse_column_3 / B) / abs(true_prms_3) * 100) (rrmse_column_4 <- sqrt(rmse_column_4 / B) / abs(true_prms_3) * 100) (rrmse_column_5 <- sqrt(rmse_column_5 / B) / abs(true_prms_4) * 100) (rrmse_column_6 <- sqrt(rmse_column_6 / B) / abs(true_prms_4) * 100) (rrmse_column_7 <- sqrt(rmse_column_7 / B) / abs(true_prms_4) * 100) (res_table1 <- data.frame(col_0a = rrmse_column_0a, col_0b = rrmse_column_0b, col_1 = rrmse_column_1, col_2 = rrmse_column_2, col_3 = rrmse_column_3, col_4 = rrmse_column_4, col_5 = rrmse_column_5, col_6 = rrmse_column_6, col_7 = rrmse_column_7 )) save(res_biais, res_table1, file = "res_table1.RData") ############################################################################## ############################ Another Covariance matrix ############################################################################## sigma_bar_d <- matrix(c(0.7, 0.09, 0.09, 0.1), nrow = M_simu, ncol = M_simu) sigma_d <- matrix(c(0.7, 0, 0, 0.1), nrow = M_simu, ncol = M_simu) B <- 1000 P <- 4 library(parallel) cl <- parallel::makeCluster(P) clusterExport(cl, c("x_simu", "beta_true", "sigma_bar_d", "sigma_d", "RHO_bar_d", "RHO_d", "W_simu", "simu_spatial_multi_y", "estimate_spatial_multi_gen_N", "W1.listw")) clusterEvalQ(cl, library("MASS")) cut_B <- split(1:1000, sort(rank(1:B) %% P)) system.time( res_par <- clusterApply(cl, cut_B, fun = my_func) ) stopCluster(cl) # merge the results column_0a <- res_par[[1]]$column_0a column_0b <- res_par[[1]]$column_0b column_1 <- res_par[[1]]$column_1 column_2 <- res_par[[1]]$column_2 column_3 <- res_par[[1]]$column_3 column_4 <- res_par[[1]]$column_4 column_5 <- res_par[[1]]$column_5 column_6 <- res_par[[1]]$column_6 column_7 <- res_par[[1]]$column_7 for (k in 2:P) { column_0a <- cbind(column_0a, res_par[[k]]$column_0a) column_0b <- cbind(column_0b, res_par[[k]]$column_0b) column_1 <- cbind(column_1, res_par[[k]]$column_1) column_2 <- cbind(column_2, res_par[[k]]$column_2) column_3 <- cbind(column_3, res_par[[k]]$column_3) column_4 <- cbind(column_4, res_par[[k]]$column_4) column_5 <- cbind(column_5, res_par[[k]]$column_5) column_6 <- cbind(column_6, res_par[[k]]$column_6) column_7 <- cbind(column_7, res_par[[k]]$column_7) } # compute the RRMMSE # The true parameters in the 4 models true_prms_0a <- c(as.vector(beta_true), c(0, 0, 0, 0), as.vector(sigma_bar_d)) true_prms_0b <- c(as.vector(beta_true), c(0, 0, 0, 0), as.vector(sigma_d)) true_prms_1 <- c(as.vector(beta_true), as.vector(RHO_bar_d), as.vector(sigma_bar_d)) true_prms_2 <- c(as.vector(beta_true), as.vector(RHO_bar_d), as.vector(sigma_d)) true_prms_3 <- c(as.vector(beta_true), as.vector(RHO_d), as.vector(sigma_bar_d)) true_prms_4 <- c(as.vector(beta_true), as.vector(RHO_d), as.vector(sigma_d)) # initialization # biais biais_column_0a <- 0 biais_column_0b <- 0 biais_column_1 <- 0 biais_column_2 <- 0 biais_column_3 <- 0 biais_column_4 <- 0 biais_column_5 <- 0 biais_column_6 <- 0 biais_column_7 <- 0 # rmse rmse_column_0a <- 0 rmse_column_0b <- 0 rmse_column_1 <- 0 rmse_column_2 <- 0 rmse_column_3 <- 0 rmse_column_4 <- 0 rmse_column_5 <- 0 rmse_column_6 <- 0 rmse_column_7 <- 0 # beta_matrix <- matrix(0, B, 8) # RHO_matrix <- matrix(0, B, 4) # sigma_matrix <- matrix(0, B, 4) for (k in 1:B) { # column 0 parms_col_0a <- column_0a[, k] biais_column_0a <- biais_column_0a + (parms_col_0a - true_prms_0a) rmse_column_0a <- rmse_column_0a + (parms_col_0a - true_prms_0a) ^ 2 parms_col_0b <- column_0b[, k] biais_column_0b <- biais_column_0b + (parms_col_0b - true_prms_0b) rmse_column_0b <- rmse_column_0b + (parms_col_0b - true_prms_0b) ^ 2 # column 1 # beta_matrix[k, ] <- as.vector(column_1[[k]]$est_beta) # RHO_matrix[k, ] <- as.vector(column_1[[k]]$est_RHO) # sigma_matrix[k, ] <- as.vector(column_1[[k]]$est_SIGMA) parms_col_1 <- column_1[, k] biais_column_1 <- biais_column_1 + (parms_col_1 - true_prms_1) rmse_column_1 <- rmse_column_1 + (parms_col_1 - true_prms_1) ^ 2 # column 2 parms_col_2 <- column_2[, k] biais_column_2 <- biais_column_2 + (parms_col_2 - true_prms_2) rmse_column_2 <- rmse_column_2 + (parms_col_2 - true_prms_2) ^ 2 # column 3 parms_col_3 <- column_3[, k] biais_column_3 <- biais_column_3 + (parms_col_3 - true_prms_3) rmse_column_3 <- rmse_column_3 + (parms_col_3 - true_prms_3) ^ 2 # column 4 parms_col_4 <- column_4[, k] biais_column_4 <- biais_column_4 + (parms_col_4 - true_prms_3) rmse_column_4 <- rmse_column_4 + (parms_col_4 - true_prms_3) ^ 2 # column 5 parms_col_5 <- column_5[, k] biais_column_5 <- biais_column_5 + (parms_col_5 - true_prms_4) rmse_column_5 <- rmse_column_5 + (parms_col_5 - true_prms_4) ^ 2 # column 6 parms_col_6 <- column_6[, k] biais_column_6 <- biais_column_6 + (parms_col_6 - true_prms_4) rmse_column_6 <- rmse_column_6 + (parms_col_6 - true_prms_4) ^ 2 # column 7 parms_col_7 <- column_7[, k] biais_column_7 <- biais_column_7 + (parms_col_7 - true_prms_4) rmse_column_7 <- rmse_column_7 + (parms_col_7 - true_prms_4) ^ 2 } # Bias (res_biais <- data.frame(col_0a = biais_column_0a / B, col_0b = biais_column_0b / B, col_1 = biais_column_1 / B, col_2 = biais_column_2 / B, col_3 = biais_column_3 / B, col_4 = biais_column_4 / B, col_5 = biais_column_5 / B, col_6 = biais_column_6 / B, col_7 = biais_column_7 / B )) # RRMSE (rrmse_column_0a <- sqrt(rmse_column_0a / B) / abs(true_prms_0a) * 100) (rrmse_column_0b <- sqrt(rmse_column_0b / B) / abs(true_prms_0b) * 100) (rrmse_column_1 <- sqrt(rmse_column_1 / B) / abs(true_prms_1) * 100) (rrmse_column_2 <- sqrt(rmse_column_2 / B) / abs(true_prms_2) * 100) (rrmse_column_3 <- sqrt(rmse_column_3 / B) / abs(true_prms_3) * 100) (rrmse_column_4 <- sqrt(rmse_column_4 / B) / abs(true_prms_3) * 100) (rrmse_column_5 <- sqrt(rmse_column_5 / B) / abs(true_prms_4) * 100) (rrmse_column_6 <- sqrt(rmse_column_6 / B) / abs(true_prms_4) * 100) (rrmse_column_7 <- sqrt(rmse_column_7 / B) / abs(true_prms_4) * 100) (res_table1 <- data.frame(col_0a = rrmse_column_0a, col_0b = rrmse_column_0b, col_1 = rrmse_column_1, col_2 = rrmse_column_2, col_3 = rrmse_column_3, col_4 = rrmse_column_4, col_5 = rrmse_column_5, col_6 = rrmse_column_6, col_7 = rrmse_column_7 )) save(res_biais, res_table1, file = "res_table1_big_variance.RData")