In dmbates/RePsychLing: Data sets from Psychology and Linguistics experiments

library(lme4)
library(RePsychLing)
library(knitr)
library(parallel)
library(rstan)
opts_chunk$set(comment=NA)
options(width=92,show.signif.stars = FALSE)

Load data

mm <- model.matrix(~ sze*(spt+obj+grv)*orn, data=KKL)
spt_orn <- mm[ ,11]
obj_orn <- mm[, 12]
grv_orn <- mm[, 13]

KKL$spt_orn<-spt_orn
KKL$obj_orn<-obj_orn
KKL$grv_orn<-grv_orn

dat <- 
  c(unclass(subset(KKL, select = c(lrt,
                                   sze, spt, obj, grv, orn, 
                                   spt_orn, obj_orn, grv_orn,
                                   sze_spt, sze_obj, sze_grv, sze_orn,
                                   sze_spt_orn, sze_obj_orn, sze_grv_orn))),
    subj = list(as.integer(KKL$subj)),
    N = nrow(KKL),
    I = nlevels(KKL$subj))

Maximal model

if (file.exists("../data/KKL0_maxstanresults.rda")) {
  load("../data/KKL0_maxstanresults.rda")
} else {
KKLmaxmodel <- stan("KKLmaxmodel.stan", 
                   data = dat, 
                   chains = 0)

## Adjust cores to match your computer:
sflist <- 
  mclapply(1:4, mc.cores = 4, 
           function(i) stan(fit = KKLmaxmodel, 
                            data = dat,
                            iter=2000,
                            chains = 1, 
                            chain_id = i, 
                            refresh = -1))

KKL0_maxstan <- sflist2stanfit(sflist)

  ## save as matrix
KKL0maxmat<-as.matrix(KKL0_maxstan)
KKL0maxreduced<-KKL0maxmat[,1:89]
param_names<-colnames(KKL0maxreduced)
resultsKKL0max<-matrix(rep(NA,89*3),ncol=3)
for(i in 1:89){
  resultsKKL0max[i,]<-c(mean(resultsKKL0max[,i]),
                     quantile(resultsKKL0max[,i],
                     probs=c(0.025,0.975)))
}

resultsKKL0max<-data.frame(resultsKKL0max)
colnames(resultsKKL0max)<-c("mean","lower","upper")

Save results:

## saved in RePsychLing data dir:
save(resultsKKL0max,
     file="../data/KKL0_maxstanresults.rda",
     compress="xz")
}

Final model

if (file.exists("../data/KKL0_finstanresults.rda")) {
  load("../data/KKL0_finstanresults.rda")
} else {
KKLfinmodel <- stan("KKLfinmodel.stan", 
                   data = dat, 
                   chains = 0)

## Adjust cores to match your computer:
sflist <- 
  mclapply(1:4, mc.cores = 4, 
           function(i) stan(fit = KKLfinmodel, 
                            data = dat,
                            iter=2000,
                            chains = 1, 
                            chain_id = i, 
                            refresh = -1))

KKL0fin_stan <- sflist2stanfit(sflist)

## get only relevant part:
  ## save as matrix
KKL0finmat<-as.matrix(KKL0fin_stan)
params<- c(1:20,546:548)
KKL0finreduced<-KKL0finmat[,params]

param_names<-c("Intercept","sze",
               "spt","obj",
               "grv","orn",
               "sze_spt","sze_obj",
               "sze_grv","sze_orn",
               "spt_orn","obj_orn",
               "grv_orn","sze_spt_orn",
               "sze_obj_orn",
               "sze_grv_orn","sigma_e","sigma_u1",
               "sigma_u2","sigma_u3","sigma_u_obj",
               "sigma_u_spt_orn",
               "sigma_u_orn")


resultsKKL0fin<-matrix(rep(NA,
                           length(params)*3),
                       ncol=3)

for(i in 1:length(params)){
  resultsKKL0fin[i,]<-c(mean(KKL0finreduced[,i]),
                     quantile(KKL0finreduced[,i],
                     probs=c(0.025,0.975)))
}

rownames(resultsKKL0fin)<-param_names

resultsKKL0fin<-data.frame(resultsKKL0fin)
colnames(resultsKKL0fin)<-c("mean","lower",
                            "upper")

## saved in RePsychLing data dir:xxx
save(resultsKKL0fin,
     file="../data/KKL_finstanresults.rda",
     compress="xz")

## extract correlation parameters:
e<-extract(KKL0fin_stan)
n_samp <- dim(e[[1]])[1]
## correlations subject:
cor_listu <- lapply(seq(n_samp), 
                   function(i) {
  L <- e$L_u[i, , ]
  t(L) %*% L
})

# number of random effects
n_ran <- nrow(cor_listu[[1]])

## final model:
kkl4<-lmer(lrt ~ sze * (spt + obj + grv) * orn + (spt + grv | subj) + (0 +  
    obj | subj) + (0 + orn | subj) + (0 + spt_orn | subj),KKL)

KKL4_cor_s <- attr(VarCorr(kkl4)$subj,"correlation")

# create indices of correlation coefficients (upper triangular)
cor_idx_s <- which(upper.tri(KKL4_cor_s), arr.ind = TRUE)

cor_summary<-matrix(rep(NA,3*4),ncol=4)
rnames<-rep(NA,3)
for (i in 1:3) {
  idx <- cor_idx_s[i, , drop = FALSE]
  KKL4_cor <- KKL4_cor_s[idx]
  KKL4s_cor <- sapply(cor_listu, "[", idx)
  ## save means and quantiles, Stan:
cor_summary[i,2:4]<-c(mean(KKL4s_cor),
                   quantile(KKL4s_cor,
                           probs=c(0.025,
                                   0.975))) 
## save lmer corr. estimate:
cor_summary[i,1]<-KKL4_cor
## save rownames:
rnames[i]<-paste(rownames(KKL4_cor_s)[idx], collapse = " : ")
}
rownames(cor_summary)<-rnames
colnames(cor_summary)<-c("lmer est.",
                         "Stan est.",
                         "2.5th%ile",
                         "97.5th%ile")

save(cor_summary,file="../data/KKLfincorsubj.rda",compress="xz")

dmbates/RePsychLing documentation built on May 15, 2019, 9:19 a.m.