---
title: "Model Comparisons Analysis Code"
author: "Salomons et al., 2021, Current Biology"
date: "6/9/2021"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
## Loading libraries 
```{r}
library(readxl)
library(multcomp)
library(tidyverse)
library(ggpubr)
library(dplyr)
library(sjPlot)
library(sjlabelled)
library(sjmisc)
library(ggplot2)
library(broom)
library(stargazer)
library(tibble)
```
## Loading and setting up data

```{r results='hide', message=FALSE}
############ load in data
all_dat = read.csv("Summary_Data.csv")

all_dat$Species = factor(all_dat$Species,
                                level = c("Dog","Wolf"))
all_dat$Sex = factor(all_dat$Sex,
                         level = c("M","F"))

all_dat <- all_dat %>% 
  rowwise() %>%
  mutate(Total.Temperament=case_when(is.na(Unfamiliar.Human) ~ "NA",
                                     is.na(Familiar.Human) ~ "NA",
                                     is.na(Unfamiliar.Object) ~ "NA",
                                     is.na(Familiar.Object) ~ "NA",
                     TRUE ~ as.character(sum(Unfamiliar.Human,Familiar.Human,Unfamiliar.Object,Familiar.Object)))) #adding a column with the total temperament score - only individuals who completed all trials since it is a sum
                     
all_dat$Total.Temperament <- as.numeric(all_dat$Total.Temperament)

all_dat <- all_dat %>% 
  relocate(Total.Temperament, .before = Unfamiliar.Human )

all_dat <- all_dat %>%
  rowwise() %>% 
  mutate(Point.with.Temp = case_when(is.na(Total.Temperament) ~ "NA",
                                     TRUE ~ as.character(Pointing.Gesture)
                                     ))

all_dat$Point.with.Temp <- as.numeric(all_dat$Point.with.Temp)
all_dat <- all_dat %>% 
  relocate(Point.with.Temp, .before = Pointing.Gesture) #adding a column for the pointing gesture scores only for the individuals who also completed temperament (if didn't complete temperament, NA)

all_dat <- all_dat %>%
  rowwise() %>% 
  mutate(Mark.with.Temp = case_when(is.na(Total.Temperament) ~ "NA",
                                     TRUE ~ as.character(Marker.Gesture)
                                     ))

all_dat$Mark.with.Temp <- as.numeric(all_dat$Mark.with.Temp)
all_dat <- all_dat %>% 
  relocate(Mark.with.Temp, .before = Marker.Gesture) #adding a column for the pointing gesture scores only for the individuals who also completed temperament (if didn't complete temperament, NA)

```
## Comparing models with sex and age inluded as variables - Pointing Gesture

### Please note that all table outputs within the R console will appear as a string of html code, so we have coded to hide these results from the knitted file - to view the formatted tables, please open the .html file which will automatically be uploaded to the folder which is set as the working directory. 

```{r results='hide', message=FALSE}
lm_Point_Species <- lm(Pointing.Gesture ~ Species, data=all_dat)
lm_Point_Species$AIC <- AIC(lm_Point_Species)

lm_Point_Age <- lm(Pointing.Gesture ~ Gesture.Weeks, data=all_dat) 
lm_Point_Age$AIC <- AIC(lm_Point_Age) 
 
lm_Point_Species_Age <- lm(Pointing.Gesture ~ Species + Gesture.Weeks , data=all_dat) 
lm_Point_Species_Age$AIC <- AIC(lm_Point_Species_Age)

age_models_Point <- list(lm_Point_Species, lm_Point_Age, lm_Point_Species_Age)

lm_Point_Sex <- lm(Pointing.Gesture ~ Sex, data=all_dat) #only sex

lm_Point_Species_Sex <- lm(Pointing.Gesture ~ Species + Sex , data=all_dat)#species and sex
lm_Point_Species_Sex$AIC <- AIC(lm_Point_Species_Sex)

sex_models_Point <- list(lm_Point_Species, lm_Point_Sex, lm_Point_Species_Sex)

stargazer(lm_Point_Species, lm_Point_Species_Age, lm_Point_Species_Sex, 
          covariate.labels = c("Species = Wolf", "Age (Weeks)", "Sex = F", "Constant"), 
          dep.var.labels   = "Arm Pointing Models", 
          column.labels = c("Species Only", "Species + Age", "Species + Sex"), 
          model.numbers = FALSE, omit.table.layout = "l", type="html", out="Pointing_Age_Sex_Models.htm")
```
## Comparing models with sex and age inluded as variables - Marker Gesture

```{r results='hide', message=FALSE}
lm_Mark_Species <- lm(Marker.Gesture ~ Species, data=all_dat)
lm_Mark_Species$AIC <- AIC(lm_Mark_Species)
  
lm_Mark_Age <- lm(Marker.Gesture ~ Gesture.Weeks, data=all_dat) 
lm_Mark_Age$AIC <- AIC(lm_Mark_Age) 
 
lm_Mark_Species_Age <- lm(Marker.Gesture ~ Species + Gesture.Weeks , data=all_dat) 
lm_Mark_Species_Age$AIC <- AIC(lm_Mark_Species_Age)

age_models_Mark <- list(lm_Mark_Species, lm_Mark_Age, lm_Mark_Species_Age)

lm_Mark_Sex <- lm(Marker.Gesture ~ Sex, data=all_dat) 

lm_Mark_Species_Sex <- lm(Marker.Gesture ~ Species + Sex , data=all_dat)
lm_Mark_Species_Sex$AIC <- AIC(lm_Mark_Species_Sex)

sex_models_Mark <- list(lm_Mark_Species, lm_Mark_Sex, lm_Mark_Species_Sex)

stargazer(lm_Mark_Species, lm_Mark_Species_Age, lm_Mark_Species_Sex, 
          covariate.labels = c("Species = Wolf", "Age (Weeks)", "Sex = F", "Constant"), 
          dep.var.labels   = "Marker Models", 
          column.labels = c("Species Only", "Species + Age", "Species + Sex"), 
          model.numbers = FALSE, omit.table.layout = "l", type="html", out="Marker_Age_Sex_Models.htm")
```

## Comparing models with age and sex included as variables - total temperament score

```{r results='hide', message=FALSE}
lm_TotalTemp_Species <- lm(Total.Temperament ~ Species, data=all_dat)
lm_TotalTemp_Species$AIC <- AIC(lm_TotalTemp_Species)

lm_TotalTemp_Species_Age <- lm(Total.Temperament ~ Species + Temperament.Weeks, data=all_dat)
lm_TotalTemp_Species_Age$AIC <- AIC(lm_TotalTemp_Species_Age)

lm_TotalTemp_Sex <- lm(Total.Temperament ~ Sex, data=all_dat)

lm_TotalTemp_Species_Sex <- lm(Total.Temperament ~ Species + Sex, data=all_dat)
lm_TotalTemp_Species_Sex$AIC <- AIC(lm_TotalTemp_Species_Sex)

stargazer(lm_TotalTemp_Species, lm_TotalTemp_Species_Age, lm_TotalTemp_Species_Sex, 
          covariate.labels = c("Species = Wolf", "Age (Weeks)", "Sex = F", "Constant"), 
          dep.var.labels   = "Total Temperament Models", 
          column.labels = c("Species Only", "Species + Age", "Species + Sex"), 
          model.numbers = FALSE, omit.table.layout = "l", type="html", out="Total_Temperament_Age_Sex_Models.htm") 
```

## Comparing models with age and sex included as variables - working memory
```{r results='hide', message=FALSE}
lm_Memory_Species <- lm(Working.Memory ~ Species, data=all_dat)
lm_Memory_Species$AIC <- AIC(lm_Memory_Species)

lm_Memory_Species_Age <- lm(Working.Memory ~ Species + Working.Memory.Weeks, data=all_dat)
lm_Memory_Species_Age$AIC <- AIC(lm_Memory_Species_Age)

lm_Memory_Species_Sex <- lm(Working.Memory ~ Species + Sex, data=all_dat)
lm_Memory_Species_Sex$AIC <- AIC(lm_Memory_Species_Sex)

stargazer(lm_Memory_Species, lm_Memory_Species_Age, lm_Memory_Species_Sex, 
          covariate.labels = c("Species = Wolf", "Age (Weeks)", "Sex = F", "Constant"), 
          dep.var.labels   = "Working Memory Models", 
          column.labels = c("Species Only", "Species + Age", "Species + Sex"), 
          model.numbers = FALSE, omit.table.layout = "l", type="html", out="Memory_Age_Sex_Models.htm") 
```

## Comparing models with age and sex included as variables - inhibitory control

```{r results='hide', message=FALSE}
lm_Inhibitory_Species <- lm(Inhibitory.Control ~ Species, data=all_dat)
lm_Inhibitory_Species$AIC <- AIC(lm_Inhibitory_Species)

lm_Inhibitory_Species_Age <- lm(Inhibitory.Control ~ Species + Inhibitory.Control.Weeks, data=all_dat)
lm_Inhibitory_Species_Age$AIC <- AIC(lm_Inhibitory_Species_Age)

lm_Inhibitory_Species_Sex <- lm(Inhibitory.Control ~ Species + Sex, data=all_dat)
lm_Inhibitory_Species_Sex$AIC <- AIC(lm_Inhibitory_Species_Sex)

stargazer(lm_Inhibitory_Species, lm_Inhibitory_Species_Age, lm_Inhibitory_Species_Sex, 
          covariate.labels = c("Species = Wolf", "Age (Weeks)", "Sex = F", "Constant"), 
          dep.var.labels   = "Inhibitory Control Models", 
          column.labels = c("Species Only", "Species + Age", "Species + Sex"), 
          model.numbers = FALSE, omit.table.layout = "l", type="html", out="Inhibitory_Age_Sex_Models.htm") 
```
## Comparing models with age and sex included as variables - human eye contact 

```{r results='hide', message=FALSE}
lm_EyeContact_Species <- lm(Eye.Contact ~ Species, data=all_dat)
lm_EyeContact_Species$AIC <- AIC(lm_EyeContact_Species)

lm_EyeContact_Species_Age <- lm(Eye.Contact ~ Species + Eye.Contact.Weeks, data=all_dat)
lm_EyeContact_Species_Age$AIC <- AIC(lm_EyeContact_Species_Age)

lm_EyeContact_Species_Sex <- lm(Eye.Contact ~ Species + Sex, data=all_dat)
lm_EyeContact_Species_Sex$AIC <- AIC(lm_EyeContact_Species_Sex)

stargazer(lm_EyeContact_Species, lm_EyeContact_Species_Age, lm_EyeContact_Species_Sex, 
          covariate.labels = c("Species = Wolf", "Age (Weeks)", "Sex = F", "Constant"), 
          dep.var.labels   = "Human Eye Contact Models", 
          column.labels = c("Species Only", "Species + Age", "Species + Sex"), 
          model.numbers = FALSE, omit.table.layout = "l", type="html", out="EyeContact_Age_Sex_Models.htm") 

```


## Comparing models with total temperament score included as variable - pointing gesture 

```{r results='hide', message=FALSE}
#only include individuals who completed both pointing and temperament tests

lm_Point_Species_onlytemp <- lm(Point.with.Temp ~ Species, data=all_dat)
lm_Point_Species_onlytemp$AIC <- AIC(lm_Point_Species_onlytemp)

lm_Point_Temperament <- lm(Point.with.Temp ~ Total.Temperament, data=all_dat)
lm_Point_Temperament$AIC <- AIC(lm_Point_Temperament)

lm_Point_Species_Temperament <- lm(Point.with.Temp ~ Species + Total.Temperament, data=all_dat)
lm_Point_Species_Temperament$AIC <- AIC(lm_Point_Species_Temperament)

stargazer(lm_Point_Species_onlytemp,lm_Point_Temperament, lm_Point_Species_Temperament,
          covariate.labels = c("Species = Wolf","Total Temperament Score", "Constant"), 
          dep.var.labels   = "Arm Pointing Models", 
          column.labels = c("Species Only", "Temperament Only", "Species + Temperament"),
          model.numbers = FALSE, omit.table.layout = "l", type="html", out="Pointing_Temperament_Model.htm") 
```

## Comparing models with total temperament score included as variable - marker gesture 

```{r results='hide', message=FALSE}
#only include individuals who completed both marker and temperament tests

lm_Mark_Species_onlytemp <- lm(Mark.with.Temp ~ Species, data=all_dat)
lm_Mark_Species_onlytemp$AIC <- AIC(lm_Mark_Species_onlytemp)

lm_Mark_Temperament <- lm(Mark.with.Temp ~ Total.Temperament, data=all_dat)
lm_Mark_Temperament$AIC <- AIC(lm_Mark_Temperament)

lm_Mark_Species_Temperament <- lm(Mark.with.Temp ~ Species + Total.Temperament, data=all_dat)
lm_Mark_Species_Temperament$AIC <- AIC(lm_Mark_Species_Temperament)

stargazer(lm_Mark_Species_onlytemp,lm_Mark_Temperament, lm_Mark_Species_Temperament,
          covariate.labels = c("Species = Wolf","Total Temperament Score", "Constant"), 
          dep.var.labels   = "Marker Gesture Models", 
          column.labels = c("Species Only", "Temperament Only", "Species + Temperament"),
          model.numbers = FALSE, omit.table.layout = "l", type="html", out="Marker_Temperament_Model.htm") 
```

## Comparing models with trial number included as variable - pointing gesture 
```{r results='hide', message=FALSE}
Gesture_Dat = read.csv("Gesture_Data.csv")
Gesture_Dat$Species = factor(Gesture_Dat$Species,
                         level = c("Dog","Wolf"))
Gesture_Dat$Sex = factor(Gesture_Dat$Sex,
                     level = c("M","F"))
Gesture_Dat = Gesture_Dat[!is.na(Gesture_Dat$Correct),]

trials_pointing <- filter(Gesture_Dat, Condition == "proximal")
trials_marker <- filter(Gesture_Dat, Condition == "marker")
trials_bodyvspoint <- filter(Gesture_Dat, Condition == "asymmetrical")
trials_odor <- filter(Gesture_Dat, Condition == "odor")

logistic_mod_Point_Species = glm(Correct ~ Species,family = binomial,
                               data = trials_pointing)

logistic_mod_Point_Species_Trials = glm(Correct ~ Species + Trial,family = binomial,
                               data = trials_pointing)

stargazer(logistic_mod_Point_Species, logistic_mod_Point_Species_Trials,
          covariate.labels = c("Species = Wolf", "Trial Number", "Constant"), 
          dep.var.labels   = "Arm Pointing Models", 
          column.labels = c("Species Only", "Species + Trial"), 
          model.numbers = FALSE, omit.table.layout = "l", type="html", out="Pointing_Trial_Number_Models.htm")
          
```

## Comparing models with trial number included as variable - marker gesture 
```{r results='hide', message=FALSE}
logistic_mod_Mark_Species = glm(Correct ~ Species,family = binomial,
                               data = trials_marker)

logistic_mod_Mark_Species_Trials = glm(Correct ~ Species + Trial,family = binomial,
                               data = trials_marker)

stargazer(logistic_mod_Mark_Species, logistic_mod_Mark_Species_Trials,
          covariate.labels = c("Species = Wolf", "Trial Number", "Constant"), 
          dep.var.labels   = "Marker Gesture Models", 
          column.labels = c("Species Only", "Species + Trial"), 
          model.numbers = FALSE, omit.table.layout = "l", type="html", out="Marker_Trial_Number_Models.htm")
          
```


