#=============================================================================== # Body size & aggression in field crickets #=============================================================================== # Load the necessary packages library(rstan) library(lme4) library(brms) library(shinystan) library(bayesplot) library(rstanarm) library(fitdistrplus) library(ggplot2) library(dplyr) library(tidyr) library(stringr) #=============================================================================== # Load in the data Subject_1_data <- read.csv("Subject_1_data.csv", header = TRUE, sep = ";", dec=",") Subject_2_data <- read.csv("Subject_2_data.csv", header = TRUE, sep = ";", dec=",") datafemur <- read.csv("Data_crickets.csv", header = TRUE, sep = ";", dec=",") Data_subject1_2 <- merge(x=Subject_1_data, y=Subject_2_data, by.x="Assay", by.y="Assay") # y=focal, x=opponent Data_merged1 <- merge(x=Data_subject1_2, y=datafemur, by.x="Animal_ID.x", by.y="Animal_ID") df <- merge(x=Data_merged1, y=datafemur, by.x="Animal_ID.y", by.y="Animal_ID") #=============================================================================== # Prepare the data for the stan model n_distinct(df$Trial.x) #number of unique trials # Turn binary predictors & random effects into integer indexes df <- df %>% mutate(Focal_ID = Animal_ID.y %>% as.factor() %>% as.numeric()) df <- df %>% mutate(Opponent_ID = Animal_ID.x %>% as.factor() %>% as.numeric()) #Include unique ID for each dyad df <- df %>% unite("Dyad_ID", c("Trial.x", "Animal_ID.y", "Animal_ID.x"), sep = "_", remove = FALSE, na.rm = TRUE) df <- df %>% mutate(Dyad_ID_index = Dyad_ID %>% as.factor() %>% as.numeric()) df$Bin_focal_size <-ifelse(df$size.y=="S", -0.5, 0.5) df$Bin_opp_size <-ifelse(df$size.x=="S", -0.5, 0.5) df$Bin_focal_size <- as.numeric(df$Bin_focal_size) df$Bin_opp_size <- as.numeric(df$Bin_opp_size) #=============================================================================== # Final data table also including opponent's data as focal data stan_data_final <- list( n_obs = nrow(df)*2, n_ind = length(unique(c(df$Focal_ID, df$Opponent_ID))), n_dyad = max(df$Dyad_ID_index), individual = c(df$Focal_ID, df$Opponent_ID), opponent = c(df$Opponent_ID, df$Focal_ID), dyad = c(rep(df$Dyad_ID_index, 2)), size_focal = c(df$Bin_focal_size, df$Bin_opp_size), size_opp = c(df$Bin_opp_size, df$Bin_focal_size), z = as.vector(scale(log(c(df$Relative_1moveto.y,df$Relative_1moveto.1.x)))), z_log_sd = sd(log(c(df$Relative_1moveto.y,df$Relative_1moveto.1.x))), z_log_mu = mean(log(c(df$Relative_1moveto.y,df$Relative_1moveto.1.x))) ) #=============================================================================== # stan model #1: focal body size included write( temp <- "data { //Phenotype dataset // Number of clusters (an integer) int n_obs; // number of observations for phenotypes or number of rows int n_ind; // number of individuals int n_dyad; // number of dyads // Clusters identifiers (an integer) int individual[n_obs]; // Individual ID repeated obs int opponent[n_obs]; // Individual ID opponent repeated obs int dyad[n_obs]; // Dyad ID // Predictors vector[n_obs] size_opp; // fixed effect of opponent size vector[n_obs] size_focal; // fixed effect of focal size // Continuous variables vector[n_obs] z; // phenotypic observations // Parameters for backtransformation real z_log_sd; // log standard deviation of phenotype real z_log_mu; // log mean of phenotype } parameters { // Define parameters to estimate // Parameters for model of phenotype real B_0; //intercept real B_1; // fixed effect focal size real B_2; // interaction focal-opp size real psi; //slope; social effect partner size real B_0e; // mean residual real beta_res; // difference in variance // Random effects matrix[n_ind,3] zI; //(intercepts and slopes for each individual) vector[3] sigma_I; // sd intercepts and slopes cholesky_factor_corr[3] L; // factor to estimate covariance int-slopes vector[n_dyad] zDyad; // real sigma_dyad; // } transformed parameters{ matrix[n_ind,3] I; // Unscaled blups intercept and slope and res_impact for each individual vector[n_dyad] effect_dyad; // effect_dyad = zDyad * sigma_dyad; // get the unscaled value I = zI * diag_pre_multiply(sigma_I, L)'; // get the unscaled value } model{ // Priors B_0 ~ normal(0, 1); B_1 ~ normal(0, 1); B_2 ~ normal(0, 1); psi ~ normal(0, 1); B_0e ~ normal(0, 1); beta_res ~ normal(0, 1); // Random effects distribution to_vector(zI) ~ normal(0, 1); sigma_I ~ normal(0, 1); to_vector(zDyad) ~ normal(0, 1); // sigma_dyad ~ normal(0, 1); // L ~ lkj_corr_cholesky(3); vector[n_obs] e_z; // predicted values for phenotype e_z = B_0 + B_1 * size_focal + B_2 .* size_focal .* size_opp + I[individual, 1] + (psi + I[individual, 2]) .* size_opp + I[opponent, 3] + effect_dyad[dyad]; // Heterogeneous residual SD vector[n_obs] sigma_e = sqrt(exp(B_0e + beta_res * size_opp)); // ensures positivity z ~ normal(e_z, sigma_e); } generated quantities{ // Variances (incl bactransformation) real Sigma2_intercept = sigma_I[1]^2; real Sigma2_psi = sigma_I[2]^2; real Sigma2_res_impact = sigma_I[3]^2; real Sigma2_size_focal = variance(size_focal); real Sigma2_size_opp = variance(size_opp); real Sigma2_dyad = sigma_dyad^2; real Sigma2_res = exp(B_0e); // residual variance at midpoint (same as sigma_e^2) real Sigma2_res_small = exp(B_0e - 0.5 * beta_res); real Sigma2_res_large = exp(B_0e + 0.5 * beta_res); real Sigma2_total = Sigma2_intercept + B_1^2*Sigma2_size_focal + Sigma2_psi*Sigma2_size_opp + psi^2*Sigma2_size_opp + Sigma2_res_impact + Sigma2_dyad + Sigma2_res; // Repeatabilities real VC_dyad = Sigma2_dyad/Sigma2_total; // Dyad repeatability real VC_intercept = Sigma2_intercept/Sigma2_total; // Average aggressiveness repeatability real VC_psi = Sigma2_psi*Sigma2_size_opp/Sigma2_total; // responsiveness repeatabiltiy real VC_res_impact = Sigma2_res_impact/Sigma2_total; // Social impact repeatability // Backtransform estimates // To log scale // Main effects real B_0_e = B_0 * z_log_sd + z_log_mu; // intercept real B_1_e = B_1 * z_log_sd; // focal body size real B_2_e = B_2 * z_log_sd; // interaction real psi_e = psi * z_log_sd; // slope to opponent body size // To y-scale // Main effects real B_0_y = exp(B_0_e + 0.5*(Sigma2_res*z_log_sd^2)); // intercept real B_1_y = exp(B_0_e + 0.5*(Sigma2_res*z_log_sd^2)) * (exp(B_1_e)-1); // focal body size //real B_2_y = ; // interaction real psi_y = exp(B_0_e + 0.5*(Sigma2_res*z_log_sd^2)) * (exp(psi_e)-1); // slope to opponent body size // BLUPS matrix[n_ind,3] I_y; // Individual BLUPS I_y[,1] = exp((B_0_e + I[,1] * z_log_sd) + 0.5*(Sigma2_res*z_log_sd^2)); I_y[,2] = exp(B_0_e + 0.5*(Sigma2_res*z_log_sd^2) + I[,1]) .* (exp(psi_e + I[,2] * z_log_sd) -1); I_y[,3] = exp((B_0_e + I[,3] * z_log_sd) + 0.5*(Sigma2_res*z_log_sd^2)); // Variances real Sigma2_intercept_y = (exp(Sigma2_intercept*z_log_sd^2)-1) * exp((2*B_0_e) + Sigma2_intercept*z_log_sd^2); real Sigma2_psi_y = (exp(Sigma2_psi*z_log_sd^2)-1) * exp((2*B_0_e) + Sigma2_psi*z_log_sd^2); real Sigma2_res_impact_y = (exp(Sigma2_res_impact*z_log_sd^2)-1) * exp((2*B_0_e) + Sigma2_res_impact*z_log_sd^2); real Sigma2_dyad_y = (exp(Sigma2_dyad*z_log_sd^2)-1) * exp((2*B_0_e) + Sigma2_dyad*z_log_sd^2); // real Sigma2_res_y = (exp(Sigma2_res*z_log_sd^2)-1) * exp((2*B_0_e) + Sigma2_res*z_log_sd^2); // residual variance at midpoint real Sigma2_res_small_y = (exp(Sigma2_res_small*z_log_sd^2)-1) * exp((2*B_0_e) + Sigma2_res_small*z_log_sd^2); real Sigma2_res_large_y = (exp(Sigma2_res_large*z_log_sd^2)-1) * exp((2*B_0_e) + Sigma2_res_large*z_log_sd^2); real Sigma2_total_y = (exp(Sigma2_total*z_log_sd^2)-1) * exp((2*B_0_e) + Sigma2_total*z_log_sd^2); }", file="RR_mod1.stan") # Compile Stan model for faster fitting md1 <- rstan::stan_model(file = "RR_mod1.stan") # Fit the stan model fit_md1 <- rstan::sampling(md1, data = stan_data_final, chains = 4, iter = 6000, warmup = 2000, thin = 1, cores = 4, save_warmup = FALSE, seed = 458069796 ) rstan::get_seed(fit_md1)#458069796 params1 <- c("B_0", "B_1","B_2", "psi", "B_0_e", "B_1_e","B_2_e", "psi_e", "B_0_y", "B_1_y", "psi_y", "VC_dyad", "VC_intercept", "VC_psi", "VC_res_impact", "Sigma2_dyad", "Sigma2_intercept", "Sigma2_psi", "Sigma2_res_impact", "Sigma2_dyad_y", "Sigma2_intercept_y", "Sigma2_psi_y", "Sigma2_res_impact_y", "Sigma2_res", "Sigma2_res_small", "Sigma2_res_large", "Sigma2_total", "Sigma2_res_y", "Sigma2_res_small_y", "Sigma2_res_large_y", "Sigma2_total_y" ) round(summary(fit_md1, pars=params1)$summary[, c(4, 6, 8, 9, 10)],5) #saveRDS(fit_md1, file = paste("fit_md1.rds", sep = "")) # Save model as RDS file #=============================================================================== # stan model #2: focal body size excluded write( temp <- "data { //Phenotype dataset // Number of clusters (an integer) int n_obs; // number of observations for phenotypes or number of rows int n_ind; // number of individuals int n_dyad; // number of dyads // Clusters identifiers (an integer) int individual[n_obs]; // Individual ID repeated obs int opponent[n_obs]; // Individual ID opponent repeated obs int dyad[n_obs]; // Dyad ID // Predictors vector[n_obs] size_opp; // fixed effect of opponent size // Continuous variables vector[n_obs] z; // phenotypic observations // Parameters for backtransformation real z_log_sd; // log standard deviation of phenotype real z_log_mu; // log mean of phenotype } parameters { // Define parameters to estimate // Parameters for model of phenotype real B_0; //intercept real psi; //slope; social effect partner size real B_0e; // mean residual real beta_res; // difference in variance // Random effects matrix[n_ind,3] zI; //(intercepts and slopes for each individual) vector[3] sigma_I; // sd intercepts and slopes cholesky_factor_corr[3] L; // factor to estimate covariance int-slopes vector[n_dyad] zDyad; // real sigma_dyad; // } transformed parameters{ matrix[n_ind,3] I; // Unscaled blups intercept and slope and res_impact for each individual vector[n_dyad] effect_dyad; // effect_dyad = zDyad * sigma_dyad; // get the unscaled value I = zI * diag_pre_multiply(sigma_I, L)'; // get the unscaled value } model{ // Priors B_0 ~ normal(0, 1); psi ~ normal(0, 1); B_0e ~ normal(0, 1); beta_res ~ normal(0, 1); // Random effects distribution to_vector(zI) ~ normal(0, 1); sigma_I ~ normal(0, 1); to_vector(zDyad) ~ normal(0, 1); // sigma_dyad ~ normal(0, 1); // L ~ lkj_corr_cholesky(3); vector[n_obs] e_z; // predicted values for phenotype e_z = B_0 + I[individual, 1] + (psi + I[individual, 2]) .* size_opp + I[opponent, 3] + effect_dyad[dyad]; // Heterogeneous residual SD vector[n_obs] sigma_e = sqrt(exp(B_0e + beta_res * size_opp)); // ensures positivity z ~ normal(e_z, sigma_e); } generated quantities{ // Variances (incl bactransformation) real Sigma2_intercept = sigma_I[1]^2; real Sigma2_psi = sigma_I[2]^2; real Sigma2_res_impact = sigma_I[3]^2; real Sigma2_size_opp = variance(size_opp); real Sigma2_dyad = sigma_dyad^2; // real Sigma2_res = exp(B_0e); // residual variance at midpoint (same as sigma_e^2) real Sigma2_res_small = exp(B_0e - 0.5 * beta_res); real Sigma2_res_large = exp(B_0e + 0.5 * beta_res); real Sigma2_total = Sigma2_intercept + Sigma2_psi*Sigma2_size_opp + psi^2*Sigma2_size_opp + Sigma2_res_impact + Sigma2_dyad + Sigma2_res; // Repeatabilities real VC_dyad = Sigma2_dyad/Sigma2_total; // Dyad repeatability real VC_intercept = Sigma2_intercept/Sigma2_total; // Average aggressiveness repeatability real VC_psi = Sigma2_psi*Sigma2_size_opp/Sigma2_total; // responsiveness repeatabiltiy real VC_res_impact = Sigma2_res_impact/Sigma2_total; // Social impact repeatability // Backtransform estimates // To log scale // Main effects real B_0_e = B_0 * z_log_sd + z_log_mu; // intercept real psi_e = psi * z_log_sd; // slope to opponent body size // To y-scale // Main effects real B_0_y = exp(B_0_e + 0.5*(Sigma2_res*z_log_sd^2)); // intercept real psi_y = exp(B_0_e + 0.5*(Sigma2_res*z_log_sd^2)) * (exp(psi_e)-1); // slope to opponent body size // BLUPS matrix[n_ind,3] I_y; // Individual BLUPS I_y[,1] = exp((B_0_e + I[,1] * z_log_sd) + 0.5*(Sigma2_res*z_log_sd^2)); I_y[,2] = exp(B_0_e + 0.5*(Sigma2_res*z_log_sd^2) + I[,1]) .* (exp(psi_e + I[,2] * z_log_sd) -1); I_y[,3] = exp((B_0_e + I[,3] * z_log_sd) + 0.5*(Sigma2_res*z_log_sd^2)); // Variances real Sigma2_intercept_y = (exp(Sigma2_intercept*z_log_sd^2)-1) * exp((2*B_0_e) + Sigma2_intercept*z_log_sd^2); real Sigma2_psi_y = (exp(Sigma2_psi*z_log_sd^2)-1) * exp((2*B_0_e) + Sigma2_psi*z_log_sd^2); real Sigma2_res_impact_y = (exp(Sigma2_res_impact*z_log_sd^2)-1) * exp((2*B_0_e) + Sigma2_res_impact*z_log_sd^2); real Sigma2_dyad_y = (exp(Sigma2_dyad*z_log_sd^2)-1) * exp((2*B_0_e) + Sigma2_dyad*z_log_sd^2); // real Sigma2_res_y = (exp(Sigma2_res*z_log_sd^2)-1) * exp((2*B_0_e) + Sigma2_res*z_log_sd^2); // residual variance at midpoint real Sigma2_res_small_y = (exp(Sigma2_res_small*z_log_sd^2)-1) * exp((2*B_0_e) + Sigma2_res_small*z_log_sd^2); real Sigma2_res_large_y = (exp(Sigma2_res_large*z_log_sd^2)-1) * exp((2*B_0_e) + Sigma2_res_large*z_log_sd^2); real Sigma2_total_y = (exp(Sigma2_total*z_log_sd^2)-1) * exp((2*B_0_e) + Sigma2_total*z_log_sd^2); }", file="RR_md2.stan") # Compile Stan model for faster fitting md2 <- rstan::stan_model(file = "RR_md2.stan") # Fit the stan model fit_md2 <- rstan::sampling(md2, data = stan_data_final, chains = 4, iter = 6000, warmup = 2000, thin = 1, cores = 4, save_warmup = FALSE, seed = 458069796) rstan::get_seed(fit_md2)#458069796 params2 <- c( "B_0", "psi", "B_0_e", "psi_e", "B_0_y", "psi_y", "VC_dyad", "VC_intercept", "VC_psi", "VC_res_impact", "Sigma2_dyad", "Sigma2_intercept", "Sigma2_psi", "Sigma2_res_impact", "Sigma2_dyad_y", "Sigma2_intercept_y", "Sigma2_psi_y", "Sigma2_res_impact_y", "Sigma2_res", "Sigma2_res_small", "Sigma2_res_large", "Sigma2_total", "Sigma2_res_y", "Sigma2_res_small_y", "Sigma2_res_large_y", "Sigma2_total_y" ) round(summary(fit_md2, pars=params2)$summary[, c(4, 6, 8, 9, 10)],5) #saveRDS(fit_md_NoSkewExcBS, file = paste("fit_md2.rds", sep = "")) # Save model as RDS file