Date rendered: Mon Sep 30 08:49:29 2019

To Do

-What does the DHARMA error mean??? -I’d like some feedback on how to deal with betweenness data. The QQ plot is ugly. -Why is the model singular for clustering coefficient? And is this treatment effect real? Seems like means are the same. -I can also run these analyses for male-male networks and female-female networks. -What is the story that I want to tell with this data? Simply that resource distribution does not affect social networks?

Objectives

The goal of this project is to compare individual- and population-level social network metrics between the two resource dispersion treatments.

NOTE (1) This data is not fully error-checked. (2) The social network metrics are calcualted from social networks that include beetles that died and moved (3MZ) during the experiment. But the analyses do not include beetles that died or 3MZ. (3) The beetle size measurements are elytra measurements done by Ruth and Patrick (work study students). In the future, we will have horn measurements as well and can calculate PCA. (4) These analyses ignore the lack of indepdence in network metrics. In the future, I need to incorporate permutations.

Set-Up and Data

#clear memory
rm(list=ls())

#libraries
library(lme4)
## Loading required package: Matrix
library(car)
## Loading required package: carData
## Registered S3 methods overwritten by 'car':
##   method                          from
##   influence.merMod                lme4
##   cooks.distance.influence.merMod lme4
##   dfbeta.influence.merMod         lme4
##   dfbetas.influence.merMod        lme4
library(glmmTMB)
library(DHARMa)
## Registered S3 method overwritten by 'DHARMa':
##   method        from   
##   refit.glmmTMB glmmTMB
library(effects)
## lattice theme set by effectsTheme()
## See ?effectsTheme for details.
library(tidyverse)
## ── Attaching packages ──────────────────────────────────────────────────────────────────────────── tidyverse 1.2.1 ──
## ✔ ggplot2 3.2.0     ✔ purrr   0.3.2
## ✔ tibble  2.1.3     ✔ dplyr   0.8.3
## ✔ tidyr   0.8.3     ✔ stringr 1.4.0
## ✔ readr   1.3.1     ✔ forcats 0.4.0
## ── Conflicts ─────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ tidyr::expand() masks Matrix::expand()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ✖ dplyr::recode() masks car::recode()
## ✖ purrr::some()   masks car::some()
#function that excludes
`%sans%` <- Negate(`%in%`)

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

#FUNCTION: check.assumptions
check.assumptions <- function(mod){
  # Multiple panels
  par(mfrow=c(2,2))
  
  # Normality
  hist(resid(mod))
  qqnorm(resid(mod))
  qqline(resid(mod))
  
  # Linearity
  # --> Plot residual versus fitted values with loess line
  plot(fitted(mod), resid(mod), main = "Residual vs fitted values")
  fr.fit <- loess(resid(mod)~fitted(mod))
  fr.vec <- seq(min(fitted(mod)), max(fitted(mod)), length.out=100)
  lines(fr.vec,predict(fr.fit,fr.vec), col="red")
  
  # Homoscedasticity
  # --> Plot the absolute value of the residuals against the fitted values
  plot(fitted(mod), sqrt(abs(resid(mod))), main="Scale-location")
  sl.fit <- loess(sqrt(abs(resid(mod)))~fitted(mod))
  sl.vec <- seq(min(fitted(mod)), max(fitted(mod)), length.out=100)
  lines(sl.vec,predict(sl.fit,sl.vec), col="red")
}

#upload individual-level social network data
IndividSN<-read.csv("socialnetworks_20190513_.csv")

#upload population-level social network data
PopSN<-read.csv("PopSN_20190513_.csv")

#upload surveys dataset
surveys_2018<-read.csv("20190422_OffMountainBeetleBase2.0.csv")
#NOTE: This is not a fully cleaned dataset.

#upload beetle measurements
elytra<-read.csv("BeetleMeasurements_20190324.csv")
#NOTE: This elytra data is not completely cleaned.
#(1) Need to check with Hannah Donald about how to measure elytra 3FM and 2MX. These are dirty scans. Ruth and Patrick measured these scans very differently. I already emailed Hannah about this (20190221) and just need to follow-up by looking at these measurements. Right now, these are wrong measurements!
#(2) Missing dorsal scan for 3NH. For now, I am using Rachel's measurement of the elytra. I do not have a pronotum measurement for 3NH. Luckily, this beetle died during Period 1. Rachel's measurement: 7.281

#create objects for ease of coding
DeadBeetles <- c("2BO", "2BX", "2H2", "2M6", "2MS", "2N6", "2N7", "2OH", "2OW", "2PK", "2PW", "2T9", "2TS", "2VB", "2VO", "2WV", "3A6", "3AB", "3AW", "3HK", "3HX", "3KA", "3LV", "3LX", "3MN", "3NH", "3ON", "3VO", "3W9", "3XM", "3XN", "2FF", "2HF", "3F3", "3KL", "3YF")

Dispersed <- c("1A_Period1", "2A_Period1", "3A_Period1", "4A_Period1", "5A_Period1", "6A_Period1", "1B_Period2", "2B_Period2", "3B_Period2", "4B_Period2", "5B_Period2", "6B_Period2")
Clumped <- c("1B_Period1", "2B_Period1", "3B_Period1", "4B_Period1", "5B_Period1", "6B_Period1", "1A_Period2", "2A_Period2", "3A_Period2", "4A_Period2", "5A_Period2", "6A_Period2")

Clean Up Individual Social Network Dataset

Separate out columns to create Condo, Period, ID columns

Social network metrics were not calculated for UK beetles.

IndividSN %>%
  #create Condo column
  mutate(Condo = as.factor(substr(Condo_Period, 1, 2))) %>%
  #create Period column
  mutate(Period = as.factor(substr(Condo_Period, 4, 10))) %>%
  #create ID column
  mutate(ID = as.factor(substr(ids, 1, 3))) -> IndividSN

Treatment

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

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

Remove beetles that died or moved (3MV/3XY) from the analyses

These beetles are currently included in the social networks, just removed from the analyses

3MZ seen in 1B starting on 20180702_A. It is supposed to be in 6A. We know 3MZ moved to 1B. We collected 3MZ from 1B on 20180706. 3XY moved from 1A to 1B. 3XY was seen in both 1A and 1B on 20180715_C. 3XY was only in 1B for 4 scans. Moved 3XY back to 1A.

IndividSN %>%
  #remove dead beetles
  filter(ID %sans% DeadBeetles) %>%
  #remove 3MZ
  filter(ID != "3MZ") %>%
  #remove 3XY from 1B
  filter(ids != "3XY_1B") -> IndividSN

Sex

Assign a Sex to each ID

IndividSN %>%
  #assign sex based on first letter. 2 = M, 3 = F
  mutate(Sex = ifelse(substr(ID, 1, 1) == "2", "M", "F")) %>%
  #make Sex a factor
  mutate(Sex = as.factor(Sex)) -> IndividSN

Elytra

Merge elytra dataset In the future, can use a PCA for body size. For now, use elytra measurements done by Ruth and Patrick (work-study students).

#remove pronotum and horn measurements
elytra %>%
  filter(Item.Measured == "Elytra") -> elytra

#add Rachel's measurement for 3NH. 3NH is missing a dorsal scan.
MissingDorsal <- c("3NH-D.jpg")
elytra %>%
  mutate(Average = replace(Average, Image.Name %in% MissingDorsal, 7.281)) -> elytra

#create simplified elytra dataset: just ID and Elytra Length
elytra %>%
  #create ID column
  mutate(ID = as.factor(substr(Image.Name, 1, 3))) %>%
  #rename elytra column
  dplyr::rename(Elytra = Average) %>%
  #simplify the dataset to just ID and Elytra Length
  dplyr::select(ID, Elytra) -> elytra

#merge body size measures to raw dataset
IndividSN <- merge(IndividSN, elytra, by=("ID"), all.x=TRUE, all.y=FALSE)

Merge Individual and Population SN metrics

SN <- merge(IndividSN, PopSN, by="Condo_Period", all.x=TRUE, all.y=FALSE)

Explore Individual Social Network Metrics

NOTE: How should I deal with betweenness data? It is very left-skewed.

#what does the data look like?
#degree
hist(SN$degree)

#output
hist(SN$output)

#strength
hist(SN$alpha)

#betweenness
hist(SN$betweenness)

#very left-skewed distribution. what do I do about this? taking the log produces -infinity values

#cc
hist(SN$am)

Correlation Among SN Metrics

#subset SN for SN Metrics
SN %>%
  dplyr::select(degree, output, alpha, betweenness, am) -> SNMetrics

#calculate correlation among SN metrics
cor(SNMetrics)
##                degree    output     alpha betweenness am
## degree      1.0000000 0.9005541 0.9697288   0.5420912 NA
## output      0.9005541 1.0000000 0.9791138   0.7510443 NA
## alpha       0.9697288 0.9791138 1.0000000   0.6711591 NA
## betweenness 0.5420912 0.7510443 0.6711591   1.0000000 NA
## am                 NA        NA        NA          NA  1
#clustering coefficient (am) is 0 for some individuals. get correlations with am and other SN metrics
cor.test(SN$degree, SN$am, method=c("pearson"))
## 
##  Pearson's product-moment correlation
## 
## data:  SN$degree and SN$am
## t = 1.228, df = 781, p-value = 0.2198
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.02624406  0.11361319
## sample estimates:
##        cor 
## 0.04389964
cor.test(SN$output, SN$am, method=c("pearson"))
## 
##  Pearson's product-moment correlation
## 
## data:  SN$output and SN$am
## t = 2.6084, df = 781, p-value = 0.00927
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.02302015 0.16194180
## sample estimates:
##        cor 
## 0.09293325
cor.test(SN$alpha, SN$am, method=c("pearson"))
## 
##  Pearson's product-moment correlation
## 
## data:  SN$alpha and SN$am
## t = 2.0302, df = 781, p-value = 0.04268
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.002403876 0.141797826
## sample estimates:
##        cor 
## 0.07245466
cor.test(SN$betweenness, SN$am, method=c("pearson"))
## 
##  Pearson's product-moment correlation
## 
## data:  SN$betweenness and SN$am
## t = -2.1019, df = 781, p-value = 0.03588
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.144304192 -0.004962568
## sample estimates:
##         cor 
## -0.07499948

Clean Up Population Social Network Dataset

Note: These Pop SN Metrics are calculated from SRI matrices that include all males and females.

PopSN %>%
  mutate(Condo = substr(Condo_Period, 1, 2)) %>%
  mutate(Period = substr(Condo_Period, 4, 10)) %>%
  mutate(Treatment = ifelse(Condo_Period %in% Dispersed, "Dispersed", "Clumped")) %>%
  mutate(Condo = as.factor(Condo)) %>%
  mutate(Period = as.factor(Period)) %>%
  mutate(Treatment = as.factor(Treatment)) -> PopSN

Correlation Among Pop SN Metrics

NOTE: Given these correlations, I am choosing to analyze density, elytra assortment, sex assortment, global clustering coefficient. Density, average strength, average shortest path, and CV all capture network connectedness.

#subset PopSN for Population Metrics
PopSN %>%
  ungroup() %>%
  dplyr::select(avg.strength, avg.shortest.path, global.cc, elytra.assortment, sex.assortment, density, CV) -> PopSNMetrics

#calculate correlation among Pop SN metrics
cor(PopSNMetrics)
##                   avg.strength avg.shortest.path  global.cc
## avg.strength         1.0000000        -0.9011753  0.7727890
## avg.shortest.path   -0.9011753         1.0000000 -0.6289503
## global.cc            0.7727890        -0.6289503  1.0000000
## elytra.assortment    0.2795325        -0.3308289  0.3327614
## sex.assortment       0.4851818        -0.4259942  0.5037023
## density              0.9267853        -0.9041021  0.7441608
## CV                  -0.8813793         0.8689923 -0.6833799
##                   elytra.assortment sex.assortment    density         CV
## avg.strength             0.27953253     0.48518175  0.9267853 -0.8813793
## avg.shortest.path       -0.33082887    -0.42599422 -0.9041021  0.8689923
## global.cc                0.33276142     0.50370235  0.7441608 -0.6833799
## elytra.assortment        1.00000000     0.02419616  0.2098711 -0.2082098
## sex.assortment           0.02419616     1.00000000  0.4415824 -0.6244768
## density                  0.20987112     0.44158238  1.0000000 -0.9339359
## CV                      -0.20820980    -0.62447681 -0.9339359  1.0000000

Restructure Population SN Metric Dataframe for Paired T-Tests

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

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

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

Do Individual Social Network Metrics Differ Between Treatments?

Strength

alpha.comparison <- lmer(alpha~Treatment*Sex+Treatment*Elytra+(1|Condo/ID), data=SN)
summary(alpha.comparison)
## Linear mixed model fit by REML ['lmerMod']
## Formula: alpha ~ Treatment * Sex + Treatment * Elytra + (1 | Condo/ID)
##    Data: SN
## 
## REML criterion at convergence: 1908.2
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -3.10298 -0.61044  0.01178  0.61747  2.55517 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID:Condo (Intercept) 0.07815  0.2795  
##  Condo    (Intercept) 0.01102  0.1050  
##  Residual             0.53887  0.7341  
## Number of obs: 800, groups:  ID:Condo, 400; Condo, 12
## 
## Fixed effects:
##                   Estimate Std. Error t value
## (Intercept)       0.998592   0.472266   2.114
## Treatment1        0.010494   0.414632   0.025
## Sex1              0.411785   0.029895  13.774
## Elytra            0.145958   0.067745   2.155
## Treatment1:Sex1   0.037944   0.026305   1.442
## Treatment1:Elytra 0.003134   0.059602   0.053
## 
## Correlation of Fixed Effects:
##             (Intr) Trtmn1 Sex1   Elytra Tr1:S1
## Treatment1   0.000                            
## Sex1        -0.163  0.000                     
## Elytra      -0.996  0.000  0.163              
## Trtmnt1:Sx1  0.000 -0.163  0.000  0.000       
## Trtmnt1:Ely  0.000 -0.998  0.000  0.000  0.163
Anova(alpha.comparison, type=3)
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: alpha
##                     Chisq Df Pr(>Chisq)    
## (Intercept)        4.4710  1    0.03448 *  
## Treatment          0.0006  1    0.97981    
## Sex              189.7278  1    < 2e-16 ***
## Elytra             4.6419  1    0.03120 *  
## Treatment:Sex      2.0806  1    0.14918    
## Treatment:Elytra   0.0028  1    0.95806    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(allEffects(alpha.comparison))

check.assumptions(alpha.comparison)

#test for zero inflation and overdispersion with DHARMa
alpha.comparison.resid<-simulateResiduals(alpha.comparison, n=250)
## Model family was recognized or set as continuous, but duplicate values were detected in the response. Consider if you are fitting an appropriate model.
plot(alpha.comparison.resid)

testDispersion(alpha.comparison.resid)
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted
##  vs. simulated
## 
## data:  simulationOutput
## ratioObsSim = 0.99562, p-value = 0.816
## alternative hypothesis: two.sided
testZeroInflation(alpha.comparison.resid)
## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = Inf, p-value < 2.2e-16
## alternative hypothesis: two.sided
#the model is VERY zero-inflated 
#BUT I am shocked that this is zero-inflated. What is going on here??
#when I create residuals, I get the message: Model family was recognized or set as continuous, but duplicate values were detected in the response. Consider if you are fitting an appropriate model. WHAT DOES THIS ERROR MESSAGE MEAN?

Betweenness

NOTE: how deal with betweenness data?

betweenness.comparison <- lmer(betweenness~Treatment*Sex+Treatment*Elytra+(1|Condo/ID), data=SN)
## boundary (singular) fit: see ?isSingular
summary(betweenness.comparison)
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## betweenness ~ Treatment * Sex + Treatment * Elytra + (1 | Condo/ID)
##    Data: SN
## 
## REML criterion at convergence: 7456.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.2782 -0.4926 -0.2886  0.2624  5.7618 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID:Condo (Intercept)   6.571   2.563  
##  Condo    (Intercept)   0.000   0.000  
##  Residual             663.067  25.750  
## Number of obs: 800, groups:  ID:Condo, 400; Condo, 12
## 
## Fixed effects:
##                   Estimate Std. Error t value
## (Intercept)         4.4820    14.6880   0.305
## Treatment1         18.3070    14.5445   1.259
## Sex1               10.0903     0.9318  10.828
## Elytra              2.2243     2.1113   1.054
## Treatment1:Sex1    -0.1435     0.9227  -0.155
## Treatment1:Elytra  -2.7030     2.0907  -1.293
## 
## Correlation of Fixed Effects:
##             (Intr) Trtmn1 Sex1   Elytra Tr1:S1
## Treatment1   0.000                            
## Sex1        -0.163  0.000                     
## Elytra      -0.998  0.000  0.163              
## Trtmnt1:Sx1  0.000 -0.163  0.000  0.000       
## Trtmnt1:Ely  0.000 -0.998  0.000  0.000  0.163
## convergence code: 0
## boundary (singular) fit: see ?isSingular
Anova(betweenness.comparison, type=3)
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: betweenness
##                     Chisq Df Pr(>Chisq)    
## (Intercept)        0.0931  1     0.7603    
## Treatment          1.5843  1     0.2081    
## Sex              117.2539  1     <2e-16 ***
## Elytra             1.1099  1     0.2921    
## Treatment:Sex      0.0242  1     0.8764    
## Treatment:Elytra   1.6714  1     0.1961    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(allEffects(betweenness.comparison))

check.assumptions(betweenness.comparison)

#error: boundary (singular) fit. betweenness is extremely 0-inflated and non-normal. 

#specify 0-inflation and negative binomial family with glmmTMB
betweenness.comparison.0inflation.nbinom2<-glmmTMB(betweenness~Treatment*Sex+Treatment*Elytra+(1|Condo/ID), data=SN, family=nbinom2, ziformula=~1)
summary(betweenness.comparison.0inflation.nbinom2)
##  Family: nbinom2  ( log )
## Formula:          
## betweenness ~ Treatment * Sex + Treatment * Elytra + (1 | Condo/ID)
## Zero inflation:               ~1
## Data: SN
## 
##      AIC      BIC   logLik deviance df.resid 
##   5933.5   5980.3  -2956.7   5913.5      790 
## 
## Random effects:
## 
## Conditional model:
##  Groups   Name        Variance  Std.Dev. 
##  ID:Condo (Intercept) 1.861e-02 1.364e-01
##  Condo    (Intercept) 2.140e-09 4.626e-05
## Number of obs: 800, groups:  ID:Condo, 400; Condo, 12
## 
## Overdispersion parameter for nbinom2 family (): 0.71 
## 
## Conditional model:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        0.944990   0.735307   1.285  0.19874    
## Treatment1         0.447261   0.714506   0.626  0.53133    
## Sex1               0.521367   0.058824   8.863  < 2e-16 ***
## Elytra             0.294170   0.105308   2.793  0.00522 ** 
## Treatment1:Sex1    0.003603   0.052265   0.069  0.94504    
## Treatment1:Elytra -0.068374   0.102257  -0.669  0.50372    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.7547     0.2005   -8.75   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(betweenness.comparison.0inflation.nbinom2, type=3)
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: betweenness
##                    Chisq Df Pr(>Chisq)    
## (Intercept)       1.6516  1   0.198735    
## Treatment         0.3918  1   0.531333    
## Sex              78.5562  1  < 2.2e-16 ***
## Elytra            7.8032  1   0.005216 ** 
## Treatment:Sex     0.0048  1   0.945038    
## Treatment:Elytra  0.4471  1   0.503719    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(allEffects(betweenness.comparison.0inflation.nbinom2))
## Warning in Effect.glmmTMB(predictors, mod, vcov. = vcov., ...): overriding
## variance function for effects: computed variances may be incorrect

## Warning in Effect.glmmTMB(predictors, mod, vcov. = vcov., ...): overriding
## variance function for effects: computed variances may be incorrect

check.assumptions(betweenness.comparison.0inflation.nbinom2)

#test for zero inflation and overdispersion with DHARMa
betweenness.comparison.0inflation.nbinom2.resid<-simulateResiduals(betweenness.comparison.0inflation.nbinom2, n=250)
## It seems you are diagnosing a glmmTBM model. There are still a few minor limitations associatd with this package. The most important is that glmmTMB doesn't implement an option to create unconditional predictions from the model, which means that predicted values (in res ~ pred) plots include the random effects. With strong random effects, this can sometimes create diagonal patterns from bottom left to top right in the res ~ pred plot. All other tests and plots should work as desired. Please see https://github.com/florianhartig/DHARMa/issues/16 for further details.
plot(betweenness.comparison.0inflation.nbinom2.resid)

testDispersion(betweenness.comparison.0inflation.nbinom2.resid)
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted
##  vs. simulated
## 
## data:  simulationOutput
## ratioObsSim = 0.87723, p-value = 0.064
## alternative hypothesis: two.sided
testZeroInflation(betweenness.comparison.0inflation.nbinom2.resid)
## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 1.0699, p-value = 0.384
## alternative hypothesis: two.sided
#no longer over/underdispersed nor zero inflated. BUT QQ plot residuals look awful

Clustering Coefficient

Remove individuals that have no clustering coefficient metrics

NoCC <- c("2A3", "2BS", "2FH", "2FM", "2FN", "2FP", "2FZ", "2NN", "2NT", "2SW", "2T6", "2TX", "3F2", "3PN", "3PW", "3TP")

SN %>%
  #remove beetles that do not have clustering coefficient metrics
  filter(ID %sans% NoCC) %>%
  droplevels() -> SNCC
am.comparison <- lmer(am~Treatment*Sex+Treatment*Elytra+(1|Condo/ID), data=SNCC)
## boundary (singular) fit: see ?isSingular
summary(am.comparison)
## Linear mixed model fit by REML ['lmerMod']
## Formula: am ~ Treatment * Sex + Treatment * Elytra + (1 | Condo/ID)
##    Data: SNCC
## 
## REML criterion at convergence: -1005.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.9424 -0.5208  0.0252  0.5444  4.3417 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. 
##  ID:Condo (Intercept) 1.187e-10 1.089e-05
##  Condo    (Intercept) 6.028e-04 2.455e-02
##  Residual             1.464e-02 1.210e-01
## Number of obs: 768, groups:  ID:Condo, 384; Condo, 12
## 
## Fixed effects:
##                    Estimate Std. Error t value
## (Intercept)        0.449161   0.070787   6.345
## Treatment1        -0.160244   0.070275  -2.280
## Sex1              -0.008486   0.004442  -1.910
## Elytra             0.003757   0.010121   0.371
## Treatment1:Sex1   -0.002560   0.004438  -0.577
## Treatment1:Elytra  0.023548   0.010099   2.332
## 
## Correlation of Fixed Effects:
##             (Intr) Trtmn1 Sex1   Elytra Tr1:S1
## Treatment1   0.000                            
## Sex1        -0.179  0.000                     
## Elytra      -0.993  0.000  0.179              
## Trtmnt1:Sx1  0.000 -0.179  0.000  0.000       
## Trtmnt1:Ely  0.000 -0.998  0.000  0.000  0.178
## convergence code: 0
## boundary (singular) fit: see ?isSingular
Anova(am.comparison, type=3)
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: am
##                    Chisq Df Pr(>Chisq)    
## (Intercept)      40.2622  1  2.221e-10 ***
## Treatment         5.1996  1    0.02259 *  
## Sex               3.6494  1    0.05609 .  
## Elytra            0.1378  1    0.71048    
## Treatment:Sex     0.3327  1    0.56407    
## Treatment:Elytra  5.4370  1    0.01971 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(allEffects(am.comparison))

check.assumptions(am.comparison)

#error: boundary (singular) fit. WHY? I removed individuals that did not have a clustering coefficient metric. 

#specify 0inflation with glmmTMB
am.comparison.0inflation<-glmmTMB(am~Treatment*Sex+Treatment*Elytra+(1|Condo/ID), data=SNCC, ziformula=~1)
summary(am.comparison.0inflation)
##  Family: gaussian  ( identity )
## Formula:          
## am ~ Treatment * Sex + Treatment * Elytra + (1 | Condo/ID)
## Zero inflation:      ~1
## Data: SNCC
## 
##      AIC      BIC   logLik deviance df.resid 
##  -1051.6  -1005.2    535.8  -1071.6      758 
## 
## Random effects:
## 
## Conditional model:
##  Groups   Name        Variance  Std.Dev. 
##  ID:Condo (Intercept) 3.903e-11 6.247e-06
##  Condo    (Intercept) 4.123e-04 2.030e-02
##  Residual             1.282e-02 1.132e-01
## Number of obs: 768, groups:  ID:Condo, 384; Condo, 12
## 
## Dispersion estimate for gaussian family (sigma^2): 0.0128 
## 
## Conditional model:
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        0.4731649  0.0664755   7.118  1.1e-12 ***
## Treatment1        -0.1748917  0.0659791  -2.651  0.00803 ** 
## Sex1              -0.0105778  0.0042004  -2.518  0.01179 *  
## Elytra             0.0008703  0.0095114   0.091  0.92710    
## Treatment1:Sex1    0.0016338  0.0041992   0.389  0.69722    
## Treatment1:Elytra  0.0253833  0.0094783   2.678  0.00741 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -4.7645     0.4108   -11.6   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(am.comparison.0inflation, type=3)
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: am
##                    Chisq Df Pr(>Chisq)    
## (Intercept)      50.6643  1  1.096e-12 ***
## Treatment         7.0263  1   0.008032 ** 
## Sex               6.3418  1   0.011793 *  
## Elytra            0.0084  1   0.927096    
## Treatment:Sex     0.1514  1   0.697221    
## Treatment:Elytra  7.1719  1   0.007405 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(allEffects(am.comparison.0inflation))

plot(effect("Treatment", am.comparison.0inflation))
## NOTE: Treatment is not a high-order term in the model

check.assumptions(am.comparison.0inflation)

#I don't understand how this treatment effect is significant? Calculate means to compare.
#Clumped: 0.4779163
#Dispersed: 0.4714959
SNCC %>%
  filter(Treatment == "Clumped") %>%
  summarise(mean = mean(am))
##        mean
## 1 0.4779163
SNCC %>%
  filter(Treatment == "Dispersed") %>%
  summarise(mean = mean(am))
##        mean
## 1 0.4714959
#test for zero inflation and overdispersion with DHARMa
am.comparison.0inflation.resid<-simulateResiduals(am.comparison.0inflation, n=250)
## It seems you are diagnosing a glmmTBM model. There are still a few minor limitations associatd with this package. The most important is that glmmTMB doesn't implement an option to create unconditional predictions from the model, which means that predicted values (in res ~ pred) plots include the random effects. With strong random effects, this can sometimes create diagonal patterns from bottom left to top right in the res ~ pred plot. All other tests and plots should work as desired. Please see https://github.com/florianhartig/DHARMa/issues/16 for further details.
## Model family was recognized or set as continuous, but duplicate values were detected in the response. Consider if you are fitting an appropriate model.
## Model family was recognized or set as continuous, but duplicate values were detected in the simulation - changing to integer residuals (see ?simulateResiduals for details)
plot(am.comparison.0inflation.resid)

testDispersion(am.comparison.0inflation.resid)
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted
##  vs. simulated
## 
## data:  simulationOutput
## ratioObsSim = 0.97349, p-value = 0.432
## alternative hypothesis: two.sided
testZeroInflation(am.comparison.0inflation.resid)
## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 1.0651, p-value = 0.976
## alternative hypothesis: two.sided
#not underdispersed nor zero inflated!
#when I create residuals, I get the message: Model family was recognized or set as continuous, but duplicate values were detected in the response. Consider if you are fitting an appropriate model. - changing to integer residuals. WHAT IS THIS?

Do Population Social Network Metrics Differ Across Treatments?

Avg. Strength

#paired t-test
t.test(PopSNLong$avg.strength.clump, PopSNLong$avg.strength.disperse, paired=TRUE, conf.level=0.95)
## 
##  Paired t-test
## 
## data:  PopSNLong$avg.strength.clump and PopSNLong$avg.strength.disperse
## t = 0.09286, df = 11, p-value = 0.9277
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.2586448  0.2814306
## sample estimates:
## mean of the differences 
##              0.01139292
#check assumptions
avg.strength.difference = PopSNLong$avg.strength.clump - PopSNLong$avg.strength.disperse
hist(avg.strength.difference)

#visualize result
ggplot(PopSN, aes(x=Treatment, y=avg.strength)) +
  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 = 0.93082, df = 11, p-value = 0.3719
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.02380960  0.05870655
## sample estimates:
## mean of the differences 
##              0.01744847
#check assumptions
global.cc.difference = PopSNLong$global.cc.clump - PopSNLong$global.cc.disperse
hist(global.cc.difference)

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

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.5647, df = 11, p-value = 0.5836
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.07880996  0.04662704
## sample estimates:
## mean of the differences 
##             -0.01609146
#check assumptions
elytra.assortment.difference = PopSNLong$elytra.assortment.clump - PopSNLong$elytra.assortment.disperse
hist(elytra.assortment.difference)

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

Sex Assortment

This is the only resource distribution effect. Sexes are more negatively assorted when resources are dispersed. This effect is not signficant when account for bonferonni correction.

#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 = 2.3494, df = 11, p-value = 0.03853
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.002011522 0.061684715
## sample estimates:
## mean of the differences 
##              0.03184812
#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(PopSN, aes(x=Treatment, y=sex.assortment)) +
  geom_boxplot()

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.1208, df = 11, p-value = 0.906
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.06035546  0.05407477
## sample estimates:
## mean of the differences 
##            -0.003140346
#check assumptions
density.difference = PopSNLong$density.clump - PopSNLong$density.disperse
hist(density.difference)

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