# path to project directory
knitr::opts_knit$set(root.dir = normalizePath(".."))

# fig quality
knitr::opts_chunk$set(dpi = 50, fig.width = 12)

Exploratory graphs

# read data
dat <- read_excel("./data/Datos de respuestas_cap2.xlsx")

# convert to regular data frame
dat <- as.data.frame(dat)

# create new variable abou
dat$est_repr[dat$estado_repr == 0] <- "Inactivo"
dat$est_repr[dat$estado_repr == 1] <- "Activo"

# aggregate total number of calls
agg_dat <- aggregate(n_llamadas ~ ID + sexo_consulta + sexo_respuesta + est_repr, data = dat, FUN = sum)

# plot
ggplot(agg_dat, aes(fill = sexo_consulta, y = n_llamadas, x = sexo_respuesta)) +
geom_boxplot() +
  scale_fill_viridis_d(alpha = 0.7, begin = 0.4) +
theme_classic(base_size = 24) +
  labs(x = "Sexo respuesta", y = "Numero de llamadas")+
facet_wrap(~ est_repr) + ggtitle("Total de llamadas")

# create binary variable for calling
dat$n_llam_bin <- ifelse(dat$n_llamadas > 0, 1, 0)

# aggregate
agg_dat2 <- aggregate(n_llam_bin ~ ID + sexo_consulta + sexo_respuesta + est_repr, data = dat, FUN = sum)

#plot 
ggplot(agg_dat2, aes(fill = sexo_consulta, y = n_llam_bin, x = sexo_respuesta)) +
  geom_boxplot() +
  scale_fill_viridis_d(alpha = 0.7, begin = 0.4) +
  theme_classic(base_size = 24) +
  labs(x = "Sexo respuesta", y = "Cuenta de experimentos con respuesta")+
facet_wrap(~ est_repr) + ggtitle("Experimentos con respuesta")

Multinomial mixed effect models (MCMCglmm)

# define parmeters for MCMCglmm models
itrns <- 300000
burnin <- 3000
thin <- 1000

# prior for effect models
pr <- list(B = list(mu = rep(0, 8), V = diag(8) * (1 + pi^2/3)), R = list(V = 1, fix = 1),  G = list(G1 = list(V = 1, fix = 1)))

# run model
md <- MCMCglmm(n_llam_bin ~ sexo_respuesta:est_repr:sexo_consulta  - 1, random = ~ ID, data = dat, family = "categorical", prior = pr, verbose = FALSE, nitt = itrns, start = list(QUASI = FALSE), burnin = burnin, thin = thin)

Caculating p-values:

The output contains the posterior distribution of the parameter estimates. These parameter distributions can be used to test specific hypothesis about differences between sexes/stages/inquiry sexes. Column names in md$Sol (solutions) refer to the combination of levels from the 3 interacting variables:

# simplify names
colnames(md$Sol) <- gsub("sexo_respuesta|est_repr", "", colnames(md$Sol))

colnames(md$Sol)
## [1] "H:Activo:sexo_consultaH"   "M:Activo:sexo_consultaH"  
## [3] "H:Inactivo:sexo_consultaH" "M:Inactivo:sexo_consultaH"
## [5] "H:Activo:sexo_consultaM"   "M:Activo:sexo_consultaM"  
## [7] "H:Inactivo:sexo_consultaM" "M:Inactivo:sexo_consultaM"
# stack posteriors
Y <- stack(as.data.frame(md$Sol))

# plot posteriors
ggplot(Y, aes(x=values)) + 
  geom_vline(xintercept = 0, col = "red", lty = 2) +
  geom_density() + 
  labs(y = "Density", x = "Posterior") +
  facet_wrap(~ ind, ncol = 2) +
  theme_classic(base_size = 24)

Estimate of the overlap of posteriors can be used a statistical significance test. For instance we can compare calling activity between sexes during the active stage:

# get p value
p.val <- sum(md$Sol[, "M:Activo:sexo_consultaH"] - md$Sol[, "H:Activo:sexo_consultaH"] < 0) / nrow(md$Sol)

p.val
## [1] 0.01010101
# plot distributions
ggplot(Y[Y$ind %in% c("H:Activo:sexo_consultaH", "M:Activo:sexo_consultaH"), ], aes(x=values)) + geom_density(aes(group = ind, colour = ind, fill = ind), alpha=0.3) + 
      scale_color_viridis_d(alpha = 0.7, begin = 0.4, end = 0.8) +
    scale_fill_viridis_d(alpha = 0.7, begin = 0.4, end = 0.8) +
  labs(y = "Density", x = "Posterior") +
  theme_classic(base_size = 24)

Or whether males responde more to females during the reproductive stage:

# p value
p.val <- sum(md$Sol[, "M:Activo:sexo_consultaH"] - md$Sol[, "M:Inactivo:sexo_consultaH"] < 0) / nrow(md$Sol)

p.val
## [1] 0.04040404
# plot distributions
ggplot(Y[Y$ind %in% c("M:Activo:sexo_consultaM", "M:Inactivo:sexo_consultaM"), ], aes(x=values)) + geom_density(aes(group = ind, colour = ind, fill = ind), alpha=0.3) + 
    scale_color_viridis_d(alpha = 0.7, begin = 0.4, end = 0.8) +
    scale_fill_viridis_d(alpha = 0.7, begin = 0.4, end = 0.8) +
    labs(y = "Density", x = "Posterior") +
  theme_classic(base_size = 24)

Diagnostic plots

plot(md$Sol)

R session information

## R version 3.6.1 (2019-07-05)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.4 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/atlas/libblas.so.3.10.3
## LAPACK: /usr/lib/x86_64-linux-gnu/atlas/liblapack.so.3.10.3
## 
## locale:
##  [1] LC_CTYPE=es_ES.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=es_CR.UTF-8        LC_COLLATE=es_ES.UTF-8    
##  [5] LC_MONETARY=es_CR.UTF-8    LC_MESSAGES=es_ES.UTF-8   
##  [7] LC_PAPER=es_CR.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=es_CR.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] MCMCglmm_2.29     ape_5.4           coda_0.19-3       Matrix_1.2-18    
## [5] viridis_0.5.1     viridisLite_0.3.0 readxl_1.3.1      ggplot2_3.3.1    
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.4.6     pillar_1.4.4     compiler_3.6.1   cellranger_1.1.0
##  [5] tools_3.6.1      digest_0.6.25    nlme_3.1-142     lattice_0.20-41 
##  [9] evaluate_0.14    lifecycle_0.2.0  tibble_3.0.1     gtable_0.3.0    
## [13] pkgconfig_2.0.3  rlang_0.4.6      parallel_3.6.1   yaml_2.2.1      
## [17] xfun_0.14        gridExtra_2.3    withr_2.2.0      stringr_1.4.0   
## [21] dplyr_1.0.0      knitr_1.28       generics_0.0.2   vctrs_0.3.1     
## [25] grid_3.6.1       tidyselect_1.1.0 glue_1.4.1       R6_2.4.1        
## [29] rmarkdown_2.2    tensorA_0.36.1   farver_2.0.3     corpcor_1.6.9   
## [33] purrr_0.3.4      magrittr_1.5     scales_1.1.1     ellipsis_0.3.1  
## [37] htmltools_0.4.0  cubature_2.0.4   colorspace_1.4-1 labeling_0.3    
## [41] stringi_1.4.6    munsell_0.5.0    crayon_1.3.4