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

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)
#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"))
#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

#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]
##       min       max      mean   std.dev 
## 14.190850 34.555117 24.377813  5.500061
rm(df_ToT)

repair

#overall
freq_T0_bypair <- df %>% subset(repair_part=="T0") %>% group_by(pairnr, .drop=FALSE) %>% count() 
stat.desc(freq_T0_bypair$n)[stat_desc_output]
##       min       max      mean   std.dev 
##  6.000000 45.000000 18.900000  9.920208
#types
df %>% subset(repair_part=="T0") %>% group_by(T0_type) %>% count()
## # A tibble: 3 x 2
## # Groups:   T0_type [3]
##   T0_type               n
##   <fct>             <int>
## 1 openrequest          24
## 2 restrictedrequest    39
## 3 restrictedoffer     315
# 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]
##       min       max      mean   std.dev 
## 0.6374085 2.7936181 1.5332414 0.6089943
rm(repair_by_time, ToT_bypair)

speech (orthographic length)

tapply(df$speech_count, list(df$T0_type, df$repair_part), mean)
##                       Tmin1       T0    Tplus1
## openrequest       156.35000 20.12500 116.62500
## restrictedrequest  65.56250 27.05128  63.07692
## restrictedoffer    64.55642 34.37778  16.69524
tapply(df$speech_count, list(df$T0_type, df$repair_part), sd)
##                      Tmin1       T0   Tplus1
## openrequest       90.59062 12.21301 87.17838
## restrictedrequest 40.69076 16.63216 50.74518
## restrictedoffer   38.13173 20.18324 27.28575
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)
## 
##  Pearson's product-moment correlation
## 
## data:  df_characters_vs_duration$speech_count and df_characters_vs_duration$duration
## t = 70.197, df = 754, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.9211489 0.9401580
## sample estimates:
##       cor 
## 0.9312846

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))
##   sum_total sum_iconic sum_deictic sum_other
## 1       479        430          15        34
#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)
##   prop_iconic prop_deictic prop_other
## 1    89.77035     3.131524   7.098121
#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)
##     n
## 2 281
  #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)
## # A tibble: 12 x 5
## # Groups:   T0_type, repair_part, "total_count != 0" [6]
##    `"total_count != 0"` repair_part T0_type               n prop_gesture_in_turn
##    <chr>                <fct>       <fct>             <int>                <dbl>
##  1 total_count != 0     T0          openrequest          23               0.958 
##  2 total_count != 0     T0          openrequest           1               0.0417
##  3 total_count != 0     T0          restrictedrequest    34               0.872 
##  4 total_count != 0     T0          restrictedrequest     5               0.128 
##  5 total_count != 0     T0          restrictedoffer     161               0.511 
##  6 total_count != 0     T0          restrictedoffer     154               0.489 
##  7 total_count != 0     Tplus1      openrequest           5               0.208 
##  8 total_count != 0     Tplus1      openrequest          19               0.792 
##  9 total_count != 0     Tplus1      restrictedrequest    13               0.333 
## 10 total_count != 0     Tplus1      restrictedrequest    26               0.667 
## 11 total_count != 0     Tplus1      restrictedoffer     239               0.759 
## 12 total_count != 0     Tplus1      restrictedoffer      76               0.241
  #per sequence
  df_div_labour %>% subset(gesture_subm_T0!=0|gesture_subm_Tplus1!=0) %>% nrow()
## [1] 238

gesture submovements

tapply(df$sum_subm, list(df$T0_type,df$repair_part), mean)
##                      Tmin1         T0    Tplus1
## openrequest       6.600000 0.04166667 3.6250000
## restrictedrequest 2.000000 0.15384615 2.6666667
## restrictedoffer   2.155642 1.01587302 0.5269841
tapply(df$sum_subm, list(df$T0_type, df$repair_part), sd)
##                      Tmin1        T0   Tplus1
## openrequest       8.500155 0.2041241 3.842695
## restrictedrequest 2.651841 0.4315493 3.148377
## restrictedoffer   2.657371 1.5668954 1.741696

multimodal effort

tapply( (df_div_labour$multi_T0 + df_div_labour$multi_Tplus1), df_div_labour$T0_type, mean)
##       openrequest restrictedrequest   restrictedoffer 
##          3.514744          2.454998          1.372160
tapply( (df_div_labour$multi_T0 + df_div_labour$multi_Tplus1), df_div_labour$T0_type, sd)
##       openrequest restrictedrequest   restrictedoffer 
##          2.582550          1.741437          1.209902

Figures

#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)

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

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)

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)

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

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

#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)
#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))
## # A tibble: 3 x 2
##   T0_type           `sum(n)`
##   <fct>                <dbl>
## 1 openrequest             24
## 2 restrictedrequest       39
## 3 restrictedoffer        315
#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)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: n ~ T0_type + (1 | pairnr) + (1 | target)
##    Data: df_freq_per_type
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e+05))
## 
## REML criterion at convergence: 2158.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.8737 -0.3676 -0.0830  0.1453  6.4975 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  pairnr   (Intercept) 0.03179  0.1783  
##  target   (Intercept) 0.01840  0.1357  
##  Residual             0.52436  0.7241  
## Number of obs: 960, groups:  pairnr, 20; target, 16
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)   0.39375    0.05732  24.97946   6.869 3.38e-07 ***
## T0_type1      0.04688    0.05725 923.00000   0.819    0.413    
## T0_type2      0.86250    0.05725 923.00000  15.066  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##          (Intr) T0_ty1
## T0_type1  0.000       
## T0_type2  0.000 -0.500
#descriptives
df_freq_per_type %>% group_by(T0_type) %>% summarise(n=sum(n)) %>% mutate(prop=n/sum(n))
## # A tibble: 3 x 3
##   T0_type               n   prop
##   <fct>             <dbl>  <dbl>
## 1 openrequest          24 0.0635
## 2 restrictedrequest    39 0.103 
## 3 restrictedoffer     315 0.833
#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)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: multi_T0_prop ~ T0_type + (1 | pairnr)
##    Data: df_div_labour
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e+05))
## 
## REML criterion at convergence: 5.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9357 -0.5828  0.2918  0.7073  3.0639 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  pairnr   (Intercept) 0.004016 0.06337 
##  Residual             0.055182 0.23491 
## Number of obs: 378, groups:  pairnr, 20
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)   0.37824    0.02551  55.89161  14.828   <2e-16 ***
## T0_type1      0.14482    0.06173 368.01536   2.346   0.0195 *  
## T0_type2      0.46669    0.04086 374.49178  11.422   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##          (Intr) T0_ty1
## T0_type1 -0.181       
## T0_type2 -0.410 -0.583
#negative correlation between effort in initiation and effort in solution
cor.test(df_div_labour$multi_T0, df_div_labour$multi_Tplus1)
## 
##  Pearson's product-moment correlation
## 
## data:  df_div_labour$multi_T0 and df_div_labour$multi_Tplus1
## t = -2.7257, df = 376, p-value = 0.006717
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.23674107 -0.03887504
## sample estimates:
##        cor 
## -0.1391971

Replicate the Dingemanse et al. finding, using nr. of orthographic characters

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)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: speech_T0_prop ~ T0_type + (1 | pairnr)
##    Data: df_div_labour
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e+05))
## 
## REML criterion at convergence: -58.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9803 -0.5962  0.2004  0.7999  3.0384 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  pairnr   (Intercept) 0.002943 0.05425 
##  Residual             0.046653 0.21599 
## Number of obs: 378, groups:  pairnr, 20
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercep