Main.Rmd

Libraries
```{r}
  library(purrr)
  library(readr)
  library(ggplot2)
  library(langcog)
  library(boot)
#  library(lazyeval)
  library(dplyr)
  library(tidyr)
  library(wordbankr)
  library(directlabels)
#  library(scales)
  library(stringr)
  library(lmtest)
  library(rwebppl)
  library(jsonlite)
  library(nlme)
  library(feather)
  library(broom)
  library(HDInterval)
  library(BBmisc)
```

Writings


Import functions 

```{r}
source(paste(getwd(),"/helper_functions/all_helper.r",sep = ""), chdir = T)
```

Make a data set that combines all predictores

```{r, message=FALSE, warning=FALSE}

languages = c("Croatian","Danish","English (American)", "Italian", "Norwegian", "Russian","Spanish","Turkish")

#This command create combined data (all predcitores) and save them in the out_files folder (under: lang_name_combo.csv)

#create_combo(languages = languages)

#Now to save time, we will load already saved files 
Cro_combo <- read.csv("out_files/Croatian_combo.csv")
Dan_combo <- read.csv("out_files/Danish_combo.csv")
Eng_combo <- read.csv("out_files/English (American)_combo.csv")
Ita_combo <- read.csv("out_files/Italian_combo.csv")
Nor_combo <- read.csv("out_files/Norwegian_combo.csv")
Rus_combo <- read.csv("out_files/Russian_combo.csv")
Spa_combo <- read.csv("out_files/Spanish_combo.csv")
Tur_combo <- read.csv("out_files/Turkish_combo.csv")


language_dict = data.frame(language=c(1,2,3,4,5,6,7,8),
                           language_name=c("Croatian","Danish","English", "Italian", "Norwegian", "Russian","Spanish","Turkish"))
```


Select data for first month 
```{r}

data_original <- list(Cro_combo, Dan_combo, Eng_combo, Ita_combo, Nor_combo, Rus_combo, Spa_combo, Tur_combo)

data_first_month <- data.frame()

for (i in 1:length(data_original)) {
  
   #get first month
  first_month <- (data_original[[i]] %>%
    distinct(age))$age[1]
  
  #select the intersction of all predictors
  data_lang <- data_original[[i]] %>%
    filter(age == first_month) %>%
    mutate(language = i) %>%
    select(-age)
  
  data_first_month <- bind_rows(data_first_month, data_lang)
  
}

data_first_month <- data_first_month %>%
  left_join(language_dict) %>%
  select(-language) %>%
  rename(language = language_name) %>%
  select(-PAT_assoc, -PAT_phono_t1, -PAT_phono_t2, -PAT_phono_nt)

```

Here do the power-law analysis and the log-log plots

```{r}
#Data for plot
degreeDist <- data.frame(matrix(ncol = 4, nrow = 0))
dist_names <- c("x", "y", "dimension", "language")
colnames(degreeDist) <- dist_names

#Parameters and test
degreeTest <- data.frame(matrix(ncol = 5, nrow = 0))
test_names <- c("xMin", "alpha", "pVal", "dimension","language")
colnames(degreeTest) <- test_names


for (i in language_dict$language_name) {
  #Select data for the language
  data_lang <- data_first_month %>%
    filter(language == i) %>%
    select(language, PAC_assoc, PAC_phono_t2) %>%
    rename(Sem = PAC_assoc, Phono = PAC_phono_t2) %>%
    filter (!is.na(Sem),!is.na(Phono)) 
  
  #fit and derive parameters for power law
  semList <- data_lang$Sem[data_lang$Sem != 0]
  sem_pl = displ$new(semList)
  sem_est = estimate_xmin(sem_pl)
  sem_pl$setXmin(sem_est)
  sem_dist = plot(sem_pl) %>%
    mutate(dimension='Sem', language  =  i)
  
  #bootstrap to get p-value
  sem_boot = bootstrap_p(sem_pl, no_of_sims=1000, threads=2)
  
  sem_test <- data.frame(as.numeric(sem_pl$xmin), as.numeric(sem_pl$pars), as.numeric(sem_boot$p), 'Sem', i)
  colnames(sem_test) <- test_names
  
  phonoList <- data_lang$Phono[data_lang$Phono != 0]
  phono_pl = displ$new(phonoList)
  phono_est = estimate_xmin(phono_pl)
  phono_pl$setXmin(phono_est)
  phono_dist = plot(phono_pl) %>%
    mutate(dimension='Phono', language = i)
  
  phono_boot = bootstrap_p(phono_pl, no_of_sims=1000, threads=2)
  
  phono_test <- data.frame(as.numeric(phono_pl$xmin), as.numeric(phono_pl$pars), as.numeric(phono_boot$p), 'Phono', i)
  colnames(phono_test) <- test_names
  
  dist_lang <- bind_rows(sem_dist, phono_dist)
  test_lang <- bind_rows(sem_test, phono_test)
  
  degreeDist <- bind_rows(degreeDist, dist_lang) 
  degreeTest <- bind_rows(degreeTest, test_lang) 
  
}

#feather::write_feather(degreeDist, "saved_data/degreeDist.feather")
degreeDist <- feather::read_feather("saved_data/degreeDist.feather")

feather::write_feather(degreeTest, "saved_data/degreeTest.feather")
degreeTest <- feather::read_feather("saved_data/degreeTest.feather")

#plot cumulative distributions
ggplot(data = degreeDist,  aes(x=x, y=y, col=dimension))+
  facet_grid(. ~ language)+#, scales = "free") +
  geom_point(#aes(colour = lexical_class),
    #colour = solarized_palette(1),
    size = 0.5, alpha = 0.5) +
  scale_y_log10() + scale_x_log10() +
  theme(aspect.ratio = 1)
  #scale_x_continuous(limits=c(-2,5))+
  #geom_smooth(method = "lm", colour = "grey3", se=FALSE)+
  #scale_colour_solarized(name = "") 


```

Here do the statistics 

```{r}

#Total nounds used
total_nouns <- data_first_month %>%
  group_by(language) %>%
  summarise(total=n())

#Nounns that have unilemma
total_unilemas <- data_first_month %>%
  filter(!is.na(uni_lemma)) %>%
  group_by(language) %>%
  summarise(translated=n())

#Nouns with unilemas that have entries in Nelson's association norms 
total_association <- data_first_month %>%
  filter(!is.na(PAC_assoc)) %>%
  group_by(language) %>%
  summarise(normed=n())

total_nouns <- total_nouns %>%
  left_join(total_unilemas) %>%
  left_join(total_association)

feather::write_feather(total_nouns, "saved_data/statistics.feather")
total_nouns <- feather::read_feather("saved_data/statistics.feather")


  
```


Prepare data for age of acquisition analysis (regressions)

```{r}
data_original <- list(Cro_combo, Dan_combo, Eng_combo, Ita_combo, Nor_combo, Rus_combo, Spa_combo, Tur_combo)

data_regression <- data.frame()

for (i in 1:length(data_original)) {
  
   #get first month
  first_month <- (data_original[[i]] %>%
    distinct(age))$age[1]
  
  #select the intersction of all predictors
  data_lang <- data_original[[i]] %>%
    filter(age == first_month) %>% #the first month contains all the data 
    filter (#!is.na(freq), 
            !is.na(length), 
            !is.na(PAC_assoc),
            !is.na(PAC_phono_t1),
            !is.na(PAC_phono_t2)) %>%
    mutate(language = i) %>%
    select(-age)
  # replace NA in PAT with O (NA is when the utility cannot be computed, e.g., a candidate word is not linked to any previousely learned word)
  
  data_regression <- bind_rows(data_regression, data_lang)
  
}

data_regression <- data_regression %>%
  left_join(language_dict) %>%
  select(-language) %>%
  rename(language = language_name) %>%
  select(-PAT_assoc, -PAT_phono_t1, -PAT_phono_t2, -PAT_phono_nt)

```


Derive the statistics from Barginsky et al. 2016, and add to our data

```{r}

uni_joined <- feather::read_feather("saved_data/uni_joined.feather")
summ <- feather::read_feather("saved_data/uni_model_data.feather") 

counts_mika <- uni_joined %>%
  distinct(language, uni_lemma, frequency) %>%
  rename(freq_mika = frequency)

data_regression_all <- data_regression %>%
  left_join(counts_mika)


#Here I should count how much data we missed when using the mika_freq

wordn_byLang <- data_regression_all %>%
  group_by(language) %>%
  dplyr::summarise(total=n())

missing <- data_regression_all %>%
  filter(is.na(freq_mika)) %>%
  group_by(language)  %>%
  dplyr::summarise(missing=n()) %>%
  left_join(wordn_byLang) %>%
  mutate(prop = missing/total) %>%
  mutate(remain = total-missing )

#There are large missing data in English and Russian (about 30%). For English, I will use instead the Ping Li of Penn State frequnecy count (availabale from CHILDES online). The problem remains in Russian

freq_eng <- data_regression_all %>%
  filter(language=='English') %>%
  select(-freq_mika) %>%
  mutate(freq_mika = freq)
  
  
data_regression_all_bis <- data_regression_all %>%
  filter(language != 'English')

data_regression_all_bis <- bind_rows(data_regression_all_bis, freq_eng)




```

Compute AOA  (50% acquisition rate)

```{r}
make_aoa <- function(lang, lang_form = "WS", lex_class = "nouns") {
  item_data <- get_lang_item_data(lang = lang,
                                  lang_form = lang_form,
                                  lex_class = lex_class)
  admin_data <-
    get_lang_admin_data(lang = lang, lang_form = lang_form)
  instr_data <-
    get_lang_instr_data(lang = lang, lang_form = lang_form)
  word_aoa <-
    get_lang_aoa(item_data = item_data,
                 admin_data = admin_data,
                 instr_data = instr_data)
  
  return(word_aoa)
}

aoa = data.frame()

for (lang in languages) {
aoa_lang <- make_aoa(lang) %>%
  select(item, age) %>%
  mutate(language = lang)

aoa <- bind_rows(aoa, aoa_lang)

}

aoa <- aoa %>%
  rename(AoA = age) 

aoa_all <- aoa
aoa_all$language[aoa_all$language == 'English (American)'] <- 'English'


#all data with AoA
data_all <- data_regression_all %>%
  left_join(aoa_all) 

#Scale all the predictors here (Wait I should scale for each language separately)
data_all_scaled <- data_all %>%
  group_by(language) %>%
  mutate_at(c("length", "freq", "PAC_assoc", "PAC_phono_t1", "PAC_phono_t2", "PAC_phono_nt", "freq_mika"), funs(as.numeric(scale(.)))) %>%
  select(item, definition, uni_lemma, learned,   IPA,    length,        freq,   PAC_assoc, PAC_phono_t1, PAC_phono_t2, PAC_phono_nt, freq_mika, AoA, language)

#############
#all data with AoA
data_all_bis <- data_regression_all_bis %>%
  left_join(aoa_all) 

#Scale all the predictors here
data_all_scaled_bis <- data_all_bis %>%
  group_by(language) %>%
  mutate_at(c("length", "freq", "PAC_assoc", "PAC_phono_t1", "PAC_phono_t2", "PAC_phono_nt", "freq_mika"), funs(as.numeric(scale(.)))) %>%
  select(item, definition, uni_lemma, learned,   IPA,    length,        freq,   PAC_assoc, PAC_phono_t1, PAC_phono_t2, PAC_phono_nt, freq_mika, AoA, language)
#############


#Save at this point 

feather::write_feather(data_all_scaled, "saved_data/data_all_scaled.feather")
data_all_scaled <- feather::read_feather("saved_data/data_all_scaled.feather")

#With English frequency
feather::write_feather(data_all_scaled_bis, "saved_data/data_all_scaled_bis.feather")
data_all_scaled_bis <- feather::read_feather("saved_data/data_all_scaled_bis.feather")
```


Plot data for all languages
```{r fig.width=4, fig.height=2}
#tidy data
data_long <- data_all_scaled %>%
  gather(predictor, value, length:freq_mika)

data_long %>%
  #filter(predictor %in% c("PAC_assoc", "PAC_phono_t2", "freq_mika", "length")) %>%
  filter(predictor %in% c("PAC_assoc", "PAC_phono_t2")) %>%
ggplot(aes(x=value, y=AoA))+
  facet_grid(predictor ~ language)+#, scales = "free") +
  geom_jitter(#aes(colour = lexical_class),
    colour = solarized_palette(1),
    size = 0.5, alpha = 0.5) +
  #scale_x_continuous(limits=c(-2,5))+
  geom_smooth(method = "lm", colour = "grey3", se=FALSE)+
  scale_colour_solarized(name = "") 
```


Regression Analysis 
```{r}

#apply lm for each langauge 
models <- data_all_scaled_bis %>%
  group_by(language) %>%
  do (fit_aoa = lm(AoA ~  freq_mika + length + PAC_assoc+ PAC_phono_t2, data=.) )
stat_estim <- tidy(models, fit_aoa)

models_confint <- data_all_scaled_bis %>%
  group_by(language) %>%
  do (fit_ci = confint(lm(AoA ~  freq_mika + length + PAC_assoc+ PAC_phono_t2, data=.) )) 
stat_conf <- tidy(models_confint, fit_ci) %>%
  rename(term = .rownames)

stat_all <- stat_estim %>%
  left_join(stat_conf) %>%
  filter(term != '(Intercept)') %>%
  rename(predictor = term)

stat_all$predictor <- mapValues(stat_all$predictor, 
                                 from = c("freq_mika", "length","PAC_assoc","PAC_phono_t2"), 
                                 to = c("frequency", "length", "semPAC", "phonoPAC"))

stat_all$predictor <- factor(stat_all$predictor, levels = c("frequency", "length", "semPAC", "phonoPAC"))

feather::write_feather(stat_all, "saved_data/regressions_bis.feather")
stat_all <- feather::read_feather("saved_data/regressions_bis.feather")


#General regression with language as a random variable (with length and frequency as covariate)
general <- lme(AoA ~  freq_mika + length + PAC_assoc+ PAC_phono_t2 , random=~1|language,  data=data_all_scaled_bis, na.action=na.omit)
general_ci <- data.frame(intervals(general)$fixed) 
general_ci$term <- rownames(general_ci)
general_ci <- general_ci %>%
  filter(term != '(Intercept)') %>%
  rename(predictor = term,
         estimate = est.)


general_ci$predictor <- mapValues(general_ci$predictor, 
                                 from = c("freq_mika", "length","PAC_assoc","PAC_phono_t2"), 
                                 to = c("frequency", "length", "semPAC", "phonoPAC"))

general_ci$predictor <- factor(general_ci$predictor, levels = c("frequency", "length", "semPAC", "phonoPAC"))

feather::write_feather(general_ci, "saved_data/regression_all_bis.feather")
general_ci <- feather::read_feather("saved_data/regression_all_bis.feather")


#General regression with language as a random variable (without the length and frequency as covariate)
general_noCov <- lme(AoA ~ PAC_assoc+ PAC_phono_t2 , random=~1|language,  data=data_all_scaled_bis, na.action=na.omit)
general_ci_noCov <- data.frame(intervals(general_noCov)$fixed) 
general_ci_noCov$term <- rownames(general_ci_noCov)
general_ci_noCov <- general_ci_noCov %>%
  filter(term != '(Intercept)') %>%
  rename(predictor = term,
         estimate = est.)


general_ci_noCov$predictor <- mapValues(general_ci_noCov$predictor, 
                                 from = c("PAC_assoc","PAC_phono_t2"), 
                                 to = c("semPAC", "phonoPAC"))

general_ci_noCov$predictor <- factor(general_ci_noCov$predictor, levels = c("frequency", "length", "semPAC", "phonoPAC"))

feather::write_feather(general_ci_noCov, "saved_data/regression_all_bis_noCov.feather")
general_ci_noCov <- feather::read_feather("saved_data/regression_all_bis_noCov.feather")

```


Plot coefficients

```{r fig.width=4, fig.height=3}
ggplot(stat_all, aes(x = predictor, y = estimate)) +
  geom_pointrange(aes(ymin = X2.5.., ymax = X97.5.., y = estimate, col = predictor))+
  geom_hline(yintercept = 0, color = "grey", linetype = "dashed")+
  facet_wrap(~language, ncol=4)  +
  coord_flip() +
  guides(colour=FALSE)+
  scale_colour_solarized() +
  theme_bw()+
  theme(aspect.ratio = 0.7)


ggplot(general_ci, aes(x = predictor, y = estimate)) +
  geom_pointrange(aes(ymin = lower, ymax = upper, col = predictor))+
  geom_pointrange(data= general_ci_noCov, aes(ymin = lower, ymax = upper, col = predictor), linetype=2)+
  geom_hline(yintercept = 0, color = "grey", linetype = "dashed")+
  coord_flip() +
  guides(colour=FALSE)+
  scale_colour_solarized() +
  theme_bw()+
  theme(aspect.ratio = 0.7)

ggplot(general_ci_noCov, aes(x = predictor, y = estimate)) +
  geom_pointrange(aes(ymin = lower, ymax = upper, col = predictor))+
  geom_hline(yintercept = 0, color = "grey", linetype = "dashed")+
  coord_flip() +
  guides(colour=FALSE)+
  scale_colour_solarized() +
  theme_bw()+
  theme(aspect.ratio = 0.7)


```
Pre-process data for network growth models

```{r}

data_original <- list(Cro_combo, Dan_combo, Eng_combo, Ita_combo, Nor_combo, Rus_combo, Spa_combo, Tur_combo)

data <- vector("list", length = length(data_original))

for (i in 1:length(data_original)) {
  
  #get first month
  first_month <- (data_original[[i]] %>%
    distinct(age))$age[1]
  
  #select the intersction of all predictors
  data_lang <- data_original[[i]] %>%
    filter (#!is.na(freq), 
            !is.na(length), 
            !is.na(PAC_assoc),
            !is.na(PAC_phono_t1),
            !is.na(PAC_phono_t2)) %>%
   #Restrict data to up to 26 months in order to avoid overfitting to final months (which have small datapoints). In this we followed Hills et al. 2009.
  #Remove first month (16) in order to compare PAT to PAC (PAT utility can only be computed based on words learned in a previous months)
    filter (age < 27, 
            age > first_month) %>%
    mutate(language = i) %>%
    left_join(language_dict) %>%
    select(-language) %>%
    rename(language = language_name) %>%
    left_join(counts_mika) %>% #(but I am not using frequency at this stage, so this is not necessary)
    select(-language) #%>%
    #filter(!is.na(freq_mika))
  

  # replace NA in PAT with O (NA is when the utility cannot be computed, e.g., a candidate word is not linked to any previousely learned word)
  data_lang$PAT_assoc[is.na(data_lang$PAT_assoc)] <- 0
  data_lang$PAT_phono_t1[is.na(data_lang$PAT_phono_t1)] <- 0
  data_lang$PAT_phono_t2[is.na(data_lang$PAT_phono_t2)] <- 0
  
  data[[i]] <- data_lang
  
}

```

Compute z-score and test network growth 

```{r}

zscores <- data.frame(matrix(ncol = 7, nrow = 0))
x <- c("mean", "parameter", "p.value", "lower","upper", "model", "language")
colnames(zscores) <- x

for (j in 1:length(data)) {

  age_list <- (data[[j]] %>% 
                 distinct(age))$age

  data_z <- data.frame()
  
  for (i in age_list) {
    data_age <- data[[j]] %>%
      filter(age==i)  %>%
      mutate(PAC_assoc_z=scale(PAC_assoc),
             PAT_assoc_z=scale(PAT_assoc),
             
             PAC_phono_t1_z=scale(PAC_phono_t1),
             PAT_phono_t1_z=scale(PAT_phono_t1),
             
             PAC_phono_t2_z=scale(PAC_phono_t2),
             PAT_phono_t2_z=scale(PAT_phono_t2)) %>%
      
      select(age, item, learned, PAC_assoc_z, PAT_assoc_z, PAC_phono_t1_z, PAT_phono_t1_z, PAC_phono_t2_z, PAT_phono_t2_z) 
    
    data_z <- bind_rows(data_z, data_age) 
  }
  
  #Replace NA z-zcore with 0 (this happens when variance is 0, i.e., all data is a constant)
  data_z$PAT_assoc_z[is.na(data_z$PAT_assoc_z)] <- 0
  data_z$PAT_phono_t1_z[is.na(data_z$PAT_phono_t1_z)] <- 0
  data_z$PAT_phono_t2_z[is.na(data_z$PAT_phono_t2_z)] <- 0
  
  data_z <- data_z %>%
    filter(learned == 1)

  #Make table for z-score statistsics
  
  test <- t.test(data_z$PAC_assoc_z, mu=0)
  x1 <- data.frame(as.numeric(test$estimate), as.numeric(test$parameter), as.numeric(test$p.value), as.numeric(test$conf.int[1]), as.numeric(test$conf.int[2]),  'PAC_assoc', j)
  colnames(x1) <- x
  
  test <- t.test(data_z$PAT_assoc_z, mu=0)
  x2 <- data.frame(as.numeric(test$estimate), as.numeric(test$parameter), as.numeric(test$p.value), as.numeric(test$conf.int[1]), as.numeric(test$conf.int[2]), 'PAT_assoc', j)
  colnames(x2) <- x
  
  test <- t.test(data_z$PAC_phono_t1_z, mu=0)
  x3 <- data.frame(as.numeric(test$estimate), as.numeric(test$parameter), as.numeric(test$p.value), as.numeric(test$conf.int[1]), as.numeric(test$conf.int[2]), 'PAC_phono_t1', j)
  colnames(x3) <- x
  
  test <- t.test(data_z$PAT_phono_t1_z, mu=0)
  x4 <- data.frame(as.numeric(test$estimate), as.numeric(test$parameter), as.numeric(test$p.value), as.numeric(test$conf.int[1]),as.numeric(test$conf.int[2]), 'PAT_phono_t1', j)
  colnames(x4) <- x
  
  test <- t.test(data_z$PAC_phono_t2_z, mu=0)
  x5 <- data.frame(as.numeric(test$estimate), as.numeric(test$parameter), as.numeric(test$p.value), as.numeric(test$conf.int[1]), as.numeric(test$conf.int[2]), 'PAC_phono_t2', j)
  colnames(x5) <- x
  
  test<- t.test(data_z$PAT_phono_t2_z, mu=0)
  x6 <- data.frame(as.numeric(test$estimate), as.numeric(test$parameter), as.numeric(test$p.value), as.numeric(test$conf.int[1]), as.numeric(test$conf.int[2]), 'PAT_phono_t2', j)
  colnames(x6) <- x
  
  
 zscores <- bind_rows(zscores, x1, x2, x3, x4, x5, x6)
 
}

zscores <- zscores %>%
  left_join(language_dict) %>% 
  select(-language) %>% 
  rename(language = language_name) %>%
  #Since only phono t2 works, I filter out phono t1
  filter(model != 'PAC_phono_t1',
         model != 'PAT_phono_t1')


```

Plot z-score data 

```{r}
ggplot(zscores, aes(x = model, y = mean)) +
  geom_pointrange(aes(ymin = lower, ymax = upper, y = mean, col = model))+
  geom_hline(yintercept = 0, color = "grey", linetype = "dashed")+
  coord_flip() +
  facet_wrap(~language, ncol=4)  +
  theme(aspect.ratio = 0.7)

```

This section is about combining predictors in a network growth model.

Import models

```{r}
source(paste(getwd(),"/models/all_models.r",sep = ""), chdir = T)
```


Scale data before input to the model
```{r}

data_scaled <- vector("list", length = length(data))

for (i in 1:length(data)) {

   data_lang <- data[[i]] %>%
     mutate_at(vars(length:freq_mika), funs(as.numeric(scale(.)))) 
   
   data_scaled[[i]] <- data_lang
} 

```


Combine all data in one data frame (which we will input to the model)
```{r}
data_combined <- data.frame()
#x <- c("mean", "parameter", "p.value", "lower","upper", "model", "language")
#colnames(zscores) <- x

for (i in 1:length(data_scaled)) {

   data_lang <- data_scaled[[i]] %>%
     mutate (language = i)
   
   data_combined <- bind_rows(data_combined, data_lang)
} 

```

Make models


```{r}
#Isolateed predictors
model_semPAT <-  paste(helper, sem_PAT, optimize, sep = '\n')
model_semPAC <-  paste(helper, sem_PAC, optimize, sep = '\n')
model_phonoPAT <-  paste(helper, phono_PAT, optimize, sep = '\n')
model_phonoPAC <-  paste(helper, phono_PAC, optimize, sep = '\n')

#Combined predictors
model_semPAC_semPAT <-  paste(helper, semPAC_semPAT, optimize, sep = '\n')
model_phonoPAC_phonoPAT <-  paste(helper, phonoPAC_phonoPAT, optimize, sep = '\n')
model_semPAC_phonoPAC <-  paste(helper, semPAC_phonoPAC, optimize, sep = '\n')
model_semPAT_phonoPAT <-  paste(helper, semPAT_phonoPAT, optimize, sep = '\n')
model_semPAT_phonoPAC <-  paste(helper, semPAT_phonoPAC, optimize, sep = '\n')
model_semPAC_phonoPAT <-  paste(helper, semPAC_phonoPAT, optimize, sep = '\n')

model_semPAC_phonoPAC_semPAT_phonoPAT <-  paste(helper, semPAC_phonoPAC_semPAT_phonoPAT, optimize, sep = '\n')

model_freq_len_model_semPAC_phonoPAC_semPAT_phonoPAT <- paste(helper, freq_len_semPAC_phonoPAC_semPAT_phonoPAT, optimize, sep = '\n')


#Make a list of theses models
model_list = list(model_semPAT, model_semPAC, model_phonoPAT, model_phonoPAC,
                  model_semPAC_phonoPAC_semPAT_phonoPAT,
                  model_freq_len_model_semPAC_phonoPAC_semPAT_phonoPAT
                  #model_semPAC_semPAT, 
                  #model_phonoPAC_phonoPAT,
                  #model_semPAC_phonoPAC,
                  #model_semPAT_phonoPAT,
                  #model_semPAT_phonoPAC,
                  #model_semPAC_phonoPAT
                  )

model_dict = data.frame(model_index=seq(1, length(model_list)),
                           model_name=c('semPAT', 'semPAC', 'phonoPAT', 'phonoPAC',
                                        'model_semPAC_phonoPAC_semPAT_phonoPAT',
                                        'model_freq_len_semPAC_phonoPAC_semPAT_phonoPAT'
                                        #'model_semPAC_semPAT', 
                                        #'model_phonoPAC_phonoPAT',
                                        #'model_semPAC_phonoPAC',
                                        #'model_semPAT_phonoPAT',
                                        #'model_semPAT_phonoPAC',
                                        #'model_semPAC_phonoPAT'
                                        ))

```

Compute the posterior (looping over languages and models)

```{r}

posteriors_all <- vector("list", length = length(data_scaled))

for (j in 1:length(data_scaled)) {

  posteriors_lang <- vector("list", length = length(model_list))
  
  for (i in 1:length(model_list)) {
    
    posterior_model <- NULL
    
    while (typeof(posterior_model)=="NULL") {
      posterior_model <- webppl(model_list[[i]],
                                data = data_scaled[[j]],
                                data_var = "data",
                                inference_opts = list(
                                  method = "MCMC",
                                  samples = 200, 
                                  burn = 50),
                                model_var = "AOA",
                                output_format = "webppl")
    }
    
    posteriors_lang[[i]] <- posterior_model
  }
  
  posteriors_all[[j]] <- posteriors_lang

}

#feather::write_feather(posteriors_all, "saved_data/posterior.feather")
```

Model for all data combined
```{r}

posteriors_combined <- vector("list", length = length(model_list))

for (i in 1:length(model_list)) {
    
    posterior_model <- NULL
    
    while (typeof(posterior_model)=="NULL") {
      
      posterior_model <- webppl(model_list[[i]],
                                data = data_combined,
                                data_var = "data",
                                inference_opts = list(
                                  method = "MCMC",
                                  samples = 2000, 
                                  burn = 500),
                                model_var = "AOA",
                                output_format = "webppl")
    }
    
    posteriors_combined[[i]] <- posterior_model
  

}

#feather::write_feather(posteriors_all, "saved_data/posterior.feather")
```


Process parameters

```{r}

parameters_all = data.frame()

for (j in 1:length(posteriors_all)) {
  
  parameters_lang = data.frame()
  
  for (i in 1:length(posteriors_all[[j]])) {
    
    param <- posteriors_all[[j]][[i]] %>%
      
      select(starts_with("value.parameters")) %>%
      gather(alpha, value, value.parameters.alpha1:value.parameters.alpha4) %>%
      mutate(alpha = gsub("value.parameters.", "", alpha)) %>%
      mutate(model_index = i) #%>%
      #left_join(model_dict)
    
    parameters_lang <- bind_rows(parameters_lang, param)
  }
  
   parameters_lang <- parameters_lang %>%
     mutate(language = j ) 
   
  parameters_all <- bind_rows(parameters_all, parameters_lang) 
  
}
 
parameters_all <- parameters_all %>%
  left_join(model_dict) %>%
  left_join(language_dict) %>%
  select(-model_index, -language) %>%
  rename(language = language_name,
         model = model_name)

feather::write_feather(parameters_all, "saved_data/parameters_all.feather")
parameters_all <- feather::read_feather("saved_data/parameters_all.feather")
```

Compute parameters for data combined:

```{r}
parameters_combined = data.frame()
  
  for (i in 1:length(posteriors_combined)) {
    
    param <- posteriors_combined[[i]] %>%
      
      select(starts_with("value.parameters")) %>%
      gather(alpha, value, value.parameters.alpha1:value.parameters.alpha4) %>%
      mutate(alpha = gsub("value.parameters.", "", alpha)) %>%
      mutate(model_index = i) #%>%
      #left_join(model_dict)
    
    parameters_combined <- bind_rows(parameters_combined, param)
}
  
 
parameters_combined <- parameters_combined %>%
  left_join(model_dict) %>%
  rename(model= model_name)
```


Compute point estimate and credible intervals
```{r}

parameters_sum <- parameters_all %>%
  group_by(language, model, alpha) %>%
  summarise(mean = mean(value),
            median = median(value),
            quantile_lower = quantile(value, c(0.025, 0.975))["2.5%"],
            quantile_upper = quantile(value, c(0.025, 0.975))["97.5%"],
            hdi_lower = hdi(value)['lower'],
            hdi_upper = hdi(value)['upper']
  )

  
post_combined <- parameters_sum %>%
  filter(model == "model_semPAC_phonoPAC_semPAT_phonoPAT") %>%
  ungroup() %>%
  select(-model,  -median, -quantile_lower, -quantile_upper) %>%
  rename(model = alpha)

post_combined$model <- mapValues(post_combined$model, 
                                 from = c("alpha1", "alpha2","alpha3","alpha4"), 
                                 to = c("semPAC", "phonoPAC", "semPAT", "phonoPAT"))

post_combined$model <- factor(post_combined$model, levels = c("phonoPAC", "phonoPAT", "semPAC", "semPAT"))

feather::write_feather(post_combined, "saved_data/combined_net.feather")
post_combined <- feather::read_feather("saved_data/combined_net.feather")


post_individual <- parameters_sum %>%
  filter(model != "model_semPAC_phonoPAC_semPAT_phonoPAT",
         alpha == "alpha1") 

post_individual$model <- factor(post_individual$model, levels = c("semPAC", "phonoPAC","semPAT", "phonoPAT"))

##Prepare plot: join posterior on individual predictors and t-test in one data.frame
post_individual_toJoin <- post_individual %>%
  select(-alpha, -median, -quantile_lower, -quantile_upper) %>%
  rename(lower = hdi_lower,
         upper = hdi_upper) %>%
  mutate(Evaluation = 'Fit')
  

zscores_toJoin <- zscores
zscores_toJoin$model <- mapValues(zscores_toJoin$model, 
                                 from = c("PAC_assoc", "PAC_phono_t2", "PAT_assoc", "PAT_phono_t2"), 
                                 to = c("semPAC", "phonoPAC", "semPAT", "phonoPAT"))

zscores_toJoin <- zscores_toJoin %>%
  select(-parameter, -p.value) %>%
  mutate(Evaluation = 't-test')

individual_both <- bind_rows(post_individual_toJoin, zscores_toJoin)



feather::write_feather(individual_both, "saved_data/individual_both.feather")
individual_both <- feather::read_feather("saved_data/individual_both.feather")

######

ggplot(individual_both, aes(x = model , y = mean)) +
  geom_pointrange(aes(ymin = lower, ymax = upper, col = model, linetype=Evaluation),
                  position = position_dodge(width = .5),
                  size = 0.5,
                  fatten = 0.5)+
  geom_hline(yintercept = 0, color = "grey", linetype = "dashed")+
  coord_flip() +
  facet_wrap(~language, ncol=4)  +
  guides(colour=FALSE, linetype = guide_legend(override.aes = list(size=0.3)))+
  scale_colour_solarized() +
  theme_bw()+
  theme(aspect.ratio = 0.7)

ggplot(post_combined, aes(x = alpha, y = mean)) +
  geom_pointrange(aes(ymin = hdi_lower, ymax = hdi_upper, y = mean, col = alpha))+
  geom_hline(yintercept = 0, color = "grey", linetype = "dashed")+
  coord_flip() +
  facet_wrap(~language, ncol=4)  +
  guides(colour=FALSE)+
  scale_colour_solarized() +
  theme_bw()+
  theme(aspect.ratio = 0.7)

individ <- individual_both %>%
  


```