# Statistical analysis of population differences in post-winter development # parameters # This script contains statistical analyses of post-winter development in the # orange tip butterfly, Anthocharis cardamines, from six different populations # using linear models with planned comparisons. Locality and temperature # treatment are both used as factors for this analysis. Initial pupal weight is # included in all models as a covariate. # Parameters from Bayesian MCMC model fitting (usind Stan) of two-phase # post-winter development dat <- read.delim("postdia_dev_r1_r2.txt") # Construct locality and temperature treatment factors dat$flocality <- factor(dat$locality) dat$ftemp <- factor(dat$temp_treatment) # Subset data to females dat <- subset(dat, sex==1) # Clear empty levels (males) dat <- droplevels(dat) # Weight loss data for all individuals ind_dat <- read.delim("postdia_data.txt") # Subset data to females ind_dat <- subset(ind_dat, sex==1) # Select first weight recorded dat$w0_raw <- ind_dat$weight[!duplicated(ind_dat$ind_id)] # Define comparison contrasts # order of locality levels is: # 1 - N. Swe # 2 - S. UK # 3 - C. UK # 4 - N. UK # 5 - C. Swe # 6 - S. Swe # Define contrasts. These comparisons equate to: # Swe vs. UK # C. Swe vs. N. & S. Swe # N. Swe vs. S. Swe # C. UK vs. N. & S. UK # N. UK vs. S. UK contr <- cbind(c(-1, 1, 1, 1, -1, -1), c(-1, 0, 0, 0, 2, -1), c(1, 0, 0, 0, 0, -1), c(0, -1, 2, -1, 0, 0), c(0, -1, 0, 1, 0, 0)) # Load needed packages library(car) # Log transform overall rate of pupal development mod_tr_r <- lm(log(r) ~ flocality + ftemp + w0_raw, contrast=list(flocality=contr), data=dat) print(summary(mod_tr_r)) print(Anova(mod_tr_r)) # Log transform rate of developmentalin first phase (r1) mod_r1_tr <- lm(log(r1) ~ flocality + ftemp + w0_raw, contrast=list(flocality=contr), data=dat) print(summary(mod_r1_tr)) print(Anova(mod_r1_tr)) # Log transform rate of developmental in second phase (r2) mod_r2_tr <- lm(log(r2) ~ flocality + ftemp + w0_raw, contrast=list(flocality=contr), data=dat) print(summary(mod_r2_tr)) print(Anova(mod_r2_tr)) # Proportion of development in phase 1 (transformed to achieve homogeneous # variances) mod_p1_tr <- lm(p1_ave^4 ~ flocality * ftemp + w0_raw, contrast=list(flocality=contr), data=dat) print(summary(mod_p1_tr)) print(Anova(mod_p1_tr))