RELIVAL/CODE/z_outdated/simulations.Rmd

---
title: "Simulation test"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)

library(simstudy)
library(lme4)
library("performance") # ICC
library(ggplot2)

RECALC=TRUE

data_sets = c('aclew', 'lena')


check_small_var<-function(x,y,i) round(x[i],3)==round(y[i],3) & round(y[i],3) == 0

fit_child_model<-function(dataframe, metric){
    # Fit formula where experiment is removed
    formula <- as.formula(paste0(metric, "~  (1|child_id)")) #removed age_s +
    model <- lmer(formula, data=dataframe)
    
    return (model)
}

extract_chi_variables<-function(model){
    icc.result.mixed <- c(icc(model)$ICC_adjusted,icc(model)$ICC_conditional)
    icc.result.split <- c(as.data.frame(icc(model, by_group=TRUE))$ICC, NA)
    
    ranefs_vars <- t(as.data.frame(VarCorr(model))["vcov"])
    ranefs_stdv <- t(as.data.frame(VarCorr(model))["sdcor"])
    # chi, NA, residual
    ranefs_vars <-c(ranefs_vars[1],NA,ranefs_vars[2])
    ranefs_stdv <-c(ranefs_stdv[1],NA,ranefs_stdv[2])
    # chi, NA
    ns<-c(unlist(summary(model)$n),NA)
    return (c(#coefficients(summary(model))["age_s",],
                icc.result.mixed,
                icc.result.split,
                ranefs_vars,
                ranefs_stdv,
                nobs(model),ns))
}

extract_full_variables<-function(model){
    icc.result.mixed <- c(icc(model)$ICC_adjusted,icc(model)$ICC_conditional)
    icc.result.split <- t(as.data.frame(icc(model, by_group=TRUE))$ICC)

    ranefs_vars <- t(as.data.frame(VarCorr(model))["vcov"])
    ranefs_stdv <- t(as.data.frame(VarCorr(model))["sdcor"])
      
    ns<-t(data.frame(summary(model)$n))
    
    return (c( # coefficients(summary(model))["age_s",],
                icc.result.mixed,
                icc.result.split,
                ranefs_vars,
                ranefs_stdv,
                nobs(model),ns))
}

new_fit_models<-function(dataframe, data_set, metric,  fit_full = TRUE){ #age,
  iqr = quantile(dataframe[,metric],.75,na.rm=T)-quantile(dataframe[,metric],.25,na.rm=T)
  
  if(fit_full){
    # Fit full model
    formula <- as.formula(paste0(metric, "~  (1|experiment/child_id)")) #removed age_s +
    model <- lmer(formula, data=dataframe)
  }
  
  
  if(fit_full & !isSingular(model)) # Fitted full model
  {
    form="full"
    
    sw=shapiro.test(resid(model))$p
    mod_variables <- extract_full_variables(model)
    
    # Build line
  } else {
    model <- fit_child_model(dataframe, metric)
    if(!isSingular(model)){
      
      form = "no_exp"
      sw=shapiro.test(resid(model))$p
      mod_variables <- extract_chi_variables(model)
      
    } else {
      form='no_chi_effect'
      sw = NA
      mod_variables = c(
      #  NA,NA,NA, # c(coefficients(summary(model))["age_s",]
        NA,NA,    # icc.result.mixed
        NA,NA,    # icc.result.split
        NA,NA,NA, # ranefs_vars
        NA,NA,NA, # raners_std
        NA,NA,NA) # nobs (child, corpus), nobs
    }
  }
  
  icc.row = c(data_set,  metric, iqr, mod_variables, form, sw) #age,
  return (icc.row)
}

```


## Assumptions behind the simulation

We will inform the simulation by the data we have as follows:

- we'll have the same N of corpora, and of children in each corpus
- we'll have the same variables for each -- and these variables will have the same mean & SD for day 1 of recordings as in observed data

We'll depart from reality as follows:

- we will not consider the r across multiple days observed in the data, but instead generate data points to vary r from a small correlation (r=.1), a moderate one (r=.3), and a larger one (r=.5) -- higher values of r do not seem reasonable given what we know about infant measures
- we will not consider child age, nor variable re-recording periods
- we will have a single pair of recordings (rather than variable N of re-recordings)

## Implementation approach

We use simstudy, a package created for such simulations, following the vignette https://cran.r-project.org/web/packages/simstudy/vignettes/correlated.html to create correlated data providing a correlation matrix



```{r extract-parameters, eval=RECALC}

# Columns that should not be scaled or taken into account as metrics
no.scale.columns <- c('experiment', 'session_id', 'child_id','child_id_unique','age_s',
                      'date_iso', 'child_dob', 'missing_audio',"age_bin","duration","usession_id",
                      "normative","age","duration_alice", "duration_vcm" ,  "duration_vtc","duration_its" )


#create matrix to hold info in
df.vars.cols = c("data_set","experiment","n","metric", "mean","sd")
df.vars = data.frame(matrix(ncol=length(df.vars.cols),nrow=0, dimnames=list(NULL, df.vars.cols)), stringsAsFactors=FALSE)

for (data_set in data_sets){   # data_set = "aclew"
  mydat <- read.csv(paste0('../data_output/', data_set,'_metrics_scaled.csv')) # TO DISCUSS: scaled or unscaled?

  metrics = colnames(mydat)[!is.element(colnames(mydat), no.scale.columns)]
  
  #select down to first recording by child
  mydat$uchild_id=paste(mydat$experiment,mydat$child_id)
  mydat$child_id_age=paste(mydat$experiment,mydat$child_id,mydat$age)
  mydat=mydat[order(mydat$experiment,mydat$child_id,mydat$age),] #sort by child ID & age
  mydat=mydat[!duplicated(mydat$uchild_id),] #keep only the first line for each child
 
  
  means=stack(aggregate(mydat[,metrics],by=list(mydat$experiment),mean,na.rm=T)[,-1])
  sds=stack(aggregate(mydat[,metrics],by=list(mydat$experiment),sd,na.rm=T)[-1])
  
  df.vars=rbind(df.vars,
                cbind(data_set,levels(factor(mydat$experiment)),data.frame(table(mydat$experiment))$Freq,means[,c(2,1)],sds[,-2]))
    

}
colnames(df.vars)<-df.vars.cols
```


```{r generate-data, eval=RECALC}

alldays=NULL

for(i in 1:nrow(df.vars)) for(myr in c(.1,.3,.5, .7, .9)){#i=2;myr=.5

    #use a while loop to make sure data generated are close to the target r 
    C <- matrix(c( 1, myr,myr,1), nrow=2)
    simulation_unsatisfactory <- TRUE
    while(simulation_unsatisfactory==TRUE){
      try.this  <- as.data.frame(
        genCorData(df.vars$n[i], mu = c(df.vars$mean[i],df.vars$mean[i]), sigma = df.vars$sd[i], corMatrix = C) )
      try.this <- try.this[,-1] #remove the ID column, since it's useless
      sim.cor <- cor.test(try.this$V1,try.this$V2)$estimate
      simulation_unsatisfactory = !(abs(sim.cor-myr)<0.01)
    }
    thisdat=cbind(df.vars$data_set[i],df.vars$experiment[i],1:df.vars$n[i],stack(try.this),as.character(df.vars$metric[i]),myr)
    colnames(thisdat)<-c("data_set","experiment","child_id","value","day","metric","myr")
   alldays=rbind(alldays,thisdat)

}#end i

write.csv(alldays,"../output/simulated_correlated_2day.csv",row.names=F)

```


## Create ICC dataframe

```{r, eval=RECALC}
read.csv("../output/simulated_correlated_2day.csv")->alldays

df.icc.mixed.cols = c("data_set", "metric", "iqr", # removed "age_bin",
                      #"age_b","age_se","age_t", # beta, standard error, T #removed age
                      "icc_adjusted", "icc_conditional", 
                      "icc_child_id", "icc_corpus",
                      "child_id_var","corpus_var","residual_var",
                      "child_id_sd","corpus_sd","residual_sd",
                      "nobs","nchi", "ncor",
                      "formula","sw","myr")

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

for(myr in levels(factor(alldays$myr)))for(data_set in levels(factor(alldays$data_set))) for(metric in levels(factor(alldays$metric[alldays$data_set==data_set]))) {   # myr=.5 ; data_set = "aclew"; metric="wc_adu_ph"
  #select data
  mydat=alldays[alldays$myr==myr & alldays$data_set==data_set & alldays$metric==metric,] 
  #reshape mydat so that it fits expectation from the following function
  colnames(mydat)[4]<-metric
    icc.row <- new_fit_models(mydat, data_set, metric,  TRUE) #removed age NA,
    df.icc.mixed[nrow(df.icc.mixed) + 1,] <- cbind(icc.row,myr)
  
}
write.csv(df.icc.mixed,"../output/df.icc.simu.csv",row.names=F)


```




## Analyze ICCs

What is the relationship between ICC and r? It looks like they are similar, but ICC undershoots. But just a little - nothing like squaring.

```{r}
read.csv("../output/df.icc.simu.csv")->df.icc.mixed

plot(df.icc.mixed$icc_child_id~df.icc.mixed$myr)

```

What is the relationship between iqr & ICC? None

```{r}


ggplot(df.icc.mixed, aes(iqr,icc_child_id, color=myr)) + geom_point()

```

What is the relationship between iqr & SW? None

```{r}


ggplot(df.icc.mixed, aes(icc_child_id, sw, color=myr)) + geom_point()

```