SHARE
TWEET

Untitled

a guest Sep 4th, 2019 89 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. library(devtools)
  2. library(rjags)
  3. library(ggplot2)
  4. library(dplyr)
  5.  
  6. ur <- "https://raw.githubusercontent.com/PerceptionCognitionLab/data0/master/contexteffects/FlankerStroopSimon/cleaning.R"
  7. source_url(ur)
  8.  
  9.  
  10.  
  11.  
  12. ICC_lsm <- "model{
  13.  
  14. for(j in 1:J){
  15.  
  16. # latent betas
  17. beta_raw_l[j] ~  dnorm(0, 1)
  18.  
  19. # random effect
  20. beta_l[j] <- fe_mu + tau_mu * beta_raw_l[j]
  21.  
  22. # cholesky
  23. z2[j] ~ dnorm(0, 1)
  24. beta_raw_s[j] = rho12 * beta_raw_l[j] + sqrt(1 - rho12^2) * z2[j]
  25.  
  26. beta_s[j] <- fe_sd + tau_sd * beta_raw_s[j]
  27.  
  28. }
  29.  
  30. for(i in 1:N){
  31.  
  32. # likelihood
  33. y[i] ~ dnorm(beta_l[ID[i]], 1/exp(beta_s[ID[i]])^2)
  34.  
  35. }
  36.  
  37.  
  38. # fixed effects priors
  39. fe_mu ~ dnorm(0, 1)
  40. fe_sd ~ dnorm(0, 1)
  41.  
  42. # random effects priors
  43. tau_mu ~ dt(0, pow(prior_scale,-2), 10)T(0,)
  44. tau_sd ~ dt(0, pow(prior_scale,-2), 10)T(0,)
  45.  
  46.  
  47.  
  48. # prior for RE correlation
  49. fz ~ dnorm(0, 1)
  50. rho12 = tanh(fz)
  51.  
  52. }"
  53.  
  54.  
  55.  
  56. incongruent <- subset(stroop, congruency == "congruent")
  57.  
  58.  
  59. model <- jags.model(textConnection(ICC_lsm),
  60.                     inits = list(fe_mu = mean(incongruent$rt)),
  61.                     data = list(y = as.numeric((incongruent$rt)),
  62.                                 ID = incongruent$ID,
  63.                                 J = length(unique(incongruent$ID)),  
  64.                                 N = length(incongruent$rt),
  65.                                 prior_scale = 1))
  66.  
  67. # fit for motivating example
  68. fit_1 <- coda.samples(model,
  69.                       n.iter = 5000,
  70.                       variable.names=c("fe_mu",
  71.                                        "fe_sd",
  72.                                        "tau_mu",
  73.                                        "tau_sd",
  74.                                        "beta_s",
  75.                                        "beta_l"))
  76.  
  77. # samples
  78. fit_1_samples <- do.call(rbind.data.frame, fit_1)
  79.  
  80. # mean ICC
  81. mu_icc <- mean(fit_1_samples$tau_mu^2 / (fit_1_samples$tau_mu^2 + exp(fit_1_samples$fe_sd)^2))
  82.  
  83. # colnames
  84. col_names <- colnames(fit_1_samples)
  85.  
  86. # person specific SD
  87. beta_s <- exp(fit_1_samples[, grep("beta_s", col_names)])
  88.  
  89. # person specific ICCs
  90. posterior_icc <- apply(beta_s, MARGIN = 2, FUN = function(x){ fit_1_samples$tau_mu^2  /  (fit_1_samples$tau_mu^2  + x^2)}  )
  91.  
  92. # 90 % intervals for ICCs
  93. credible_intervals <- t(apply(posterior_icc, 2, FUN = function(x)  {quantile(x, c(0.05, 0.95)) }))
  94.  
  95. posterior_means <- (colMeans(posterior_icc))
  96.  
  97.  
  98. C <- cbind(posterior_means, credible_intervals) %>%
  99.   as.data.frame() %>%
  100.   arrange(posterior_means) %>%
  101.   mutate(index = as.factor(1:121),
  102.          sig = as.factor(ifelse(.$`5%` < mu_icc & .$`95%` < mu_icc, 1,  
  103.                                 ifelse(.$`5%` > mu_icc & .$`95%` > mu_icc, 2, 3)))) %>%
  104.   ggplot(aes(x = index,y  = posterior_means)) +
  105.   geom_hline(yintercept = mu_icc, alpha = 0.5, linetype = "dotted") +
  106.   geom_errorbar(aes(ymin = `5%`,
  107.                     ymax = `95%`, color = sig)) +
  108.   geom_point(color = "black") +
  109.   scale_color_manual(values = c("#56B4E9", "#009E73", "grey70")) +
  110.   theme_bw(base_family = "Times") +
  111.   theme(panel.grid = element_blank(),
  112.         panel.grid.major.x = element_blank(),
  113.         panel.grid.minor.y = element_blank(),
  114.         axis.text.x = element_blank(),
  115.         axis.title = element_text(size = 14),
  116.         title = element_text(size = 14),
  117.         strip.background = element_rect(fill = "grey94"),
  118.         strip.text = element_text(size = 14), legend.position = "none") +
  119.   xlab("") +
  120.   scale_x_discrete(expand = c(0.015, 0.015)) +
  121.   scale_y_continuous(limits = c(.05, .6),
  122.                      labels = seq(.1, .6, .1) ,
  123.                      breaks = seq(.1, .6, .1)) +
  124.   xlab("Ascending Index") +
  125.   ylab("Intraclass Correlation Coefficient(s)")
  126. C
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top