updated 17.1-17.3 models with R file

Author: PeterLi2016

ARM/Ch.17/17.1_Varying-intercepts, varying-slopes models.R

@@ -0,0 +1,61 @@
+library(rstan)
+
+srrs2 <- read.table ("srrs2.dat", header=T, sep=",")
+mn <- srrs2$state=="MN"
+radon <- srrs2$activity[mn]
+log.radon <- log (ifelse (radon==0, .1, radon))
+floor <- srrs2$floor[mn]       # 0 for basement, 1 for first floor
+n <- length(radon)
+y <- log.radon
+x <- floor
+
+# get county index variable
+county.name <- as.vector(srrs2$county[mn])
+uniq <- unique(county.name)
+J <- length(uniq)
+county <- rep (NA, J)
+for (i in 1:J){
+  county[county.name==uniq[i]] <- i
+}
+
+ # no predictors
+ybarbar = mean(y)
+
+sample.size <- as.vector (table (county))
+sample.size.jittered <- sample.size*exp (runif (J, -.1, .1))
+cty.mns = tapply(y,county,mean)
+cty.vars = tapply(y,county,var)
+cty.sds = mean(sqrt(cty.vars[!is.na(cty.vars)]))/sqrt(sample.size)
+cty.sds.sep = sqrt(tapply(y,county,var)/sample.size)
+
+ # radon varying intercept and slope model
+dataList.1 <- list(N=length(y), y=y, county=county, J=J, x=x)
+radon_vary_inter_slope.sf1 <- stan(file='17.1_radon_vary_inter_slope.stan', data=dataList.1,
+                            iter=1000, chains=4)
+print(radon_vary_inter_slope.sf1)
+
+ # radon correlation model
+radon_correlation.sf1 <- stan(file='17.1_radon_correlation.stan', data=dataList.1,
+                            iter=1000, chains=4)
+print(radon_correlation.sf1)
+
+ # radon multiple varying coefficients model
+X <- cbind(1,x)
+W <- diag(2)
+dataList.1 <- list(N=length(y), y=y, county=county, J=J, X=X, K=2, W=W)
+radon_multi_varying_coef.sf1 <- stan(file='17.1_radon_multi_varying_coef.stan', data=dataList.1,
+                            iter=1000, chains=4)
+print(radon_multi_varying_coef.sf1)
+
+W <- diag (2)
+dataList.2 <- list(N=length(y), y=y, county=county, J=J, x=x, W=W)
+
+ # radon Scaled inverse-Wishart model
+radon_wishart.sf1 <- stan(file='17.1_radon_wishart.stan', data=dataList.2,
+                            iter=1000, chains=4)
+print(radon_wishart.sf1)
+
+ # radon Scaled inverse-Wishart model 2
+radon_wishart2.sf1 <- stan(file='17.1_radon_wishart2.stan', data=dataList.2,
+                            iter=1000, chains=4)
+print(radon_wishart2.sf1)

ARM/Ch.17/17.1_radon_correlation.stan

@@ -12,14 +12,16 @@ parameters {
   real mu_a;
   real mu_b;
   real<lower=-1,upper=1> rho;
-  matrix[J,2] B;
+  vector[2] B_temp;
 }
 model {
   vector[N] y_hat;
   vector[J] a;
   vector[J] b;
-  matrix[J,2] B_hat;
+  matrix[2,J] B_hat;
   matrix[2,2] Sigma_b;
+  vector[2] B_hat_temp;
+  matrix[2,J] B;
 
   mu_a ~ normal(0, 100);
   mu_b ~ normal(0, 100);
@@ -31,18 +33,22 @@ model {
   Sigma_b[2,1] <- Sigma_b[1,2];
 
   for (j in 1:J) {
-    B_hat[j,1] <- mu_a;
-    B_hat[j,2] <- mu_b;
-    B[j,1:2] ~ multi_normal(B_hat[j,],Sigma_b);
+    B_hat[1,j] <- mu_a;
+    B_hat[2,j] <- mu_b;
+    B_hat_temp[1] <- mu_a;
+    B_hat_temp[2] <- mu_b;
+    B_temp ~ multi_normal(B_hat_temp, Sigma_b);
+    B[1,j] <- B_temp[1];
+    B[2,j] <- B_temp[2];
   }
 
   for (j in 1:J) {
-    a[j] <- B[j,1];
-    b[j] <- B[j,2];
+    a[j] <- B[1,j];
+    b[j] <- B[2,j];
   }
 
   for (i in 1:N)
-    y_hat[i] <- a[county[i]] + b[county[i]] * x[i]
+    y_hat[i] <- a[county[i]] + b[county[i]] * x[i];
 
   y ~ normal(y_hat, sigma);
 }

ARM/Ch.17/17.1_radon_multi_varying_coef.stan

@@ -1,49 +1,53 @@
 data {
   int<lower=0> N;
   int<lower=0> J;
+  int<lower=0> K;
   vector[N] y;
-  int<lower=0,upper=1> x[N];
+  matrix[N,K] X;
   int county[N];
+  matrix[K,K] W;
 }
 parameters {
   real<lower=0> sigma;
-  vector<lower=0>[K] sigma_b;
   vector[K] mu_raw;
   vector[K] xi;
-  matrix[J,K] B_raw;
-  matrix[K,K] W;
-  matrix[K,K] Tau_b_raw;
+  corr_matrix[K] Tau_b_raw;
+  vector[K] B_raw_temp;
 }
 model {
   vector[N] y_hat;
   matrix[J,K] B;
   matrix[K,K] rho_b;
   matrix[K,K] Sigma_b_raw;
-  vector[K] mu_raw;
+  vector[K] mu;
+  matrix[J,K] B_raw;
+  vector[K] sigma_b;
 
   mu_raw ~ normal(0, 100);
   xi ~ uniform(0, 100);
 
   mu <- xi .* mu_raw;
 
-  Tau_b_raw ~ wishart(W, K+1);
+  Tau_b_raw ~ wishart(K+1, W);
   Sigma_b_raw <- inverse(Tau_b_raw);
 
   for (j in 1:J) {
-    for (k in 1:K)
+    B_raw_temp ~ multi_normal(mu_raw, Sigma_b_raw);    
+    for (k in 1:K) {
+      B_raw[j,k] <- B_raw_temp[k];
       B[j,k] <- xi[k] * B_raw[j,k];
-    B_raw[j,] ~ multi_normal(mu_raw, Sigma_b_raw);    
+    }
   }
 
   for (k in 1:K) {
     for (k_prime in 1:K)
       rho_b[k,k_prime] <- Sigma_b_raw[k,k_prime] 
         / sqrt(Sigma_b_raw[k,k] * Sigma_b_raw[k_prime,k_prime]);
-    sigma_b[k] <- abs(xi[k]) * sqrt(Sigma_b_raw[k,k]);
+    sigma_b[k] <- fabs(xi[k]) * sqrt(Sigma_b_raw[k,k]);
   }
 
   for (i in 1:N)
-    y_hat[i] <- dot_product(B[county[i],],X[i,]);
+    y_hat[i] <- dot_product(row(B,county[i]),row(X,i));
 
   y ~ normal(y_hat, sigma);
 }

ARM/Ch.17/17.1_radon_vary_inter_slope.stan

@@ -18,7 +18,7 @@ transformed parameters {
   vector[N] y_hat;
 
   for (i in 1:N)
-    y_hat[i] <- a[county[i]] + b[county[i]] * x[i]
+    y_hat[i] <- a[county[i]] + b[county[i]] * x[i];
 }
 model {
   mu_a ~ normal(0, 100);

ARM/Ch.17/17.1_radon_wishart.stan

@@ -4,56 +4,57 @@ data {
   vector[N] y;
   int<lower=0,upper=1> x[N];
   int county[N];
+  matrix[2,2] W;
 }
 parameters {
   real<lower=0> sigma;
-  real<lower=0> sigma_a;
-  real<lower=0> sigma_b;
   real mu_a_raw;
   real mu_b_raw;
-  real xi_a;
-  real xi_b;
-  real<lower=-1,upper=1> rho;
-  matrix[J,2] B_raw;
-  matrix[2,2] W;
-  matrix[2,2] Tau_b_raw;
+  real<lower=0> xi_a;
+  real<lower=0> xi_b;
+  vector[2] B_raw_temp;
+  corr_matrix[2] Tau_b_raw;
 }
 model {
+  real sigma_a;
+  real sigma_b;
   vector[N] y_hat;
   vector[J] a;
   vector[J] b;
   matrix[J,2] B_raw_hat;
-  matrix[J,2] B_raw;
   matrix[2,2] Sigma_b_raw;
+  matrix[J,2] B_raw;
   real mu_a;
   real mu_b;
+  real rho;
 
   mu_a_raw ~ normal(0, 100);
   mu_b_raw ~ normal(0, 100);
   xi_a ~ uniform(0, 100);
   xi_b ~ uniform(0, 100);
-  rho ~ uniform(-1, 1);
 
   mu_a <- xi_a * mu_a_raw;
   mu_b <- xi_b * mu_b_raw;
 
-  Tau_b_raw ~ wishart(W, 3);
+  Tau_b_raw ~ wishart(3, W);
   Sigma_b_raw <- inverse(Tau_b_raw);
 
   sigma_a <- xi_a * sqrt(Sigma_b_raw[1,1]);
-  sigma_b <- xi_b * sqrt(Sigma_b_raw[2,2]):
+  sigma_b <- xi_b * sqrt(Sigma_b_raw[2,2]);
   rho <- Sigma_b_raw[1,2] / sqrt(Sigma_b_raw[1,1] * Sigma_b_raw[2,2]);
 
   for (j in 1:J) {
     B_raw_hat[j,1] <- mu_a_raw;
     B_raw_hat[j,2] <- mu_b_raw;
-    B_raw[j,1:2] ~ multi_normal(B_raw_hat[j,], Sigma_b_raw);
+    B_raw_temp ~ multi_normal(transpose(row(B_raw_hat,j)), Sigma_b_raw);
+    B_raw[j,1] <- B_raw_temp[1];
+    B_raw[j,2] <- B_raw_temp[2];
     a[j] <- xi_a * B_raw[j,1];
     b[j] <- xi_b * B_raw[j,2];
   }
 
   for (i in 1:N)
-    y_hat[i] <- a[county[i]] + b[county[i]] * x[i]
+    y_hat[i] <- a[county[i]] + b[county[i]] * x[i];
 
   y ~ normal(y_hat, sigma);
 }

ARM/Ch.17/17.1_radon_wishart2.stan

@@ -4,21 +4,17 @@ data {
   vector[N] y;
   int<lower=0,upper=1> x[N];
   int county[N];
+  matrix[2,2] W;
 }
 parameters {
   real<lower=0> sigma;
-  real<lower=0> sigma_a;
-  real<lower=0> sigma_b;
   real mu_a_raw;
   real mu_b_raw;
-  real xi_a;
-  real xi_b;
-  real<lower=-1,upper=1> rho;
-  matrix[J,2] B_raw;
-  matrix[2,2] W;
-  matrix[2,2] Tau_b_raw;
-  vector[K] b_0;
-  matrix[N,K] X_0;
+  real<lower=0> xi_a;
+  real<lower=0> xi_b;
+  vector[2] B_raw_temp;
+  corr_matrix[2] Tau_b_raw;
+  real b_0;
 }
 model {
   vector[N] y_hat;
@@ -29,34 +25,38 @@ model {
   matrix[2,2] Sigma_b_raw;
   real mu_a;
   real mu_b;
+  real sigma_a;
+  real sigma_b;
+  real rho;
 
   mu_a_raw ~ normal(0, 100);
   mu_b_raw ~ normal(0, 100);
   xi_a ~ uniform(0, 100);
   xi_b ~ uniform(0, 100);
-  rho ~ uniform(-1, 1);
   b_0 ~ normal(0, 100);
 
   mu_a <- xi_a * mu_a_raw;
   mu_b <- xi_b * mu_b_raw;
 
-  Tau_b_raw ~ wishart(W, 3);
+  Tau_b_raw ~ wishart(3, W);
   Sigma_b_raw <- inverse(Tau_b_raw);
 
   sigma_a <- xi_a * sqrt(Sigma_b_raw[1,1]);
-  sigma_b <- xi_b * sqrt(Sigma_b_raw[2,2]):
+  sigma_b <- xi_b * sqrt(Sigma_b_raw[2,2]);
   rho <- Sigma_b_raw[1,2] / sqrt(Sigma_b_raw[1,1] * Sigma_b_raw[2,2]);
 
   for (j in 1:J) {
     B_raw_hat[j,1] <- mu_a_raw;
     B_raw_hat[j,2] <- mu_b_raw;
-    B_raw[j,1:2] ~ multi_normal(B_raw_hat[j,], Sigma_b_raw);
+    B_raw_temp ~ multi_normal(transpose(row(B_raw_hat,j)), Sigma_b_raw);
+    B_raw[j,1] <- B_raw_temp[1];
+    B_raw[j,2] <- B_raw_temp[2];
     a[j] <- xi_a * B_raw[j,1];
     b[j] <- xi_b * B_raw[j,2];
   }
 
   for (i in 1:N)
-    y_hat[i] <- dot_product(b_0,X_0[i,]) + a[county[i]] + b[county[i]] * x[i]
+    y_hat[i] <- b_0 + a[county[i]] + b[county[i]] * x[i];
 
   y ~ normal(y_hat, sigma);
 }

ARM/Ch.17/17.2_radon_correlation.stan

@@ -4,6 +4,7 @@ data {
   vector[N] y;
   int<lower=0,upper=1> x[N];
   int county[N];
+  vector[N] u;
 }
 parameters {
   real<lower=0> sigma;
@@ -14,14 +15,15 @@ parameters {
   real g_b_0;
   real g_b_1;
   real<lower=-1,upper=1> rho;
-  matrix[J,2] B;
+  vector[2] B_temp;
 }
 model {
   vector[N] y_hat;
   vector[J] a;
   vector[J] b;
   matrix[J,2] B_hat;
   matrix[2,2] Sigma_b;
+  matrix[J,2] B;
 
   g_a_0 ~ normal(0, 100);
   g_a_1 ~ normal(0, 100);
@@ -37,16 +39,15 @@ model {
   for (j in 1:J) {
     B_hat[j,1] <- g_a_0 + g_a_1 * u[j];
     B_hat[j,2] <- g_b_0 + g_b_1 * u[j];
-    B[j,1:2] ~ multi_normal(B_hat[j,],Sigma_b);
-  }
-
-  for (j in 1:J) {
+    B_temp ~ multi_normal(transpose(row(B_hat,j)),Sigma_b);
+    B[j,1] <- B_temp[1];
+    B[j,2] <- B_temp[2];
     a[j] <- B[j,1];
     b[j] <- B[j,2];
   }
 
   for (i in 1:N)
-    y_hat[i] <- a[county[i]] + b[county[i]] * x[i]
+    y_hat[i] <- a[county[i]] + b[county[i]] * x[i];
 
   y ~ normal(y_hat, sigma);
 }