---
title: "The multimodal nature of communicative efficiency in social interaction [main analysis]"
author: "Marlou Rasenberg"
date: '(this version: `r format(Sys.Date())`)'
output:
  html_document:
    toc: true
    toc_float: true
editor_options:
  chunk_output_type: console
---

```{r global_options, include=FALSE}
knitr::opts_chunk$set(echo=TRUE, warning=FALSE, message=FALSE, fig.path='4_out/')
```

```{r packages, results='hide',include=F}
# Packages
list.of.packages <- c("tidyr", "dplyr", "ggplot2","cowplot", "lmerTest","ggdist","ggh4x", "pastecs")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)>0) install.packages(new.packages)
lapply(list.of.packages, require, character.only=T)
rm(list.of.packages,new.packages)

```

This script supports the paper *The Multimodal Nature of Communicative Efficiency* by Marlou Rasenberg, Wim Pouw, Asli Özyürek and Mark Dingemanse (for the CABB team). To increase readability, not all code chunks present in the .Rmd source are shown in the output.

# Data
The 'raw' data for this study are the speech and gesture annotations (made in ELAN), as well the Kinect motion tracking data (see folder 1_data)
Those have been processed (see folder 2_processing), resulting in various dataframes usable for analyses (stored in folder 2_processing > processing_output).

The dataframes used in this script are the following:  
* df_trials (1_trials.csv): dataframe with *one trial per row* (contains information about time on task and accuracy)
* df (6_final_df_main.csv): the main dataframe with *one repair turn per row*  
* df_div_labour (6_final_df_division_labour): dataframe with *one repair sequence per row*  

Repair turns are labelled as follows:  
* Tmin1: trouble source turn  
* T0: repair initiation  
* Tplus1: repair solution

Repair sequences always consist of a T0 and Tplus1. Usually also a Tmin1 (but not if it the non-first T0 in a non-minimal sequence).

The first three dataframes are linked through 'ids':  
* repair_id: this refers to the repair turn (numbered 1 to 1065)  
* repair_seq_id: this refers to the repair sequence (numbered 1 to 378)

Note that the terms 'labour' and 'effort' are used interchangeably throughout.

# Intitial steps
```{r read data}
df <- read.csv2("2_processing/processing_output/6_final_df_main.csv", header=T, colClasses = c(pairnr="factor", round="factor", speech_original="character", speech_clean="character", onset="character", offset="character"), stringsAsFactors=T)

df_kinematics <- read.csv2("2_processing/processing_output/5_gestures_kinematics.csv", header=T, colClasses = c(pairnr="factor"), stringsAsFactors=T)

df_div_labour <- read.csv2("2_processing/processing_output/6_final_df_division_labour.csv", header=T, colClasses = c(pairnr="factor"), stringsAsFactors=T)

df_trials <- read.csv2("2_processing/processing_output/1_trials.csv", header=T, stringsAsFactors=T, colClasses = c(round="factor")) #pairnrs are '04' etc. so read in as numeric, then set as factor
df_trials$pairnr <- as.factor(df_trials$pairnr)

```

```{r process data}
#reorder levels
df$T0_type <- factor(df$T0_type, levels = c("openrequest", "restrictedrequest", "restrictedoffer"))
df$repair_part <- factor(df$repair_part, levels = c("Tmin1", "T0", "Tplus1"))

df_div_labour$T0_type <- factor(df_div_labour$T0_type, levels = c("openrequest", "restrictedrequest", "restrictedoffer"))
```

```{r stat desc selection}
#when using stat.desc to get descriptive values, show only min [4], max[5], mean [9] and std.dev [13]
stat_desc_output <- c(4,5,9,13)
```


# Descriptives

```{r task info}
#compute trial durations by dyad
df_trials$duration_msec <- df_trials$offset_msec - df_trials$onset_msec
df_ToT <- df_trials %>% group_by(pairnr) %>% summarize(total_duration_msec = sum(duration_msec), total_duration_min=sum(duration_msec)/60000)

stat.desc(df_ToT$total_duration_min)[stat_desc_output]

rm(df_ToT)
```


### repair

```{r repair frequency}
#overall
freq_T0_bypair <- df %>% subset(repair_part=="T0") %>% group_by(pairnr, .drop=FALSE) %>% count() 
stat.desc(freq_T0_bypair$n)[stat_desc_output]

#types
df %>% subset(repair_part=="T0") %>% group_by(T0_type) %>% count()
# 1 openrequest          24
# 2 restrictedrequest    39
# 3 restrictedoffer     315

#repair rate relative to time
ToT_bypair <- df_trials %>%
  group_by(pairnr) %>%
  summarise(ToT_msec = sum(duration_msec), ToT_min = sum(duration_msec)/60000)

repair_by_time <- merge(freq_T0_bypair, ToT_bypair, by="pairnr")
repair_by_time$repair_every_x_sec <- repair_by_time$ToT_min/repair_by_time$n

stat.desc(repair_by_time$repair_every_x_sec)[stat_desc_output]

rm(repair_by_time, ToT_bypair)

```

### speech (orthographic length)
```{r speech descriptives}
tapply(df$speech_count, list(df$T0_type, df$repair_part), mean)
tapply(df$speech_count, list(df$T0_type, df$repair_part), sd)
```

```{r orthographic length vs. with turn duration}

df_characters_vs_duration <- df %>% subset(repair_part!="Tmin1") %>% select(repair_id, speech_count, onset_msec, offset_msec)
df_characters_vs_duration$duration <- df_characters_vs_duration$offset_msec - df_characters_vs_duration$onset_msec

ggplot(data=df_characters_vs_duration, aes(x=speech_count, y=duration))+
         geom_point(alpha=.3)+
          geom_smooth(method="lm")+
         ylab("duration of repair turn in msec")+
         xlab("length of repair turn in orthographic characters")

       
cor.test(df_characters_vs_duration$speech_count, df_characters_vs_duration$duration)
```


### gesture counts
```{r gesture counts}
#total counts per type
df %>% subset(repair_part!="Tmin1") %>% summarize(sum_total=sum(total_count),sum_iconic=sum(iconic_count), sum_deictic=sum(deictic_count), sum_other=sum(other_count))
#percentages
df %>% subset(repair_part!="Tmin1") %>% summarize(prop_iconic=sum(iconic_count)/sum(total_count)*100, prop_deictic=sum(deictic_count)/sum(total_count)*100, prop_other=sum(other_count)/sum(total_count)*100)


#how many repair turns contain a gesture?
df %>% subset(repair_part!="Tmin1") %>% count(total_count!=0) %>% subset(`total_count != 0`==TRUE) %>% select(n) #281 (= 37.2% of all 756 T0/T+1 turns)

  #per type&part
  prop_gest_per_turn <- df %>% subset(repair_part!="Tmin1") %>% group_by(repair_part, T0_type)  %>% count(total_count!=0)
  prop_gest_per_turn %>% group_by(T0_type, repair_part, 'total_count != 0') %>% mutate(prop_gesture_in_turn = n/sum(n)) %>% select(repair_part, T0_type, n, prop_gesture_in_turn)

  #per sequence
  df_div_labour %>% subset(gesture_subm_T0!=0|gesture_subm_Tplus1!=0) %>% nrow()

```

### gesture submovements
```{r gesture kinematics descriptives}
tapply(df$sum_subm, list(df$T0_type,df$repair_part), mean)
tapply(df$sum_subm, list(df$T0_type, df$repair_part), sd)

```


### multimodal effort
```{r multimodal effort descriptives}
tapply( (df_div_labour$multi_T0 + df_div_labour$multi_Tplus1), df_div_labour$T0_type, mean)
tapply( (df_div_labour$multi_T0 + df_div_labour$multi_Tplus1), df_div_labour$T0_type, sd)
```


# Figures
```{r plot settings}
#colour scheme
colours_repair_part<- as.vector(unlist(inlmisc::GetColors(n = 7, scheme="vibrant")))

#use double space before \n for line break, ** for italics --> add ggtext::element_markdown()) to ggplot later on
repair_type_names <- list(
  'openrequest'="**open request**  \n(*n* = 24)",
  'restrictedrequest'="**restricted request**  \n(*n* = 39)",
  'restrictedoffer'="**restricted offer**  \n(*n* = 315)")

repair_type_labeller <- function(variable,value){return(repair_type_names[value])}


```

### division of speech, gesture and multimodal effort (Figure 1 is ms)
```{r division of speech and gesture effort [fig 1]}
p_speech_effort <- ggplot(data=subset(df, repair_part!="Tmin1"), aes(x=repair_part, y=speech_count, fill=repair_part, color=repair_part)) +
  facet_manual(~T0_type, design="ABC", widths = unit(3.5, "cm"))+
  facet_grid(~T0_type, labeller=repair_type_labeller, switch="both")+
  gghalves::geom_half_boxplot(alpha=.6, outlier.shape=NA, width=1, lwd=0.5, fatten=1, errorbar.length = 0.2, aes(colour = repair_part))+
  gghalves::geom_half_point(position= position_nudge(x=.35), size=.05, side = "l", range_scale = .3, alpha = .5, aes(colour = repair_part)) +
  scale_fill_manual(values=c(colours_repair_part[1],colours_repair_part[2]), labels=c("repair initiation", "repair solution"), name="repair_part")+
  scale_colour_manual(values=c(colours_repair_part[1],colours_repair_part[2]), labels=c("repair initiation", "repair solution"), name="repair_part")+
  expand_limits(x=c(0,0), y=c(0,400))+
  ggtitle("**Speech effort** (characters)")+ 
  labs(y="orthographic characters", x="")+
  #guides(color=guide_legend('repair_part',override.aes=list(size=2, alpha=1)))+
  theme_light()+
  theme(
        legend.title=element_blank(),
        legend.background=element_blank(), legend.key=element_rect(fill="#f6f6f6"), 
        legend.position=c(0.835,0.87),
        legend.key.size=unit(0.4, 'cm'),
        legend.text=element_text(size=7),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.background = element_rect(fill="#f6f6f6"),
        plot.title=ggtext::element_markdown(colour="black", size=10), #necessary to get only part in bold
        axis.line.y = element_line(colour="grey"),
        axis.title.y=element_text(size=7), axis.text.y=element_text(size=7),
        axis.text.x=element_blank(), axis.ticks.x=element_blank(),
        plot.margin = margin(5.5,0,5.5,5.5), #top, right, bottom, left+
        strip.background = element_blank(),
        strip.text.x = element_blank())

p_speech_effort
  

p_gesture_effort <- ggplot(data=subset(df, repair_part!="Tmin1"), aes(x=repair_part, y=sum_subm, fill=repair_part, color=repair_part)) +
  facet_manual(~T0_type, design="ABC", widths = unit(3.5, "cm"))+
  facet_grid(~T0_type, labeller=repair_type_labeller)+
  gghalves::geom_half_boxplot(alpha=.6, outlier.shape=NA, width=1, lwd=0.5, fatten=1, errorbar.length = 0.2)+
  gghalves::geom_half_point(position= position_nudge(x=.35), size=.05, side = "l", range_scale = .3, alpha = .5, aes(colour = repair_part)) +
  scale_x_discrete(breaks=c("T0","Tplus1"),labels=c("initiation", "solution"))+
  scale_fill_manual(values=c(colours_repair_part[1],colours_repair_part[2]))+
  scale_colour_manual(values=c(colours_repair_part[1],colours_repair_part[2]), labels=c("initiation", "solution"))+
  expand_limits(x=c(0,0), y=c(0,20))+
  ggtitle("**Gesture effort** (submovements)")+ 
  labs(y="gesture submovements", x="")+
  theme_light()+
  theme(
        legend.position="none",
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.background = element_rect(fill="#f6f6f6"),
        plot.title=ggtext::element_markdown(colour="black", size=10), #necessary to get only part in bold
        axis.line.y = element_line(colour="grey"),
        axis.title.y=element_text(size=7), axis.text.y=element_text(size=7),
        axis.text.x=element_blank(), axis.ticks.x=element_blank(),
        plot.margin = margin(5.5,0,5.5,5.5), #top, right, bottom, left+
        strip.background = element_blank(),
        strip.text.x = element_blank())
p_gesture_effort
```


```{r division of multimodal effort [fig 1]}
df_div_labour_prop_long <- df_div_labour %>% select(pairnr, target, repair_seq_id, T0_type, multi_T0_prop, multi_Tplus1_prop) %>% gather(repair_part, prop_multi, multi_T0_prop:multi_Tplus1_prop, factor_key=T)
df_div_labour_prop_long$repair_part <- gsub("multi_","",df_div_labour_prop_long$repair_part)



p_division_multi_labour_prop <- ggplot(df_div_labour_prop_long, aes(x=repair_part, y=prop_multi, fill=repair_part, color=repair_part)) +
  facet_manual(~T0_type, design="ABC", widths = unit(3.5, "cm"))+
  facet_grid(~T0_type, labeller=repair_type_labeller, switch="both")+
  gghalves::geom_half_boxplot(alpha=.6, outlier.shape=NA, width=1, lwd=0.5, fatten=1, errorbar.length = 0.2)+
  gghalves::geom_half_point(position= position_nudge(x=.35), size=.05, side = "l", range_scale = .3, alpha = .5, aes(colour = repair_part)) +
  scale_x_discrete(breaks=c("T0_prop","Tplus1_prop"), labels=c("initiation", "solution"))+ 
  scale_fill_manual(values=c(colours_repair_part[1],colours_repair_part[2]), guide=F)+
  scale_colour_manual(values=c(colours_repair_part[1],colours_repair_part[2]), labels=c("initiation", "solution"))+
  ggtitle("**Division of multimodal effort** (proportion of joint multimodal effort)")+ 
  labs(y="proportion of joint multimodal effort", x="")+
  scale_y_continuous(labels = scales::percent, limits=c(0,1))+
  expand_limits(x= c(0, 0))+
  geom_hline(yintercept=.5, linetype="dashed", colour="darkgrey", size=.4)+
  theme_light()+
  theme(
        legend.position="none",
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.background = element_rect(fill="#f6f6f6"),
        strip.background = element_rect(fill="grey", colour="grey"), #set colour of panel headings
        strip.placement = "outside",
        axis.line.y = element_line(colour="grey"),
        axis.title.y=element_text(size=7), axis.text.y=element_text(size=7),
        axis.text.x=element_blank(), axis.ticks.x=element_blank(),
        plot.margin = margin(5.5,0,0,5.5), #top, right, bottom, left+
        plot.title=ggtext::element_markdown(colour="black", size=10), #necessary to get only part in bold
        strip.text.x=ggtext::element_markdown(colour="black", size=8)) #necessary to get 'n' in italics


p_division_multi_labour_prop


```

### joint multimodal effort (Figure 2 in ms)
```{r joint multimodal effort [fig 2]}
p_joint_multi_labour <- ggplot(df_div_labour, aes(x=as.factor(1), y=multi_T0+multi_Tplus1, fill=T0_type)) +
  facet_manual(~T0_type, design="ABC", widths = unit(3.5, "cm"))+
  facet_grid(~T0_type, labeller=repair_type_labeller, switch="both")+
  gghalves::geom_half_boxplot(alpha=.6, outlier.shape=NA, width=1, lwd=0.5, fatten=1, errorbar.length = 0.2)+
  gghalves::geom_half_point_panel(position= position_nudge(x=.5), size=.05, side = "l", range_scale = 1.2, alpha=.5) +
  scale_fill_manual(values=c("lightgrey","lightgrey", "lightgrey"), guide="none")+ #colours of boxplots
  expand_limits(x= c(0, 0))+
  ggtitle("Joint multimodal effort")+ 
  labs(y="joint multimodal effort", x="")+
  theme_light()+
  theme(
        legend.title=element_blank(),
        legend.background=element_rect(fill="#f6f6f6"), legend.key=element_rect(fill="#f6f6f6"), 
        legend.position=c(0.83,0.852), legend.direction="horizontal",
        legend.key.height=unit(0.15,"cm"), #size of blue-red bar in legend
        legend.key.width=unit(0.28,"cm"), #size of blue-red bar in legend
        legend.text=element_text(size=7),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.background = element_rect(fill="#f6f6f6"),
        strip.background = element_rect(fill="grey", colour="grey"), #set colour of panel headings
        strip.placement = "outside",
        plot.title = element_text(size=10, face="bold"),
        axis.line.y = element_line(colour="grey"),
        axis.title.y=element_text(size=7), axis.text.y=element_text(size=7),
        axis.text.x=element_blank(), axis.ticks.x=element_blank(),
        plot.margin = margin(5.5,0,0,7), #top, right, bottom, left+
        strip.text.x=ggtext::element_markdown(colour="black", size=8)) #necessary to get 'n' in italics
  
p_joint_multi_labour


```

### density plots (Figure 4 in ms)
```{r density plots [fig 4]}

p_density_speech <- ggplot(data=df, aes(x=speech_count))+
  stat_halfeye(size=.2, fill="#a1bae6")+
  labs(y="density", x="orthographic characters per turn")+
  ylim(-0.07, 1)+
  theme_light()+
  theme(
    plot.title = element_text(hjust = 0.5, size=8), #centered title
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.border = element_blank(),
    axis.line.x = element_line(colour="grey"), axis.title.x=element_text(size=7), axis.text.x=element_text(size=6), #x-axis
    axis.line.y=element_line(colour="grey"), axis.title.y=element_text(size=7), axis.text.y=element_text(size=6))

p_density_gesture <- ggplot(data=df, aes(x=sum_subm))+
  stat_halfeye(size=.2, fill="#a1bae6")+
  labs(y="density", x="gesture submovements per turn")+
  ylim(-0.07, 1)+
  #xlim(0,20)+
  theme_light()+
  theme(
    plot.title = element_text(hjust = 0.5, size=8), #centered title
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.border = element_blank(),
    axis.line.x = element_line(colour="grey"), axis.title.x=element_text(size=7), axis.text.x=element_text(size=6), #x-axis
    axis.line.y=element_blank(), axis.title.y=element_blank(), axis.text.y=element_blank(), axis.ticks.y=element_blank()) #delete y-axis

p_densities <- plot_grid(p_density_speech, p_density_gesture,  nrow=1)
p_densities
```

```{r align and save figures}
p_div_effort_aligned <- plot_grid(p_speech_effort, NULL, p_gesture_effort, NULL, p_division_multi_labour_prop, 
                                         nrow=5, rel_heights=c(1, -0.04, 1, -0.04 ,1.17), align="v", labels=c("A","","B","","C"), label_size=10)
p_joint_effort_aligned <- plot_grid(p_joint_multi_labour, NULL, rel_widths=c(1, 0.3), nrow=1)

ggsave("3_figures/fig1_division_of_effort.pdf", p_div_effort_aligned, width=14, height=17, units="cm",dpi=900)
ggsave("3_figures/fig2_joint_effort.pdf", p_joint_effort_aligned, width=14, height=7, dpi=900, units="cm")

ggsave("3_figures/fig4_density_plots.png", p_densities, width=14, height=3.5, dpi=900, units="cm")
ggsave("3_figures/fig4_density_plots.eps", p_densities, width=14, height=3.5, dpi=900, units="cm")

```


# Statistics

```{r contrast coding}
#contrast coding
#backward difference coding, because we want to compare each level with the 'prior' level (so restr. request vs open request, and restr. offer vs. restr. request)
#https://stats.idre.ucla.edu/r/library/r-library-contrast-coding-systems-for-categorical-variables/#backward 


backward.diff = matrix(c(-2/3, 1/3, 1/3, -1/3, -1/3, 2/3), ncol = 2)



```

```{r stats frequency of repair types}
#set up df 
#with rows for all combinations of dyads (20) * targets (1:16) * T0_type (3 types) = 960 rows
pairnrs <- unique(df$pairnr)

df_empty <- expand.grid(pairnr = pairnrs, target = 1:16, T0_type = c("openrequest","restrictedrequest","restrictedoffer"))
df_freq <- df_div_labour %>% group_by(pairnr, target, T0_type) %>% count()

df_freq_per_type <- merge(df_empty, df_freq, by=c("pairnr","target","T0_type"), all.x=T)
df_freq_per_type[is.na(df_freq_per_type)] <- 0

#check if merging was OK; frequencies overall still the same?
df_freq_per_type %>% group_by(T0_type) %>% summarize(sum(n))



#assigning the backward difference coding to T0_type
contrasts(df_freq_per_type$T0_type) = backward.diff


freq_1 <- lmer(n ~ T0_type + (1+T0_type|pairnr) + (1+T0_type|target), data = df_freq_per_type, control=lmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e5)))
#singular fit

freq_2 <- lmer(n ~ T0_type + (1+T0_type|pairnr) + (1|target), data = df_freq_per_type, control=lmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e5)))
#singular fit

freq_3 <- lmer(n ~ T0_type + (1|pairnr) + (1+T0_type|target), data = df_freq_per_type, control=lmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e5)))
#singular fit

freq_4 <- lmer(n ~ T0_type + (1|pairnr) + (1|target), data = df_freq_per_type, control=lmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e5)))
summary(freq_4)


#descriptives
df_freq_per_type %>% group_by(T0_type) %>% summarise(n=sum(n)) %>% mutate(prop=n/sum(n))
```


```{r stats division of multimodal labour}
#assigning the backward difference coding to T0_type
contrasts(df_div_labour$T0_type) = backward.diff

multi_1 <- lmer(multi_T0_prop ~ T0_type + (1+T0_type|pairnr) + (1+T0_type|target), data = df_div_labour, control=lmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e5)))
#singular fit

multi_2 <- lmer(multi_T0_prop ~ T0_type + (1+T0_type|pairnr) + (1|target), data = df_div_labour, control=lmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e5)))
#singular fit

multi_3 <- lmer(multi_T0_prop ~ T0_type + (1|pairnr) + (1+T0_type|target), data = df_div_labour, control=lmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e5)))
#singular fit

multi_4 <- lmer(multi_T0_prop ~ T0_type + (1|pairnr) + (1|target), data = df_div_labour, control=lmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e5)))
#singular fit

multi_5  <- lmer(multi_T0_prop ~ T0_type + (1|pairnr), data = df_div_labour, control=lmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e5)))

summary(multi_5)

#negative correlation between effort in initiation and effort in solution
cor.test(df_div_labour$multi_T0, df_div_labour$multi_Tplus1)



```

Replicate the Dingemanse et al. finding, using nr. of orthographic characters
```{r stats division of speech labour}
speech_1 <- lmer(speech_T0_prop ~ T0_type + (1+T0_type|pairnr) + (1+T0_type|target), data = df_div_labour, control=lmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e5)))
#singular fit

speech_2 <- lmer(speech_T0_prop ~ T0_type + (1+T0_type|pairnr) + (1|target), data = df_div_labour, control=lmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e5)))
#singular fit

speech_3 <- lmer(multi_T0_prop ~ T0_type + (1|pairnr) + (1+T0_type|target), data = df_div_labour, control=lmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e5)))
#singular fit

speech_4 <- lmer(speech_T0_prop ~ T0_type + (1|pairnr) + (1|target), data = df_div_labour, control=lmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e5)))
#singular fit

speech_5  <- lmer(speech_T0_prop ~ T0_type + (1|pairnr), data = df_div_labour, control=lmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e5)))

summary(speech_5)

```

Least collaborative effort principle; most efficient repair type used most?
```{r stats least collaborative effort}
joint_effort_1 <- lmer(multi_T0+multi_Tplus1 ~ T0_type + (1+T0_type|pairnr) + (1+T0_type|target), data = df_div_labour, control=lmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e5)))
#singular fit

joint_effort_2 <- lmer(multi_T0+multi_Tplus1 ~ T0_type + (1+T0_type|pairnr) + (1|target), data = df_div_labour, control=lmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e5)))

joint_effort_3 <- lmer(multi_T0+multi_Tplus1 ~ T0_type + (1|pairnr) + (1+T0_type|target), data = df_div_labour, control=lmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e5)))
#singular fit

summary(joint_effort_2)
```


```{r session info, include=F}
sessionInfo()
```

R version 4.0.2 (2020-06-22)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19044)

Matrix products: default

locale:
[1] LC_COLLATE=Dutch_Netherlands.1252  LC_CTYPE=Dutch_Netherlands.1252    LC_MONETARY=Dutch_Netherlands.1252 LC_NUMERIC=C                      
[5] LC_TIME=Dutch_Netherlands.1252    

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] pastecs_1.3.21 ggh4x_0.2.1    ggdist_3.0.1   lmerTest_3.1-3 lme4_1.1-23    Matrix_1.2-18  cowplot_1.0.0  ggplot2_3.3.5  dplyr_1.0.0   
[10] tidyr_1.1.0   

loaded via a namespace (and not attached):
 [1] tinytex_0.23         statmod_1.4.34       inlmisc_0.5.2        tidyselect_1.1.0     xfun_0.14            purrr_0.3.4         
 [7] splines_4.0.2        lattice_0.20-41      colorspace_1.4-1     vctrs_0.3.8          generics_0.1.1       htmltools_0.4.0     
[13] yaml_2.2.1           utf8_1.1.4           rlang_0.4.11         nloptr_1.2.2.1       gridtext_0.1.1       pillar_1.6.4        
[19] glue_1.4.1           withr_2.2.0          sp_1.4-2             distributional_0.3.0 lifecycle_1.0.1      munsell_0.5.0       
[25] gtable_0.3.0         raster_3.4-13        codetools_0.2-16     evaluate_0.14        knitr_1.28           gghalves_0.1.1      
[31] fansi_0.4.1          Rcpp_1.0.7           scales_1.1.1         farver_2.0.3         digest_0.6.25        numDeriv_2016.8-1.1 
[37] grid_4.0.2           cli_2.0.2            rgdal_1.5-27         tools_4.0.2          magrittr_1.5         tibble_3.1.6        
[43] crayon_1.3.4         pkgconfig_2.0.3      ellipsis_0.3.2       MASS_7.3-51.6        data.table_1.12.8    xml2_1.3.2          
[49] assertthat_0.2.1     minqa_1.2.4          rmarkdown_2.2        rstudioapi_0.11      R6_2.4.1             boot_1.3-25         
[55] ggtext_0.1.0         igraph_1.2.6         nlme_3.1-148         compiler_4.0.2      

