Agregar modelo normal para compración

tereom · Oct 2, 2018 · 072e745 · 072e745
1 parent de6e7e6
commit 072e745
Showing 1 changed file with 55 additions and 0 deletions.
diff --git a/R/mrp_party_estimation.R b/R/mrp_party_estimation.R
@@ -419,3 +419,58 @@ model_bern_t_nb <- function(data_jags, n_chains, n_iter, n_burnin, seed_jags){
     n_votes <- apply(fit_jags$BUGSoutput$sims.list$x, 1, sum)
     return(list(fit = fit_jags, n_votes = n_votes))
 }
+
+model_norm <- function(data_jags, n_chains, n_iter, n_burnin, seed_jags){
+    model_string <-
+        "
+    model{
+    for(k in 1:N){
+    x[k] ~ dnorm(n[k] * theta[k], tau / n[k]) T(0, 750)
+    theta[k] <- ilogit(beta_0 + beta_rural * rural[k] +
+    beta_rural_tamano_md * rural[k] * tamano_md[k] +
+    beta_estrato_raw[estrato[k]] + beta_tamano_md * tamano_md[k] +
+    beta_tamano_gd * tamano_gd[k] + beta_tipo_ex * tipo_ex[k] +
+    beta_region * region[k])
+    }
+    beta_0 ~ dnorm(0, 0.5)
+    beta_rural ~ dnorm(0, 0.5)
+    beta_region ~ dnorm(0, 0.5)
+    beta_tamano_md  ~ dnorm(0, 0.5)
+    beta_tamano_gd  ~ dnorm(0, 0.5)
+    beta_tipo_ex  ~ dnorm(0, 0.5)
+    beta_rural_tamano_md  ~ dnorm(0, 0.5)
+    sigma ~ dunif(0, 1)
+    tau <- pow(sigma, -2)
+    beta_0_adj <- beta_0 + mean(beta_estrato_raw[])
+    for(j in 1:n_strata){
+    beta_estrato[j] <-  beta_estrato_raw[j] - mean(beta_estrato_raw[])
+    beta_estrato_raw[j] ~ dnorm(mu_estrato, tau_estrato)
+    }
+    mu_estrato ~ dnorm(0, 0.5)
+    sigma_estrato ~ dunif(0, 1)
+    tau_estrato <- pow(sigma_estrato, -2)
+    }
+    "
+
+    temp_file <- tempfile(pattern = "model_string", fileext = ".txt")
+    cat(model_string, file = temp_file)
+
+    data_jags[["n_regiones"]] <- NULL
+    data_jags[["region_mod"]] <- NULL
+    fit_jags <- R2jags::jags(
+        # inits = jags_inits,
+        data = data_jags,
+        parameters.to.save = c("x", "beta_0", "beta_0_adj", "beta_rural",
+            "beta_tamano_md", "beta_tamano_gd", "beta_tipo_ex",
+            "beta_estrato", "beta_estrato_raw",
+            "sigma", "sigma_estrato", "beta_rural_tamano_md",
+            "beta_region"),
+        model.file = temp_file,
+        n.chains = n_chains,
+        n.iter = n_iter,
+        n.burnin = n_burnin,
+        jags.seed = seed_jags
+    )
+    n_votes <- apply(fit_jags$BUGSoutput$sims.list$x, 1, sum)
+    return(list(fit = fit_jags, n_votes = n_votes))
+}