Date rendered: Thu May 21 10:29:22 2020

Objectives

The goal of this project is to measure the effect of resource dispersion on individual social network metrics (strength, betweenness, and clustering coefficient). This script first performs a node-level permutation to create 2000 randomizations of the dataset created in the “BeetlearySocialNetworks” script. This script then generates a null distribution of coefficients from models that use these randomized datasets. After generating a null distribution of coefficients, this script then compares the coefficients from the model that uses the observed dataset and calculates p-values.

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. The dead beetles, however, are removed from the analyses in this script. Is this how I want to deal with beetle deaths, replacements, and movements? -When performing the node-level permutations, should I shuffle network metrics within or across periods? Currently, I shuffle metrics within a period. -When performing the node-level permutations, does it matter if I shuffle the social network metric vs. the treatment? -Should I include number of scans seen, period, sex, or elytra in models? -What model assumptions need to be addressed when testing significance by permuting data? I don’t love any of these models. Do I need to run better models? How?

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.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 tidyr::expand() masks Matrix::expand()
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
## x tidyr::pack()   masks Matrix::pack()
## x dplyr::recode() masks car::recode()
## x purrr::some()   masks car::some()
## x tidyr::unpack() masks Matrix::unpack()
library(foreach)
## 
## Attaching package: 'foreach'
## The following objects are masked from 'package:purrr':
## 
##     accumulate, when
library(doParallel)
## Loading required package: iterators
## Loading required package: parallel
#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 observed individual-level social network data
IndividSN<-read.csv("SN_20200518.csv")

#upload surveys data
observeddata<-read.csv("Beetleary2018_20200124.csv")

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

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

CWW Function to Calculate P-Values

# FUNCTION: CWW (2015) Calculates one- and two-tailed p-values from a permuted distribution  
# Author Corlett Wolfe Wood

# One-tailed: The proportion of permuted values that are GREATER THAN the observed value
# Two-tailed: The proportion of permuted values that are MORE EXTREME than the observed

# NB: Calculation of two-tailed p-values on asymmetric null distributions
# On symmetric distributions, one-tailed p-values can be doubled to obtained two-tailed ones
# This doesn't work with asymmetric distributions because of skewness
# --> SOLUTION: This function below calculates two-tailed p-values as the proportion of
#     permuted.values that have a probability <= the probability of the observed value
# --> SOURCE: Pratt, J. W. and J. D. Gibbons. 1981. Concepts in Nonparametric Theory. 
#     Springer-Verlag, New York, NY pp. 29-32

# NB: p-values should NEVER be zero in a permutation test
# --> the observed data is one possible permutation, so at least one value
#     in the permuted data can be equal to the observed data (P>0)
# --> SOURCE: North, B. V., D. Curtis, and P. C. Sham. 2002. A note on the calculation of empirical 
#     p values from Monte Carlo procedures. Am. J. Hum. Genet. 71:439-441.
# --> SOURCE: Davison, A. C. and D. V. Hinkley. 1997. Bootstrap methods and their application. 
#     Cambridge University Press, Cambridge, United Kingdom
# --> SOLUTION: add 1 to the numerator and the denominator

pvals = function(observed.value, permuted.values){  
  # ONE-TAILED
  # Proportion of permuted values >= (if in upper tail) OR <= (if in lower tail) the observed value
  k.greater = length(which(permuted.values>=observed.value))
  k.less = length(which(permuted.values<=observed.value))
  n = length(permuted.values)
  # One-tailed p-value
  pval.onetail = (min(k.greater, k.less)+1)/(n+1) # See NB above for explanation of +1  
  
  # TWO-TAILED
  # Proportion of permuted values with PROBABILITIES <= the observed value
  df <- data.frame(permuted.values=permuted.values, prob=NA)
  names(df) <- c("permuted.values", "prob")
  for(p in 1:nrow(df)){
    # Calculate probability of each permuted value
    # Proportion of permuted values >= (if in upper tail) OR <= (if in lower tail) each permuted value
    k.greater = length(which(df$permuted.values>=df$permuted.values[p]))
    k.less = length(which(df$permuted.values<=df$permuted.values[p]))
    n = nrow(df)
    df$prob[p] = (min(k.greater, k.less)+1)/(n+1) # See NB above for explanation of +1
  }
  
  # Proportion of permuted.values with a have a probability <= the probability of the observed value
  prob.less <- length(which(df$prob<=pval.onetail))
  pval.twotail <- (prob.less+1)/(n+1) # See NB above for explanation of +1
  
  return(c(one.tailed=pval.onetail, 
           two.tailed=pval.twotail))
}

Clean Up Individ SN Data

Remove Dead Beetles

Currently, dead beetels are included in the calculation of social networks but removed from the analyses. I want to remove dead beetles from the IndividSN dataset before performing permutations. This way, the distribution of social network metrics is the same for all datasets.

Similarly, beetles that moved (3MZ/3XY) are also included in the calculation of social network metrics. 3MZ is removed from the analysis. 3XY is removed from 1B 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(Focal %sans% DeadBeetles) %>%
  #remove 3MZ
  filter(Focal != "3MZ") %>%
  #remove 3XY from 1B
  filter(Focal_Condo_Period != "3XY_1B_1") %>%
  #droplevels
  droplevels() -> IndividSN

Add Columns

IndividSN %>%
  #create Period column
  mutate(Period = as.factor(substr(Focal_Condo_Period, 8, 8))) %>%
  #create Condo_Period column
  mutate(Condo_Period = as.factor(substr(Focal_Condo_Period, 5, 8))) %>%
  #create Condo column
  mutate(Condo = as.factor(substr(Focal_Condo_Period, 5, 6))) %>%
  #rename Focal to ID
  dplyr::rename(ID= Focal) -> 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

Calculate Number of Scans Seen

Rename X_fkfieldID_surveys

observeddata %>%
  dplyr::rename(ID = X_fkfieldID_surveys) -> observeddata

Divide Dataset into two Periods

observeddata %>%
  mutate(Period = ifelse(Date_Surveys < 20180719, "1", "2")) -> observeddata

Period1<-subset(observeddata, Period=="1")
Period2<-subset(observeddata, Period=="2")

Number of scans seen

Period1 %>%
  ungroup() %>%
  group_by(ID) %>%
  dplyr::mutate(ScansSeen = length(unique(ScanID)))-> Period1

Period2 %>%
  ungroup() %>%
  group_by(ID) %>%
  dplyr::mutate(ScansSeen = length(unique(ScanID))) -> Period2

Collapse Part 1 and 2

obssurveydata<-bind_rows(Period1, Period2)

Collapse data by beetle

#select only one row for each beetle_period
obssurveydata %>%
  ungroup() %>%
  #make a ID_Period column
  mutate(ID_Period = as.factor(paste(ID, Period, sep="_"))) %>%
  #select only one row for each beetle_period
  distinct(ID_Period, .keep_all = TRUE) -> obssurveydataBeetleset

Remove Extra Columns

obssurveydataBeetleset %>%
  dplyr::select(Survey_Sex, ScansSeen, ID_Period) -> obssurveydataBeetleset

Merge Number Scans Seen with SN Metrics

#create ID_Period in SN dataset to merge
IndividSN %>%
  mutate(ID_Period = as.factor(paste(ID, Period, sep="_"))) -> IndividSN

IndividSN <- merge(IndividSN, obssurveydataBeetleset, by="ID_Period", all.x=TRUE, all.y=FALSE)

Node-Level Permutations

Set Number of Permutations

permutations=100

SHUFFLE

Shuffle strength (alpha), betweenness, and clustering coefficient (am) across Focal_Condo but within Period. Shuffle without replacement.

#create empty lists of permuted data
PermIndividSN<-rep(list(list()), length(permutations))

#for each permutation
for(i in 1:permutations) {
  PermIndividSN[[i]] <- IndividSN %>% group_by(Period) %>% mutate(alpha=sample(alpha)) %>% mutate(betweenness=sample(betweenness)) %>% mutate(am=sample(am))
}

#an alternative way to permute is to permute the Treatment
#for(i in 1:permutations) {
#  PermIndividSN[[i]] <- IndividSN %>% group_by(Period) %>% mutate(Treatment=sample(Treatment))
#}

#name each permutation sequentially
names(PermIndividSN) <- paste("perm", 1:length(PermIndividSN), sep="_")

Observed Models

Strength

obs.strength <- glmmTMB(alpha~Treatment+ScansSeen+(1|Condo/ID), data=IndividSN)
summary(obs.strength)
##  Family: gaussian  ( identity )
## Formula:          alpha ~ Treatment + ScansSeen + (1 | Condo/ID)
## Data: IndividSN
## 
##      AIC      BIC   logLik deviance df.resid 
##   2029.6   2057.7  -1008.8   2017.6      794 
## 
## Random effects:
## 
## Conditional model:
##  Groups   Name        Variance Std.Dev.
##  ID:Condo (Intercept) 0.28093  0.5300  
##  Condo    (Intercept) 0.01326  0.1152  
##  Residual             0.49433  0.7031  
## Number of obs: 800, groups:  ID:Condo, 400; Condo, 12
## 
## Dispersion estimate for gaussian family (sigma^2): 0.494 
## 
## Conditional model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) 0.549284   0.158561   3.464 0.000532 ***
## Treatment1  0.070685   0.024969   2.831 0.004642 ** 
## ScansSeen   0.050317   0.003488  14.427  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(obs.strength, type=3)
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: alpha
##                Chisq Df Pr(>Chisq)    
## (Intercept)  12.0006  1  0.0005318 ***
## Treatment     8.0137  1  0.0046424 ** 
## ScansSeen   208.1473  1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(allEffects(obs.strength))

check.assumptions(obs.strength)

#test for zero inflation and overdispersion with DHARMa
obs.strength.resid<-simulateResiduals(obs.strength, n=250, integerResponse=F)
## It seems you are diagnosing a glmmTMB model. There are still a few minor limitations associated 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(obs.strength.resid)

testDispersion(obs.strength.resid)
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## ratioObsSim = 0.61685, p-value < 2.2e-16
## alternative hypothesis: two.sided
testZeroInflation(obs.strength.resid)
## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = NaN, p-value = 1
## alternative hypothesis: two.sided
#the model is underdispersed

#extract estimate for treatment
obs.strength.treatment.coef<-summary(obs.strength)$coefficients$cond[2,1]

Betweenness

How best deal with betweenness data? QQ Plots look awful.

obs.betweenness <- glmmTMB(betweenness~Treatment+ScansSeen+(1|Condo/ID), data=IndividSN)
summary(obs.betweenness)
##  Family: gaussian  ( identity )
## Formula:          betweenness ~ Treatment + ScansSeen + (1 | Condo/ID)
## Data: IndividSN
## 
##      AIC      BIC   logLik deviance df.resid 
##   7140.9   7169.0  -3564.4   7128.9      794 
## 
## Random effects:
## 
## Conditional model:
##  Groups   Name        Variance  Std.Dev. 
##  ID:Condo (Intercept) 2.621e+01  5.119222
##  Condo    (Intercept) 1.450e-06  0.001204
##  Residual             4.087e+02 20.215221
## Number of obs: 800, groups:  ID:Condo, 400; Condo, 12
## 
## Dispersion estimate for gaussian family (sigma^2):  409 
## 
## Conditional model:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -11.20175    3.80904  -2.941  0.00327 ** 
## Treatment1    0.02569    0.71710   0.036  0.97142    
## ScansSeen     0.70190    0.08644   8.120 4.64e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(obs.betweenness, type=3)
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: betweenness
##               Chisq Df Pr(>Chisq)    
## (Intercept)  8.6485  1   0.003273 ** 
## Treatment    0.0013  1   0.971420    
## ScansSeen   65.9423  1  4.643e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(allEffects(obs.betweenness))

check.assumptions(obs.betweenness)

#awful QQ plots!

#specify negative binomial family with glmmTMB
obs.betweenness.nbinom2<-glmmTMB(betweenness~Treatment+ScansSeen+(1|Condo/ID), data=IndividSN, family=nbinom2)
summary(obs.betweenness.nbinom2)
##  Family: nbinom2  ( log )
## Formula:          betweenness ~ Treatment + ScansSeen + (1 | Condo/ID)
## Data: IndividSN
## 
##      AIC      BIC   logLik deviance df.resid 
##   6166.7   6194.8  -3077.3   6154.7      794 
## 
## Random effects:
## 
## Conditional model:
##  Groups   Name        Variance  Std.Dev. 
##  ID:Condo (Intercept) 5.792e-09 7.611e-05
##  Condo    (Intercept) 1.955e-10 1.398e-05
## Number of obs: 800, groups:  ID:Condo, 400; Condo, 12
## 
## Overdispersion parameter for nbinom2 family (): 0.591 
## 
## Conditional model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) 0.874548   0.260218   3.361 0.000777 ***
## Treatment1  0.027849   0.047272   0.589 0.555786    
## ScansSeen   0.046528   0.005913   7.869 3.58e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(obs.betweenness.nbinom2, type=3)
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: betweenness
##               Chisq Df Pr(>Chisq)    
## (Intercept) 11.2952  1  0.0007771 ***
## Treatment    0.3471  1  0.5557860    
## ScansSeen   61.9210  1  3.575e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(allEffects(obs.betweenness.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(obs.betweenness.nbinom2)

#test for zero inflation and overdispersion with DHARMa
obs.betweenness.nbinom2.resid<-simulateResiduals(obs.betweenness.nbinom2, n=250)
## It seems you are diagnosing a glmmTMB model. There are still a few minor limitations associated 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(obs.betweenness.nbinom2.resid)

testDispersion(obs.betweenness.nbinom2.resid)
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## ratioObsSim = 0.77081, p-value < 2.2e-16
## alternative hypothesis: two.sided
testZeroInflation(obs.betweenness.nbinom2.resid)
## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 1.4147, p-value < 2.2e-16
## alternative hypothesis: two.sided
#underdispersed and zero inflated.QQ plot residuals still look awful

#specify zero inflation and negative binomial family
obs.betweenness.0inflation.nbinom2<-glmmTMB(betweenness~Treatment+ScansSeen+(1|Condo/ID), data=IndividSN, family=nbinom2, ziformula=~1)
## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')
#error: Model convergence problem; extreme or very small eigen values detected

#extract estimate for treatment
obs.betweenness.treatment.coef<-summary(obs.betweenness.nbinom2)$coefficients$cond[2,1]

Clustering Coefficient

QQ Plot is a bit weird

obs.cc <- glmmTMB(am~Treatment+ScansSeen+(1|Condo/ID), data=IndividSN, na.action=na.omit)
summary(obs.cc)
##  Family: gaussian  ( identity )
## Formula:          am ~ Treatment + ScansSeen + (1 | Condo/ID)
## Data: IndividSN
## 
##      AIC      BIC   logLik deviance df.resid 
##  -1418.5  -1390.5    715.2  -1430.5      784 
## 
## Random effects:
## 
## Conditional model:
##  Groups   Name        Variance  Std.Dev. 
##  ID:Condo (Intercept) 1.013e-11 3.184e-06
##  Condo    (Intercept) 3.382e-04 1.839e-02
##  Residual             9.400e-03 9.695e-02
## Number of obs: 790, groups:  ID:Condo, 399; Condo, 12
## 
## Dispersion estimate for gaussian family (sigma^2): 0.0094 
## 
## Conditional model:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) 0.4942339  0.0191132  25.858  < 2e-16 ***
## Treatment1  0.0143004  0.0034623   4.130 3.62e-05 ***
## ScansSeen   0.0008995  0.0004163   2.161   0.0307 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(obs.cc, type=3)
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: am
##                Chisq Df Pr(>Chisq)    
## (Intercept) 668.6519  1  < 2.2e-16 ***
## Treatment    17.0598  1  3.622e-05 ***
## ScansSeen     4.6683  1    0.03072 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(allEffects(obs.cc))

check.assumptions(obs.cc)

#test for zero inflation and overdispersion with DHARMa
obs.cc.resid<-simulateResiduals(obs.cc, n=250, integerResponse=FALSE)
## It seems you are diagnosing a glmmTMB model. There are still a few minor limitations associated 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.
plot(obs.cc.resid)

testDispersion(obs.cc.resid)
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## ratioObsSim = 0.96565, p-value = 0.264
## alternative hypothesis: two.sided
testZeroInflation(obs.cc.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
#zero inflated! And ugly QQ plot.
#error when fitting residuals: Model family was recognized or set as continuous, but duplicate values were detected in the response. Consider if you are fitting an appropriate model.
 
#specify 0inflation with glmmTMB
obs.cc.0inflation<-glmmTMB(am~Treatment+ScansSeen+(1|Condo/ID), data=IndividSN, na.action=na.omit, ziformula=~1)
summary(obs.cc.0inflation)
##  Family: gaussian  ( identity )
## Formula:          am ~ Treatment + ScansSeen + (1 | Condo/ID)
## Zero inflation:      ~1
## Data: IndividSN
## 
##      AIC      BIC   logLik deviance df.resid 
##  -1438.2  -1405.5    726.1  -1452.2      783 
## 
## Random effects:
## 
## Conditional model:
##  Groups   Name        Variance  Std.Dev. 
##  ID:Condo (Intercept) 1.061e-11 3.257e-06
##  Condo    (Intercept) 3.221e-04 1.795e-02
##  Residual             8.785e-03 9.373e-02
## Number of obs: 790, groups:  ID:Condo, 399; Condo, 12
## 
## Dispersion estimate for gaussian family (sigma^2): 0.00878 
## 
## Conditional model:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) 0.5071460  0.0185891  27.282  < 2e-16 ***
## Treatment1  0.0128283  0.0033530   3.826  0.00013 ***
## ScansSeen   0.0006302  0.0004045   1.558  0.11921    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -5.9773     0.7087  -8.435   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(obs.cc.0inflation, type=3)
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: am
##                Chisq Df Pr(>Chisq)    
## (Intercept) 744.3042  1  < 2.2e-16 ***
## Treatment    14.6377  1  0.0001303 ***
## ScansSeen     2.4277  1  0.1192056    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(allEffects(obs.cc.0inflation))
check.assumptions(obs.cc.0inflation)

#test for zero inflation and overdispersion with DHARMa
obs.cc.0inflation.resid<-simulateResiduals(obs.cc.0inflation, n=250, integerResponse=FALSE)
## It seems you are diagnosing a glmmTMB model. There are still a few minor limitations associated 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(obs.cc.0inflation.resid)

testDispersion(obs.cc.0inflation.resid)
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## ratioObsSim = 0.9665, p-value = 0.32
## alternative hypothesis: two.sided
testZeroInflation(obs.cc.0inflation.resid)
## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 1.0121, p-value = 1
## alternative hypothesis: two.sided
#not underdispersed nor zero inflated! QQ plot is a bit weird
#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? It is changing the integerResponse to TRUE. But this is weird given that am is not an integer.

#extract estimate for treatment
obs.cc.treatment.coef<-summary(obs.cc.0inflation)$coefficients$cond[2,1]

Permuted Models

Strength Model

#write a function that runs the model
perm.strength = function(PermIndividSN){
  perm.strength.mod<-glmmTMB(alpha~Treatment+ScansSeen+(1|Condo/ID), data=PermIndividSN)
  return(perm.strength.mod)
}

#run the model on each permutation
perm.strength.results<-lapply(PermIndividSN, FUN=perm.strength)
## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')
## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')
## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')
## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')
## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')
## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')
## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')
## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')
## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')
## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')
## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')
## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')
## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')
## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')
## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')
## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')
## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')
## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')
## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')
## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')
## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')
## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')
#error: Model convergence problem; non-positive-definite Hessian matrix.
#error: Model convergence problem; singular convergence.
#error: Model convergence problem; extreme of very small eigen values detected

Strength Model Estimates

#write a function that extracts the type 3 sum of squares (written by VAF)
AnovaT3=function(x){
  setTimeLimit(cpu = Inf,  transient = FALSE)
  return(tryCatch(Anova(x, type=3), error=function(e) NULL))
}

#write a function that extracts the model summary (written by VAF)
Model.Summary=function(x){
  setTimeLimit(cpu = Inf,  transient = FALSE)
  return(tryCatch(summary(x), error=function(e) NULL))
}

#extract type 3 sum of squares for each permutation
perm.strength.Anova<-lapply(perm.strength.results, FUN=AnovaT3)

#extract model summary for each permutation
perm.strength.summary<-lapply(perm.strength.results, FUN=Model.Summary)

#extract model estimate for treatment from model summary
perm.strength.treatment.coef<-foreach(i = 1:length(perm.strength.summary)) %dopar% {
  perm.strength.summary[[i]][["coefficients"]][["cond"]][[2]]
}
## Warning: executing %dopar% sequentially: no parallel backend registered
#flatten the list of permutations into a dataframe with one column that contains all the estimates from every permutation
perm.strength.treatment.coef.df<-do.call(rbind, perm.strength.treatment.coef)

Strength P-Values

#use CWW pval function to calculate pvalues
strength.pvals<-pvals(obs.strength.treatment.coef, perm.strength.treatment.coef.df)

strength.pvals
## one.tailed two.tailed 
## 0.02970297 0.04950495

Visualize Permuted Strength Estimates

#clean up dataframe to plot
perm.strength.treatment.coef.df<-as.data.frame(perm.strength.treatment.coef.df)
perm.strength.treatment.coef.df %>%
  dplyr::rename(PermutedEstimate = V1) -> perm.strength.treatment.coef.df

#set ggplot theme
theme_set(theme_classic(base_size=24)) 

#ggplot
ggplot(perm.strength.treatment.coef.df, aes(x=PermutedEstimate)) +
  geom_histogram() +
  geom_vline(xintercept = obs.strength.treatment.coef, color = "red")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Betweenness Model

#write a function that runs the model
perm.betweenness = function(PermIndividSN){
  perm.betweenness.mod<-glmmTMB(betweenness~Treatment+ScansSeen+(1|Condo/ID), data=PermIndividSN, family=nbinom2)
  return(perm.betweenness.mod)
}

#run the model on each permutation
perm.betweenness.results<-lapply(PermIndividSN, FUN=perm.betweenness)
## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')
## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')
## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')
## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')
## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')
## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')
## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')
## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')
## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')
## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')
## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')
## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')
## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')
## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')
## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')
## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')
## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')
## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')
## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')
## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')
## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')
## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')
#error: Model convergence problem; non-positive-definite Hessian matrix
#error: Model convergence problem; extreme or very small eigen values detected

Betweenness Model Estimates

#write a function that extracts the type 3 sum of squares (written by VAF)
AnovaT3=function(x){
  setTimeLimit(cpu = Inf,  transient = FALSE)
  return(tryCatch(Anova(x, type=3), error=function(e) NULL))
}

#write a function that extracts the model summary (written by VAF)
Model.Summary=function(x){
  setTimeLimit(cpu = Inf,  transient = FALSE)
  return(tryCatch(summary(x), error=function(e) NULL))
}

#extract type 3 sum of squares for each permutation
perm.betweenness.Anova<-lapply(perm.betweenness.results, FUN=AnovaT3)

#extract model summary for each permutation
perm.betweenness.summary<-lapply(perm.betweenness.results, FUN=Model.Summary)

#extract model estimate for treatment from model summary
perm.betweenness.treatment.coef<-foreach(i = 1:length(perm.betweenness.summary)) %dopar% {
  perm.betweenness.summary[[i]][["coefficients"]][["cond"]][[2]]
}

#flatten the list of permutations into a dataframe with one column that contains all the estimates from every permutation
perm.betweenness.treatment.coef.df<-do.call(rbind, perm.betweenness.treatment.coef)

Betweenness P-Values

#use CWW pval function to calculate pvalues
betweenness.pvals<-pvals(obs.betweenness.treatment.coef, perm.betweenness.treatment.coef.df)

betweenness.pvals
## one.tailed two.tailed 
##  0.2970297  0.5841584

Visualize Permuted Betweenness Estimates

#clean up dataframe to plot
perm.betweenness.treatment.coef.df<-as.data.frame(perm.betweenness.treatment.coef.df)
perm.betweenness.treatment.coef.df %>%
  dplyr::rename(PermutedEstimate = V1) -> perm.betweenness.treatment.coef.df

#set ggplot theme
theme_set(theme_classic(base_size=24)) 

#ggplot
ggplot(perm.betweenness.treatment.coef.df, aes(x=PermutedEstimate)) +
  geom_histogram() +
  geom_vline(xintercept = obs.betweenness.treatment.coef, color = "red")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

CC Model

#write a function that runs the model
perm.cc = function(PermIndividSN){
  perm.cc.mod<-glmmTMB(am~Treatment+ScansSeen+(1|Condo/ID), data=PermIndividSN, na.action=na.omit, ziformula=~1)
  return(perm.cc.mod)
}

#run the model on each permutation
perm.cc.results<-lapply(PermIndividSN, FUN=perm.cc)
## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')
## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')
## Warning in fitTMB(TMBStruc): Model convergence problem; singular convergence
## (7). See vignette('troubleshooting')
## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')
## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')
## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')
## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')
## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')
## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')
## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')
## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')
## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')
## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')
## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')
## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')
## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')
## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')
## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')
## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')
## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')
## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')
## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')
## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')
## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')
## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')
## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')
## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')
## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')
## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')
## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')
#error: Model convergence problem; extreme or very small eigen values detected 
#error: Model convergence problem; non-positive-definite Hessian matrix
#error: Model convergence problem; singular convergence
#Na/NaN function evaluation

CC Model Estimates

#write a function that extracts the type 3 sum of squares (written by VAF)
AnovaT3=function(x){
  setTimeLimit(cpu = Inf,  transient = FALSE)
  return(tryCatch(Anova(x, type=3), error=function(e) NULL))
}

#write a function that extracts the model summary (written by VAF)
Model.Summary=function(x){
  setTimeLimit(cpu = Inf,  transient = FALSE)
  return(tryCatch(summary(x), error=function(e) NULL))
}

#extract type 3 sum of squares for each permutation
perm.cc.Anova<-lapply(perm.cc.results, FUN=AnovaT3)

#extract model summary for each permutation
perm.cc.summary<-lapply(perm.cc.results, FUN=Model.Summary)

#extract model estimate for treatment from model summary
perm.cc.treatment.coef<-foreach(i = 1:length(perm.cc.summary)) %dopar% {
  perm.cc.summary[[i]][["coefficients"]][["cond"]][[2]]
}

#flatten the list of permutations into a dataframe with one column that contains all the estimates from every permutation
perm.cc.treatment.coef.df<-do.call(rbind, perm.cc.treatment.coef)

CC P-Values

#use CWW pval function to calculate pvalues
cc.pvals<-pvals(obs.cc.treatment.coef, perm.cc.treatment.coef.df)

cc.pvals
## one.tailed two.tailed 
## 0.00990099 0.00990099

Visualize Permuted CC Estimates

#clean up dataframe to plot
perm.cc.treatment.coef.df<-as.data.frame(perm.cc.treatment.coef.df)
perm.cc.treatment.coef.df %>%
  dplyr::rename(PermutedEstimate = V1) -> perm.cc.treatment.coef.df

#set ggplot theme
theme_set(theme_classic(base_size=24)) 

#ggplot
ggplot(perm.cc.treatment.coef.df, aes(x=PermutedEstimate)) +
  geom_histogram() +
  geom_vline(xintercept = obs.cc.treatment.coef, color = "red")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.