Date rendered: Mon May 18 17:10:52 2020

Objectives

The goal of this project is to measure the effect of resource dispersion on population social network metrics (average strength, density, global clustering coefficient, elytra assortment, sex assortment). This script uses paired t-tests to answer how resource dispersion affects population-level social network metrics. The population-level social network metrics are built from networks that include all interactions and both sexes in the “BeetlearySocialNetworks” script.

Outstanding questions

-What type of network do I want to analyze? Male only? Female only? Include mating interactions? Proximity only? Currently, this analysis is on social networks created with both sexes and ALL interactions. -The social network metrics are calcualted from social networks that include beetles that died but do not include beetles that were replaced early in the experiment or beetles that moved into network 1B (3MZ and 3XY) during the experiment. Is this how I want to deal with beetle deaths, replacements, and movements? Should beetles that die be included in the global social network metrics?

Set-Up and Data

#clear memory
rm(list=ls())

#libraries
library(tidyverse)
## ── Attaching packages ───────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.2.1     ✓ purrr   0.3.3
## ✓ tibble  2.1.3     ✓ dplyr   0.8.3
## ✓ tidyr   1.0.2     ✓ stringr 1.4.0
## ✓ readr   1.3.1     ✓ forcats 0.4.0
## ── Conflicts ──────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
#function that excludes
`%sans%` <- Negate(`%in%`)

#Change R's defult to contrasts 
options(contrasts=c("contr.sum", "contr.poly")) 

#upload observed individual-level social network data
ObsPopSN<-read.csv("PopSN_20200518.csv")

#Treatments
Dispersed <- c("1A_1", "2A_1", "3A_1", "4A_1", "5A_1", "6A_1", "1B_2", "2B_2", "3B_2", "4B_2", "5B_2", "6B_2")
Clumped <- c("1B_1", "2B_1", "3B_1", "4B_1", "5B_1", "6B_1", "1A_2", "2A_2", "3A_2", "4A_2", "5A_2", "6A_2")

Clean Up Observed SN Dataset

Separate out columns to create Condo, Period

ObsPopSN %>%
  #create Condo column
  mutate(Condo = as.factor(substr(Condo_Period, 1, 2))) %>%
  #create Period column
  mutate(Period = as.factor(substr(Condo_Period, 4, 4))) -> ObsPopSN

Treatment

Assign each Condo in each Period to a Dispersed or Clumped Treatment

ObsPopSN %>%
  mutate(Treatment = ifelse(Condo_Period %in% Dispersed, "Dispersed", "Clumped")) %>%
  mutate(Treatment = as.factor(Treatment)) -> ObsPopSN

Restructure Population SN Metric Dataframe for Paired T-Tests

#first divide two dataframes into Clumped and Dispersed
ObsPopSN %>%
  filter(Treatment == "Clumped") %>%
  dplyr::select(avg.strength, global.cc, observed.elytra.assortment, observed.sex.assortment, density, Condo) %>%
  dplyr::rename(avg.strength.clump=avg.strength, global.cc.clump=global.cc, elytra.assortment.clump=observed.elytra.assortment, sex.assortment.clump=observed.sex.assortment, density.clump=density) -> PopSNClump

ObsPopSN %>%
  filter(Treatment == "Dispersed") %>%
  dplyr::select(avg.strength, global.cc, observed.elytra.assortment, observed.sex.assortment, density, Condo) %>%
  dplyr::rename(avg.strength.disperse=avg.strength, global.cc.disperse=global.cc, elytra.assortment.disperse=observed.elytra.assortment, sex.assortment.disperse=observed.sex.assortment, density.disperse=density) -> PopSNDisperse

PopSNLong <- merge(PopSNClump, PopSNDisperse, by="Condo", all.x=TRUE, all.y=FALSE)

Do Population Social Network Metrics Differ Across Treatments?

Density

#paired t-test
t.test(PopSNLong$density.clump, PopSNLong$density.disperse, paired=TRUE, conf.level=0.95)
## 
##  Paired t-test
## 
## data:  PopSNLong$density.clump and PopSNLong$density.disperse
## t = -0.10519, df = 11, p-value = 0.9181
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.05415292  0.04921283
## sample estimates:
## mean of the differences 
##            -0.002470045
#check assumptions
density.difference = PopSNLong$density.clump - PopSNLong$density.disperse
hist(density.difference)

#visualize result
ggplot(ObsPopSN, aes(x=Treatment, y=density)) +
  geom_boxplot()

Global Clustering Coefficient

#paired t-test
t.test(PopSNLong$global.cc.clump, PopSNLong$global.cc.disperse, paired=TRUE, conf.level=0.95)
## 
##  Paired t-test
## 
## data:  PopSNLong$global.cc.clump and PopSNLong$global.cc.disperse
## t = 1.1656, df = 11, p-value = 0.2684
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.02323353  0.07554584
## sample estimates:
## mean of the differences 
##              0.02615615
#check assumptions
global.cc.difference = PopSNLong$global.cc.clump - PopSNLong$global.cc.disperse
hist(global.cc.difference)

#visualize result
ggplot(ObsPopSN, aes(x=Treatment, y=global.cc)) +
  geom_boxplot()

Power Analysis

Why does individual clustering coefficient but not global clustering coefficient vary across resource dispersion treatments? power: the probability of rejecting null hypothesis when the alternative hypothesis is true, inverse of Type II Error

#power test
#calculate difference between global.cc.clump and global.cc.dispers
PopSNLong$global.cc.diff<-(PopSNLong$global.cc.clump-PopSNLong$global.cc.disperse)
#mean difference
mdiff<-mean(PopSNLong$global.cc.diff)
#sd difference 
sddiff<-sd(PopSNLong$global.cc.diff)

power.t.test(n=12, delta=mdiff, sd=sddiff, sig.level=0.05, power=NULL, type="paired", alternative="two.sided")
## 
##      Paired t test power calculation 
## 
##               n = 12
##           delta = 0.02615615
##              sd = 0.07773378
##       sig.level = 0.05
##           power = 0.1856761
##     alternative = two.sided
## 
## NOTE: n is number of *pairs*, sd is std.dev. of *differences* within pairs

Elytra Assortment

#paired t-test
t.test(PopSNLong$elytra.assortment.clump, PopSNLong$elytra.assortment.disperse, paired=TRUE, conf.level=0.95)
## 
##  Paired t-test
## 
## data:  PopSNLong$elytra.assortment.clump and PopSNLong$elytra.assortment.disperse
## t = -0.0371, df = 11, p-value = 0.9711
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.04746165  0.04588814
## sample estimates:
## mean of the differences 
##            -0.000786754
#check assumptions
elytra.assortment.difference = PopSNLong$elytra.assortment.clump - PopSNLong$elytra.assortment.disperse
hist(elytra.assortment.difference)

#visualize result
ggplot(ObsPopSN, aes(x=Treatment, y=observed.elytra.assortment)) +
  geom_boxplot()

Sex Assortment

#paired t-test
t.test(PopSNLong$sex.assortment.clump, PopSNLong$sex.assortment.disperse, paired=TRUE, conf.level=0.95)
## 
##  Paired t-test
## 
## data:  PopSNLong$sex.assortment.clump and PopSNLong$sex.assortment.disperse
## t = 1.2696, df = 11, p-value = 0.2304
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.01568820  0.05845679
## sample estimates:
## mean of the differences 
##               0.0213843
#check assumptions
sex.assortment.difference = PopSNLong$sex.assortment.clump - PopSNLong$sex.assortment.disperse
hist(sex.assortment.difference)

#bonferoni corrected p value is 0.05/5=0.01. 
#visualize result
ggplot(ObsPopSN, aes(x=Treatment, y=observed.sex.assortment)) +
  geom_boxplot()