log_india <- read.csv(here("data/original", "india_master.csv")) %>% 
  select(id, grade) %>%
  rename(subid = id) %>%
  mutate(age = grade + 5) %>%
  select(-grade)

log_korea <- read.csv(here("data/original", "korea.csv")) %>% 
  distinct(subid, .keep_all=T) %>%
  select(subid, age)

log_iran <- read.csv(here("data/original", "iran.csv")) %>%
  distinct(subid, .keep_all=T) %>%
  select(subid, age)

log_canada <- read.csv(here("data/original", "Canada with exact age.csv")) %>%
  rename(subid = sub_id,
         age = Age) %>%
  distinct(subid, .keep_all=T) %>%
  select(subid, age)

log_all <- rbind(log_india, log_korea, log_iran, log_canada)
  
d <- left_join(d, log_all)
## Joining, by = "subid"
## Warning: Column `subid` joining factors with different levels, coercing to
## character vector
summary(d)
##     subid                agebin               task               item     
##  Length:9519        grade3-4:3919   moral       :4080   namecall   :1360  
##  Class :character   grade5-6:4144   conventional:5439   push       :1360  
##  Mode  :character   grade8-9:1456                       rip        :1360  
##                                                         shoes      :1359  
##                                                         swim       :1360  
##                                                         teachername:1360  
##                                                         toy        :1360  
##       q_kind         answer           site           age        
##  base_q  :2379   Min.   :0.0000   Iran  :2184   Min.   : 5.000  
##  faraway :2380   1st Qu.:0.0000   India :4144   1st Qu.: 8.000  
##  no rule :2380   Median :0.0000   Korea : 952   Median : 9.680  
##  everyone:2380   Mean   :0.5357   Canada:2239   Mean   : 9.541  
##                  3rd Qu.:1.0000                 3rd Qu.:11.000  
##                  Max.   :6.0000                 Max.   :14.000  
##                  NA's   :578                    NA's   :252
d %>%
  group_by(site, agebin, subid) %>%
  summarise(n=n()) %>%
  group_by(site, agebin) %>%
  summarise(n=n())
## # A tibble: 9 x 3
## # Groups:   site [?]
##     site   agebin     n
##   <fctr>   <fctr> <int>
## 1   Iran grade3-4    57
## 2   Iran grade5-6    21
## 3  India grade3-4    17
## 4  India grade5-6    78
## 5  India grade8-9    52
## 6  Korea grade3-4    20
## 7  Korea grade5-6    14
## 8 Canada grade3-4    46
## 9 Canada grade5-6    34
ggplot(data=d %>%
         group_by(site, agebin, subid) %>%
         summarise(n = n()) %>%
         group_by(site, agebin) %>%
         summarise(number = n()) %>%
         filter(!is.na(agebin)), 
       
       aes(x=site, y=number, fill=agebin)) +
  geom_bar(position="dodge", stat="identity") +
  theme_few() +
  scale_fill_ptol() +
  ylab("number of participants")

# how bad is ~?
d %>% 
  filter(q_kind == "base_q", !is.na(answer)) %>%
  group_by(site, agebin, task) %>%
  multi_boot_standard(col="answer") %>%
  ggplot(., aes(x=agebin, y=mean, fill=task, col=task, group=task)) +
  # geom_bar(position="dodge", stat="identity") +
  # geom_errorbar(aes(ymin=ci_lower,ymax=ci_upper), position="dodge") +
  geom_line() +
  geom_linerange(aes(ymin = ci_lower, ymax = ci_upper), 
               position = position_dodge(width = .1)) +
  geom_hline(yintercept = 3, lty = 2) +
  facet_grid(.~site) + 
  theme_few() +
  scale_color_ptol()

d %>%
  mutate(age_group = floor(age)) %>%
  filter(q_kind == "base_q") %>%
  group_by(site, age_group, task, subid, item) %>%
  summarise(answer = mean(answer)) %>%
  group_by(site, age_group, task) %>%
  multi_boot_standard(col="answer", na.rm=T) %>%
  ggplot(., aes(x=age_group, y=mean, col=task)) +
  geom_line() +
  geom_linerange(aes(ymin = ci_lower, ymax = ci_upper), 
               position = position_dodge(width = .1)) +
  facet_grid(.~site) +
  theme_few() +
  scale_color_ptol()
## Warning: Removed 2 rows containing missing values (geom_path).
## Warning: Removed 8 rows containing missing values (geom_linerange).

ggplot(filter(d, q_kind=="base_q", !is.na(age), answer!="0"), 
       aes(x=age, y=answer, col=task)) +
  geom_jitter(width=.3, alpha=.3)+
  facet_grid(.~site) +
  geom_smooth() +
  theme_few() +
  scale_color_ptol()+
  geom_hline(yintercept =3, lty=2)
## `geom_smooth()` using method = 'loess'
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 9.98
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 3.02
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 2.9374e-16
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 9.1204
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used
## at 9.98
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 3.02
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal
## condition number 2.9374e-16
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other
## near singularities as well. 9.1204
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 9.98
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 3.02
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 3.5395e-15
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 9.1204
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used
## at 9.98
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 3.02
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal
## condition number 3.5395e-15
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other
## near singularities as well. 9.1204

summary(lmer(data=d %>%
               mutate(age = scale(age)), 
             answer ~ site * task * age + (task | subid) + (site + age | item)))
## Linear mixed model fit by REML ['lmerMod']
## Formula: answer ~ site * task * age + (task | subid) + (site + age | item)
##    Data: d %>% mutate(age = scale(age))
## 
## REML criterion at convergence: 21707.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.5180 -0.6098 -0.4168  0.5457  6.7819 
## 
## Random effects:
##  Groups   Name             Variance Std.Dev. Corr                   
##  subid    (Intercept)      0.010730 0.10359                         
##           taskconventional 0.004604 0.06785  1.00                   
##  item     (Intercept)      0.004443 0.06666                         
##           siteIndia        0.002669 0.05166   0.33                  
##           siteKorea        0.002104 0.04587   0.91 -0.08            
##           siteCanada       0.005980 0.07733   0.76 -0.36  0.96      
##           age              0.001496 0.03868   1.00  0.24  0.95  0.82
##  Residual                  0.684548 0.82737                         
## Number of obs: 8717, groups:  subid, 312; item, 7
## 
## Fixed effects:
##                                  Estimate Std. Error t value
## (Intercept)                      0.422298   0.059507   7.097
## siteIndia                        0.001202   0.069248   0.017
## siteKorea                       -0.012397   0.074030  -0.167
## siteCanada                       0.157479   0.076988   2.045
## taskconventional                 0.071390   0.075934   0.940
## age                              0.067980   0.050018   1.359
## siteIndia:taskconventional       0.104241   0.087010   1.198
## siteKorea:taskconventional       0.121292   0.093227   1.301
## siteCanada:taskconventional      0.145702   0.097722   1.491
## siteIndia:age                   -0.060651   0.058419  -1.038
## siteKorea:age                    0.010045   0.117246   0.086
## siteCanada:age                  -0.047122   0.063888  -0.738
## taskconventional:age             0.035911   0.062923   0.571
## siteIndia:taskconventional:age  -0.027893   0.072498  -0.385
## siteKorea:taskconventional:age   0.091923   0.146802   0.626
## siteCanada:taskconventional:age -0.033484   0.079313  -0.422
## 
## Correlation matrix not shown by default, as p = 16 > 12.
## Use print(x, correlation=TRUE)  or
##   vcov(x)     if you need it
# how bad is ~?
d %>% 
  filter(q_kind == "base_q", !is.na(answer)) %>%
  group_by(site, agebin, task, item) %>%
  multi_boot_standard(col="answer") %>%
  ungroup() %>%
  arrange(site, agebin, mean) %>%
  mutate(order = row_number()) %>%
  ggplot(., aes(x=order, y=mean, col=task, label=item)) +
  # geom_bar(position="dodge", stat="identity") +
  geom_text(size=3) +
  geom_linerange(aes(ymin=ci_lower,ymax=ci_upper), position="dodge", alpha=.3) +
  facet_wrap(site~agebin, ncol=4, scales="free") + 
  theme_few() +
  scale_color_ptol() +
  xlab("items") +
  ylab("mean rating (1=really bad)") +
  theme(legend.position = "bottom",
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank())
## Warning: Width not defined. Set with `position_dodge(width = ?)`

# would it be okay if ~?
d %>%
  filter(q_kind != "base_q", !is.na(answer)) %>%
  mutate(q_kind = fct_relevel(q_kind, "faraway", "no rule", "everyone")) %>%
  group_by(site, agebin, task, q_kind) %>%
  multi_boot_standard(col="answer") %>%
  ggplot(., aes(x=site, y=mean, fill=task)) +
  geom_bar(position="dodge", stat="identity") +
  geom_errorbar(aes(ymin=ci_lower,ymax=ci_upper), position="dodge") +
  facet_grid(q_kind~agebin) + 
  theme_few() +
  scale_fill_ptol() +
  ylab("mean proportion saying okay \ngiven new circumstances")

# would it be okay if ~?
d %>%
  filter(q_kind != "base_q", !is.na(answer)) %>%
  mutate(q_kind = fct_relevel(q_kind, "faraway", "no rule", "everyone"),
         age_group = floor(age)) %>%
  group_by(site, age_group, task, q_kind) %>%
  multi_boot_standard(col="answer") %>%
  ggplot(., aes(x=age_group, y=mean, col=task)) +
  geom_line() +
  geom_linerange(aes(ymin = ci_lower, ymax = ci_upper), 
               position = position_dodge(width = .1)) +
  facet_grid(q_kind~site) +
  theme_few() +
  scale_color_ptol() +
  ylab("mean proportion saying okay \ngiven new circumstances")
## Warning: Removed 2 rows containing missing values (geom_path).
## Warning: Removed 12 rows containing missing values (geom_linerange).

d %>% 
  filter(q_kind == "faraway", !is.na(answer)) %>%
  group_by(site, agebin, task, item) %>%
  multi_boot_standard(col="answer") %>%
  ungroup() %>%
  arrange(site, agebin, mean) %>%
  mutate(order = row_number()) %>%
  ggplot(., aes(x=order, y=mean, col=task, label=item)) +
  # geom_bar(position="dodge", stat="identity") +
  geom_text(size=3) +
  geom_linerange(aes(ymin=ci_lower,ymax=ci_upper), position="dodge", alpha=.3) +
  facet_wrap(site~agebin, ncol=4, scales="free") + 
  theme_few() +
  scale_color_ptol() +
  xlab("items") +
  ylab("mean proportion saying okay \nin faraway country") +
  theme(legend.position = "bottom",
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank())
## Warning: Width not defined. Set with `position_dodge(width = ?)`

d %>% 
  filter(q_kind == "no rule", !is.na(answer)) %>%
  group_by(site, agebin, task, item) %>%
  multi_boot_standard(col="answer") %>%
  ungroup() %>%
  arrange(site, agebin, mean) %>%
  mutate(order = row_number()) %>%
  ggplot(., aes(x=order, y=mean, col=task, label=item)) +
  # geom_bar(position="dodge", stat="identity") +
  geom_text(size=3) +
  geom_linerange(aes(ymin=ci_lower,ymax=ci_upper), position="dodge", alpha=.3) +
  facet_wrap(site~agebin, ncol=4, scales="free") + 
  theme_few() +
  scale_color_ptol() +
  xlab("items") +
  ylab("mean proportion saying okay \nif there's no rule") +
  theme(legend.position = "bottom",
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank())
## Warning: Width not defined. Set with `position_dodge(width = ?)`

d %>% 
  filter(q_kind == "everyone", !is.na(answer)) %>%
  group_by(site, agebin, task, item) %>%
  multi_boot_standard(col="answer") %>%
  ungroup() %>%
  arrange(site, agebin, mean) %>%
  mutate(order = row_number()) %>%
  ggplot(., aes(x=order, y=mean, col=task, label=item)) +
  # geom_bar(position="dodge", stat="identity") +
  geom_text(size=3) +
  geom_linerange(aes(ymin=ci_lower,ymax=ci_upper), position="dodge", alpha=.3) +
  facet_wrap(site~agebin, ncol=4, scales="free") + 
  theme_few() +
  scale_color_ptol() +
  xlab("items") +
  ylab("mean proportion saying okay \nif everyone did it") +
  theme(legend.position = "bottom",
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank())
## Warning: Width not defined. Set with `position_dodge(width = ?)`