# This code can be reproduced without access to the underlying datasets, unlike regenerate_data.R
# It relies on packages, functions, & variables that are called in SM.Rmd
# You may need to "get" the data -- so if you get
  #Error in file(file, "rt") : cannot open the connection
# then do datalad get FILE on terminal to get the relevant files

# Create main ICC including all corpora & all measurements
df.icc.mixed.cols = c("data_set","age_bin", "metric", "iqr",
                      "age_b","age_se","age_t", # beta, standard error, T
                      "icc_adjusted", "icc_conditional", 
                      "icc_child_id", "icc_corpus",
                      "child_id_var","corpus_var","residual_var",
                      "child_id_sd","corpus_sd","residual_sd","r2_cond", "r2_marg",
                      "nobs","nchi", "ncor",
                      "formula","sw")

df.icc.mixed = data.frame(matrix(ncol=length(df.icc.mixed.cols),nrow=0, dimnames=list(NULL, df.icc.mixed.cols)),
                          stringsAsFactors = FALSE)


for (data_set in data_sets){   # data_set = "aclew"
  # Load data
  mydat <- read.csv(paste0('../data_output/', data_set,'_metrics_scaled.csv'))
  metrics <- colnames(mydat)[!is.element(colnames(mydat), no.scale.columns)]
  
  for(metric in metrics) 
  {  # metric = "voc_chi_ph"
    icc.row <- new_fit_models(mydat, data_set, metric, NA, TRUE)
    df.icc.mixed[nrow(df.icc.mixed) + 1,] <- icc.row
  } 
}

write.csv(df.icc.mixed,"../output/df.icc.mixed.csv",row.names=F)

print("done with overall ana, tackling within corpus version")

# repeat for the version within each corpus

df.icc.corpus.cols = c("data_set","corpus","metric", "iqr",
                       "age_b","age_se","age_t",
                       "icc_adjusted", "icc_conditional",
                       "icc_child_id", "icc_corpus",
                       "child_id_var","corpus_var","residual_var",
                       "child_id_sd","corpus_sd","residual_sd","r2_cond", "r2_marg",
                       "nobs","nchi", "ncor",
                       "formula","sw")
df.icc.corpus = data.frame(matrix(ncol=length(df.icc.corpus.cols),nrow=0, dimnames=list(NULL, df.icc.corpus.cols)),
                           stringsAsFactors=FALSE)


for (data_set in data_sets){   # data_set = "aclew"
  # Load data 
  mydat <- read.csv(paste0('../data_output/', data_set,'_metrics_scaled.csv'))
  for(corpus in corpora){
    mycordat <- mydat[mydat$experiment == corpus, ]
    metrics <- colnames(mycordat)[!is.element(colnames(mycordat), no.scale.columns)] #this should be defined here bc it varies by data_set
    for(metric in metrics) {  # metric = "avg_voc_dur_mal"
      icc.row <- new_fit_models(mycordat, data_set, metric, corpus, FALSE)
      df.icc.corpus[nrow(df.icc.corpus) + 1,] <- icc.row
    } 
  }
}


write.csv(df.icc.corpus,"../output/df.icc.corpus.csv",row.names=F)

print("done with within corpus ana, tackling within age group version")

#We do this one separately because we want to standardize metrics within each age group

# repeat within age group bins
df.icc.age.cols =  c("data_set","age_bin","metric", "iqr",
                       "age_b","age_se","age_t",
                       "icc_adjusted", "icc_conditional",
                       "icc_child_id", "icc_corpus",
                       "child_id_var","corpus_var","residual_var",
                       "child_id_sd","corpus_sd","residual_sd","r2_cond", "r2_marg",
                       "nobs","nchi", "ncor",
                       "formula","sw")

df.icc.age = data.frame(matrix(ncol=length(df.icc.age.cols),nrow=0, dimnames=list(NULL, df.icc.age.cols)),
                        stringsAsFactors=FALSE)


for (data_set in data_sets){   # data_set = "aclew"
  # Load data and calculate age cuts
  mydat <- read.csv(paste0('../data_output/', data_set,'_metrics.csv')) # /!\ Do not use scaled version -> we'll scale by age later
  mydat$age_bin <- cut(mydat$age,c(0:6*6))
  
  metrics = colnames(mydat)[!is.element(colnames(mydat), no.scale.columns)]
  for(thisage in levels(mydat$age_bin))
  {   #thisage= "(0,6]"
    thiscordata <- mydat[mydat$age_bin == thisage,]
    for(metric in metrics) 
    {  # metric = "avg_voc_dur_mal"
      
      pre_scaled_metric <- (thiscordata[, metric] - mean(thiscordata[, metric], na.rm=T))/sd(thiscordata[, metric], na.rm=T)
      thiscordata[abs(pre_scaled_metric)>2.5 & !is.na(pre_scaled_metric), metric] <- NA
      thiscordata[, metric] <- (thiscordata[, metric] - mean(thiscordata[, metric], na.rm=T))/sd(thiscordata[, metric], na.rm=T)
      if(dim(thiscordata)[1] > 30 & length(levels(factor(thiscordata$experiment))) > 1)
      {
        icc.row <- new_fit_models(thiscordata, data_set, metric, thisage, TRUE)

      }else{
        icc.row <- c(data_set,thisage,metric,iqr,
                     NA,NA,NA,
                     NA,NA,
                     NA,NA,
                     NA,NA,NA,
                     NA,NA,NA,
                     NA,NA,NA,
                     "not_enough_data",NA)
      }
      df.icc.age[nrow(df.icc.age) + 1,] <- icc.row
    } 
  }
}
write.csv(df.icc.age,"../output/df.icc.age.csv",row.names=F)