Preliminary analysis

library(tidyverse)
library(brms)
library(modelr)
library(broom)
library(huxtable)
library(scales)
library(here)

CHAINS <- 4
ITER <- 4000
WARMUP <- 2000
BAYES_SEED <- 1234
options(mc.cores = parallel::detectCores())

constraints <- read_rds(here("data", "derived_data", "constraints.rds")) 

constraint_levels <- constraints %>% 
  mutate(constraint_clean = fct_inorder(constraint_clean)) %>% 
  mutate(levels = map(levels_clean, names)) %>% 
  unnest(levels) %>% 
  mutate(levels = fct_inorder(levels))

# TODO: Switch back to full BLS.rds and BHL.rds
BHL <- read_rds(here("data", "derived_data", "BHL_small.rds"))
bhl_long <- BHL %>% 
  gather(constraint, value, one_of(constraints$constraint)) %>% 
  left_join(constraints, by = "constraint") %>% 
  mutate(value = fct_relevel(value, levels(constraint_levels$levels)),
         constraint_clean = fct_relevel(constraint_clean, 
                                        levels(constraint_levels$constraint_clean)))

BLS <- read_rds(here("data", "derived_data", "BLS_small.rds"))
bls_long <- BLS %>% 
  gather(constraint, value, one_of(constraints$constraint)) %>% 
  left_join(constraints, by = "constraint") %>% 
  mutate(value = fct_relevel(value, levels(constraint_levels$levels)),
         constraint_clean = fct_relevel(constraint_clean, 
                                        levels(constraint_levels$constraint_clean)))

Single conditions, no nested interactions

Landscape fitness

High-low

ggplot(bhl_long, aes(x = value, y = avg_evenness_t, fill = repp)) +
  geom_violin(size = 0.25, position = position_dodge(width = 0.75)) +
  stat_summary(geom = "point", fun = "mean", 
               aes(group = repp), position = position_dodge(width = 0.75),
               size = 2, pch = 21, fill = "black", color = "white") +
  # scale_y_reverse(breaks = seq(0, 0.4, 0.2)) +
  # coord_cartesian(ylim = c(0, 0.5)) +
  scale_fill_viridis_d(option = "plasma", begin = 0.1, end = 0.9) +
  labs(x = NULL, y = "Landscape fitness, linked", fill = NULL, 
       caption = "Point = mean value") +
  facet_wrap(~ constraint_clean, scales = "free_x") +
  theme_bw() +
  theme(legend.position = "top",
        legend.key.size = unit(0.65, "lines"))

Latin squares

ggplot(bls_long, 
       aes(x = value, y = landscape_fitness_linked, fill = constraint_clean)) +
  geom_violin(size = 0.25) +
  stat_summary(geom = "point", fun = "mean", 
               size = 2, pch = 21, fill = "black", color = "white") +
  scale_y_reverse(breaks = seq(0, 0.4, 0.2)) +
  coord_cartesian(ylim = c(0, 0.5)) +
  scale_fill_viridis_d(option = "plasma", begin = 0.1, end = 0.9) +
  guides(fill = FALSE) +
  labs(x = NULL, y = "Landscape fitness, linked",
       caption = "Point = mean value") +
  facet_wrap(~ constraint_clean, scales = "free_x") +
  theme_bw()

Nested interactions

Landscape fitness, linked

High-low

ggplot(BHL, aes(x = compete, y = landscape_fitness_linked)) +
  geom_violin(aes(fill = repp), size = 0.1, width = 1,
              position = position_dodge(width = 1)) +
  stat_summary(aes(group = repp), geom = "point", fun = "mean", 
               position = position_dodge(width = 1),
               size = 2, pch = 21, fill = "black", color = "white") +
  scale_y_reverse(breaks = c(0, 0.2, 0.4)) +
  coord_cartesian(ylim = c(0, 0.5)) +
  labs(x = NULL, y = "Landscape fitness, linked", 
       fill = NULL, caption = "Point = mean value") +
  facet_grid(selectfor_d + create_network ~ disperse + catastrophe) +
  theme_bw() +
  theme(legend.position = "top",
        legend.key.size = unit(0.65, "lines"))

Latin squares

ggplot(BLS, aes(x = compete, y = landscape_fitness_linked)) +
  geom_violin(fill = "#FE7F2D", size = 0.25) +
  stat_summary(geom = "point", fun = "mean", 
               size = 2, pch = 21, fill = "black", color = "white") +
  scale_y_reverse(breaks = c(0, 0.2, 0.4)) +
  coord_cartesian(ylim = c(0, 0.5)) +
  labs(x = NULL, y = "Landscape fitness, linked",
       caption = "Point = mean value") +
  facet_grid(selectfor_d + create_network ~ disperse + catastrophe) +
  theme_bw()

Modeling + marginal effects

Landscape fitness, linear model

model1 <- lm(landscape_fitness_linked ~ create_network + disperse + compete + 
               catastrophe + selectfor_d + select + numlink + cdis1 + 
               var_fitness + max_intial_proportion_links + maxsp,
             data = BLS)
summary(model1)

## 
## Call:
## lm(formula = landscape_fitness_linked ~ create_network + disperse + 
##     compete + catastrophe + selectfor_d + select + numlink + 
##     cdis1 + var_fitness + max_intial_proportion_links + maxsp, 
##     data = BLS)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.28286 -0.05665 -0.00188  0.05300  0.43085 
## 
## Coefficients:
##                                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                       -3.910e-02  1.210e-03 -32.307  < 2e-16 ***
## create_networkNo network creation  8.311e-03  4.950e-04  16.789  < 2e-16 ***
## disperseNo dispersal               3.427e-05  4.950e-04   0.069   0.9448    
## competeNo competition              1.081e-01  4.950e-04 218.293  < 2e-16 ***
## catastropheNo catastrophe          1.335e-02  4.950e-04  26.965  < 2e-16 ***
## selectfor_dNo selection            2.361e-03  4.950e-04   4.769 1.86e-06 ***
## selectNo selection                 1.608e-01  4.950e-04 324.736  < 2e-16 ***
## numlink                            6.856e-04  1.396e-04   4.913 9.00e-07 ***
## cdis1                             -1.237e-04  7.155e-05  -1.729   0.0838 .  
## var_fitness                        1.328e-01  2.150e-03  61.771  < 2e-16 ***
## max_intial_proportion_links       -2.173e-04  9.540e-04  -0.228   0.8198    
## maxsp                             -9.868e-04  5.001e-05 -19.734  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.06261 on 63988 degrees of freedom
## Multiple R-squared:  0.7151, Adjusted R-squared:  0.7151 
## F-statistic: 1.46e+04 on 11 and 63988 DF,  p-value: < 2.2e-16

huxreg(model1)

	(1)
(Intercept)	-0.039 ***
	(0.001)
create_networkNo network creation	0.008 ***
	(0.000)
disperseNo dispersal	0.000
	(0.000)
competeNo competition	0.108 ***
	(0.000)
catastropheNo catastrophe	0.013 ***
	(0.000)
selectfor_dNo selection	0.002 ***
	(0.000)
selectNo selection	0.161 ***
	(0.000)
numlink	0.001 ***
	(0.000)
cdis1	-0.000
	(0.000)
var_fitness	0.133 ***
	(0.002)
max_intial_proportion_links	-0.000
	(0.001)
maxsp	-0.001 ***
	(0.000)
N	64000
R2	0.715
logLik	86527.063
AIC	-173028.127
*** p < 0.001; ** p < 0.01; * p < 0.05.

typical_values <- BLS %>% 
  summarize_all(typical) %>% 
  mutate(index = 1)

predictions1 <- crossing(numlink = 0:25, create_network = levels(BLS$create_network)) %>% 
  mutate(index = 1) %>%
  left_join(select(typical_values, -numlink, -create_network), by = "index") %>%
  select(-index) %>%
  augment(model1, newdata = ., se_fit = TRUE) %>% 
  mutate(pred = .fitted,
         pred.lower = pred + (qnorm(0.025) * .se.fit),
         pred.upper = pred + (qnorm(0.975) * .se.fit))

ggplot(predictions1, aes(x = numlink, y = pred, color = create_network)) +
  geom_ribbon(aes(ymin = pred.lower, ymax = pred.upper, fill = create_network),
              alpha = 0.3, color = NA) +
  geom_line(size = 1.5) + 
  scale_y_reverse(labels = percent) +
  scale_color_viridis_d(option = "plasma", end = 0.8) +
  scale_fill_viridis_d(option = "plasma", end = 0.8) +
  labs(y = "Predicted landscape fitness", x = "Number of links", color = NULL) + 
  guides(fill = FALSE) +
  theme_bw() +
  theme(legend.position = "bottom")

Landscape fitness, zero-inflated beta model

Histogram of landscape fitness to show zero-inflatedness + constraint within 0-1:

ggplot(BLS, aes(x = landscape_fitness_linked)) +
  geom_histogram(bins = 100, color = "white") +
  scale_x_continuous(labels = percent) +
  theme_bw()

Build model

model_zinb <- readRDS(here("data", "derived_data", "model_zinb_test.rds"))

summary(model_zinb)

##  Family: zero_inflated_beta 
##   Links: mu = logit; phi = identity; zi = identity 
## Formula: landscape_fitness_linked ~ create_network + disperse + compete + catastrophe + selectfor_d + select + numlink + cdis1 + var_fitness + max_intial_proportion_links + maxsp 
##    Data: BLS_small (Number of observations: 16000) 
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
##          total post-warmup samples = 8000
## 
## Population-Level Effects: 
##                                 Estimate Est.Error l-95% CI u-95% CI Rhat
## Intercept                          -4.32      0.02    -4.36    -4.28 1.00
## create_networkNonetworkcreation     0.12      0.01     0.10     0.13 1.00
## disperseNodispersal                -0.03      0.01    -0.04    -0.01 1.00
## competeNocompetition                1.17      0.01     1.16     1.19 1.00
## catastropheNocatastrophe            0.15      0.01     0.14     0.17 1.00
## selectfor_dNoselection              0.02      0.01     0.00     0.03 1.00
## selectNoselection                   1.99      0.01     1.97     2.00 1.00
## numlink                             0.01      0.00     0.00     0.01 1.00
## cdis1                              -0.00      0.00    -0.00     0.00 1.00
## var_fitness                         1.64      0.03     1.58     1.71 1.00
## max_intial_proportion_links        -0.04      0.01    -0.07    -0.01 1.00
## maxsp                              -0.01      0.00    -0.01    -0.01 1.00
##                                 Bulk_ESS Tail_ESS
## Intercept                           7333     6271
## create_networkNonetworkcreation    10443     6134
## disperseNodispersal                10740     6122
## competeNocompetition                7064     6428
## catastropheNocatastrophe            9366     5392
## selectfor_dNoselection             10848     5493
## selectNoselection                   5955     5751
## numlink                            15260     4808
## cdis1                              10038     4993
## var_fitness                         8490     5657
## max_intial_proportion_links        10226     6012
## maxsp                              13425     5594
## 
## Family Specific Parameters: 
##     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## phi    46.89      0.55    45.81    47.99 1.00     6684     5976
## zi      0.00      0.00     0.00     0.00 1.00     9224     5416
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

plot(marginal_effects(model_zinb))

## Warning: Method 'marginal_effects' is deprecated. Please use
## 'conditional_effects' instead.

Copyright © 2019–2020 Steven L. Peck and Andrew Heiss. Licensed under the MIT License. Code available on GitHub.