covarience_adventure_attempt.R

# linc_ig = alpha_g + beta_g * educ_i for individual i in cluster g
# want to estimate alpha_g and beta_g for each cluster

male_fit <- map2stan(                     # create output called male_fit
  alist(
    linc ~ dnorm( mu , sigma ),           # likelihood: log(income) has normal distribution
    mu <- alpha_attract[ornt] +           # mean mu is linearly related to alpha_clusterName + beta_clusterName * education
          beta_attract[ornt]*educ,
    c(alpha_attract,beta_attract)[ornt] ~ dmvnorm2( c(alpha,beta) , sigma_attract , Rho ),    # alpha_clusterName and beta_clusterName are jointly normaly distributed with means alpha and beta and standard deviations sigma_clusterName and correlation Rho 
    alpha ~ dnorm(11.235,1),              # prior: alpha is normally distributed with mean log(75738) and standard deviation 1
    beta ~ dnorm(0.2,1),                  # prior: beta is normally distributed with mean 0.2 and standard deviation 1
    sigma ~ dcauchy(0,2),                 # prior: sigma is distrubted half cauchy with centre 0 and scale 2
    sigma_attract ~ dcauchy(0,2),         # prior: sigma_clusterName are half cauchy distributed with centres 0 and scales 2
    Rho ~ dlkjcorr(2)                     # prior: Rho is distributed LKJcorr with scale parameter 2
  ) ,
data=d , warmup=1000 , iter=5000 , chains=1 , cores=1 )   # keep it simple for debugging