alignment_correlations.Rmd

---
title: "alignment_correlations"
author: "Joseph Denby"
date: "12/20/2018"
output: html_document
---

```{r}
library(knitr)
library(data.table)
library(tidyverse)
library(readr)
library(stringr)
library(DT)
library(tidytext)
library(RSQLite)
library(DBI)
library(here)
library(feather)
library(magrittr)
library(zoo)
library(RMySQL)
library(rstan)

options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)

setwd(here())
```

```{r Run Stan Model with PPVT}
long_present <- read_feather('data/ppvt_with_demos.feather')

liwc <- read_tsv("word_lists/liwc2007_converted.tsv", col_names = FALSE) %>%
  rename(word = X1, category = X2)

# subj_effects <- read_csv('subjs_effects_sum.csv')

# long_present <- long_present %>% 
#   left_join(tt_mlu, by=c('subject','session','person'))
# 
# write_feather(long_present,'data/ppvt_with_demos.feather')

# long_present <- long_present %>% 
#   left_join(subj_effects, by=c('subject'='subj')) %>% 
#   select(-age_years, -visit, -ppvt) 

## STAN DATA ##
MidAge <- median(long_present$session)
NumMarkers <- length(unique(liwc$category))
NumSubPops <- length(unique(long_present$person))
NumSpeakers <- length(unique(long_present$uid))
NumObservations <- dim(long_present)[1]
NumSex <- length(unique(long_present$sex))
NumMomEd <- length(unique(long_present$mother_education))
NumAges <- length(unique(long_present$session))
SpeakerSubPop <- long_present %>% 
  distinct(person, uid) %>%
  mutate(subpop =  as.numeric(as.factor(person))) %>%
  pull(subpop) # Child - 1; Parent - 2
SpeakerMomEd <- long_present %>% 
  mutate(mother_education = mother_education+1) %>% 
  distinct(uid, mother_education) %>%
  pull(mother_education)
SpeakerSex <- long_present %>% 
  distinct(uid, sex) %>%
  mutate(sex =  as.numeric(as.factor(sex))) %>%
  pull(sex) # Female - 1; Male - 2
SpeakerId <-long_present$uid
SpeakerAge <- long_present$session
MarkerType <- as.numeric(as.factor(long_present$category)) #see liwc_markers
NumUtterancesAB <- long_present$NumUtterancesAB
NumUtterancesNotAB <- long_present$NumUtterancesNotAB
CountsAB <- long_present$CountsAB
CountsNotAB <- long_present$CountsNotAB
StdDev <- .25

parentid <- long_present$parentid
childid <- long_present$childid
ppvt_vals <- long_present$ppvt
factor_intercepts <- long_present$intercept
factor_slopes <- long_present$slope
age_years <- long_present$age_years
mother_education <- long_present$mother_education
female <- long_present$female
# black <- long_present$black
# multi <- long_present$multi
# income_category <- long_present$income_category
# ppvt_slopes <- long_present$model
mlu <- long_present$mlu
types <- long_present$types
tokens <- long_present$tokens
###############

##############
###SCALING####
mlu <- as.vector(scale(mlu))
types <- as.vector(scale(types))
tokens <- as.vector(scale(tokens))
age_years <- as.vector(scale(age_years))
mother_education <- as.vector(scale(mother_education))
##############

# fit <- stan("stan_models/sessionalignment.stan",
#             chains = 1,
#             iter = 500)
# saveRDS(fit, "sessionestimates.rds")
# 

fit <- stan("stan_models/alignmenthierarchicaldemos.stan",
            chains = 1,
            iter = 500)

# saveRDS(fit, "hierppvt500iters.rds")

divergent <- get_sampler_params(fit, inc_warmup=FALSE)[[1]][,'divergent__']
```

```{r fit extraction}
# fit_hierarchical <- readRDS('stanfits/oldmodeldemos_short.rds')
# fit_hierarchical <- readRDS('newnewhierdemos200iters.rds')
# fit_hierarchical <- readRDS('hierppvt350iters.rds')
fit_hierarchical <- readRDS('hierppvt500iters.rds')


SE_UPPER <- .84
SE_LOWER <- .16

# Define and extract all alignment-related parameters
pars <- names(fit_hierarchical) %>%
  as_data_frame() %>%
  filter(!str_detect(value, "Marker"), 
         !str_detect(value, "observation"),
         !str_detect(value, "mu_ab"),
         !str_detect(value, "mu_notab"),
         !str_detect(value, "lp__")) %>%
  pull()

parameter_cis <- rstan::extract(fit_hierarchical, pars) %>%
  bind_rows(.id = "sample") %>%
  gather(measure, value, -sample) %>%
  group_by(measure) %>%
  summarise(mean = mean(value), se_upper = quantile(value,SE_UPPER), 
            se_lower = quantile(value,SE_LOWER)) %>%
  ungroup() %>% 
  mutate(measure = gsub("eta_ab","eta-ab", measure)) %>%
  separate(measure, into = c("parameter", "type") , sep = "_", 
           extra = "merge") %>%
  mutate(id = gsub("[^0-9]*","", str_extract(type, "\\[[^()]+\\]")),
         type = gsub("\\[[^()]+\\]", "", type),
         type = factor(type, levels = c("pop", "subpop", "speaker")),
         id = case_when(
           is.na(id) ~ "pop",
           type == "subpop" & id == 1 ~ "child",
           type == "subpop" & id == 2 ~ "parent",
           T ~ id)
         ) %>%
  left_join(distinct(long_present, subject, uid, person) %>% 
              mutate(uid = as.character(uid)), by = c("id" = "uid")) %>%
  arrange(type, parameter) %>%
  mutate(id = if_else(is.na(subject), id, as.character(subject))) %>%
  select(-subject)

# Define and extract parameters related to PPVT prediction model
demopars <- names(fit_hierarchical) %>%
  as_data_frame() %>%
  filter(str_detect(value, "mu") | str_detect(value, 'ppvt') | str_detect(value, 'sigma') | str_detect(value, 'factor'),
         !str_detect(value, "observation"),
         !str_detect(value, "mu_ab"),
         !str_detect(value, "mu_notab")) %>%
  pull()

demoparameters <- extract(fit_hierarchical, demopars) %>%
  bind_rows(.id = "sample") %>%
  gather(measure, value, -sample) %>%
  group_by(measure) %>%
  summarise(mean = mean(value), se_upper = quantile(value,SE_UPPER), 
            se_lower = quantile(value,SE_LOWER))

```