## Hargreaves et al 2019 LOCAL ADAPTATION TO BIOTIC INTERACTIONS ################### ## ## This script includes Data manipulation, Analyses, Figures ## ## code is structured for use with the R Studio document outline tab visible ## Analyses use one of 2 datasets: ## Dataset1 = experimental manipulations of biotic or abiotic enviornment (some studies) ## Dataset2 = most natural conditions from all studies (may include manipulations of env without a control) ############################################################################################## ### setwd("") getwd() library(tidyverse) library(mosaic) library(lmerTest) library(car) library(visreg) library(lattice) library(lsmeans) #have to load last or command lsmeans gets overwritten by some other package library(emmeans) #updated version of lsmeans package citation(package='lmerTest') citation(package='visreg') R.Version() # Explanation of key variables & data #################################################################################### #__experimental treatments ----- #native = author designation of which source population was native or closest source to transplant site #experiment = whether row is from a study that experimentally manipulated the biotic or abiotic environment (ie yes = Dataset1) #'treat_biotic' & 'treat_abiotic' = which studies experimentally manipulate the transplant environment and how: #none = did not experimentally manpulate the biotic (treat_biotic) or abiotic (treat_abiotic) environment for this study x taxon #yes = env was experimentally manipulated & this row of data is from the manipulated biotic (treat_biotic) or abiotic (treat_abiotic) treatment #no = env was experimentally manipulated & this row of data is from the manipulated env, but manipulation does not affect the biotic (treat_biotic) or abiotic (treat_abiotic) env # natural = env was experimentally manipulated and this row of data is from the control (most natural) treatment #'treat_effect' = how an experimental manpulation of the environment was predicted to affect transplant fitness: # positive = treatment expected to increase transplant fitness (eg removal of antagonists) # negative = treament expected to decrease transplant fitness (eg removing mutualists, increasing antagonists) # none = this row not part of expeirmental manipulation of the environment #__nature of manipulations (whether experimental (Dataset1) or not (Dataset2)) ------ #NOTE: non-experimental manipulations that mimicked natural conditions were not counted. eg if common gardens were watered or shaded to replicate naturally wet or shady habitat, irrigation and shading would be 'no'. if transplants were already in natural habitat and watered or shaded to improve survival, irrigation and shading would be 'yes' #fences_cages = were transplants fenced or caged to exclude predators/herbivores? (yes or no) #insecticide_traps = were insecticides or insect traps used to control insect pests? (yes or no) #herbicide_weeding = were neighbouring plants removed or controlled through herbicides or weeding to reduce competition? (yes or no) #shading = were transplants shaded? no=natural, yes=shaded to improve fitness, shaderemoved=shade removal treatment intended to decrease fitness #pots_garden = were transplants in pots or a garden (ie tilled soil)? (yes or no) #soil_type = was soil composition altered (eg plants placed in sterile soil or potting soil)? (yes or no) #irrigation = were transplants watered (beyond the 1st few days of establishment)? (yes or no) #fertilizer = was fertilizer applied? (yes or no) #agfield = were transplant sites in (old) agricultural fields? #__transplant details ------ #year_begin = year transplant cohort established in common garden #duration = years transplants remained in common garden before performance assessed #temp_rep = temporal replicate within study #lifestage = lifestage at which transplant placed in common garden #lifestage_all = all lifestages that contributed to the performance metric #fitness_type: germ = germination/emergence, surv=survival, germ_survival = germ x surv, repro=reproduction, comp=germ_surv x repro or surv x repro #__fitness/performance details ------ #fitness_data_source = where was fitness info obtained #fitness = numerical fitness response #n_indiv = # individuals contributing to fitness metric (eg 5 blocks of 10 plants per block, n_indiv=50) #n_block = # blocks contributing to fitness metric (eg 5 blocks of 10 plants per block, n_block=5) #fit_var = numeric variation associated with 'fitness' #fit_var_type = type of variation measurement (bvar=binomial variance, se = standard error, ci = confidence interval etc) #__locations ------ #continent_site = continent where common garden was #lat_site, long_site = decimal degree coordinates of common garden site #elevation_m_site = elevation of common garden site in masl #geo_dist = geographic distance as the crow flies between this source population and this transplant site #elev_diff = elevational difference between this source population and this transplant site. <0 means source was from lower elevation than site #lat_diff = latitudinal difference between this source population and this transplant site. <0 means source was from lower latitude than site #thermal_dist converts elevation_diffs to their km thermal equivalents, maintains their direction, and adds them to latitudinal differences (*111 to be on a km scale) #__Key data frames --------- # tp - all studiesxtaxa that have local & foreign source as originally defined (1 row per source, all treatments included) # tpexp - Dataset 1 - studies that experimentally manipulate biotic or abitoic env (1 line per source) # tpexpLA - Dataset 1 - tpexp rearranged to calculate local adaptation metrics (each line compares local vs foreign sources) # tpCTfit - Dataset 2 - most natural treatment from all studies, with standardized fitness variables calculated # tpCTLA - Dataset 2 - tpCTfit rearranged to calculate local adaptation metrics (each line compares local vs foreign sources) #DATA PREP - General (generates df tp) ########################################################## # General (generates df tp) #################################### tp.orig <- read.csv("data/Hargreaves etal 2020 AmNat transplant data.csv") summary(tp.orig) dim(tp.orig) #11088 53 #order variables tp.orig$treat_effect <- factor(tp.orig$treat_effect, levels=c('none','positive','negative')) tp.orig$fitness_type <- factor(tp.orig$fitness_type, levels=c('germ','surv','germ_surv','repro','comp')) summary(tp.orig) #check missing locations tp.orig[is.na(tp.orig$geo_dist),c('study','treat_biotic','treat_effect')] #12 NAs Thompson (no location info, author not contactable), Viveros (coordinates for some sites way off in elevation from GoogleEarth & couldnt confirm) #get rid of lifestage_all (denotes all lifestages that went into calculating each fitness metric but gets confusing with lifestage) tp.orig$lifestage_all <- NULL #tidy - reorder & get rid of columns we don't need names(tp.orig) tp0 <- tp.orig[with(tp.orig, order(study, taxon, site, source, fitness_type)), c('study','functional_group','kingdom', 'phylum_division', 'class', 'order', 'family','genus','taxon', "source","site",'native', 'experiment',"treat_biotic","treat_abiotic","treat_effect","treatment_notes", "fences_cages","insecticide_traps","herbicide_weeding","shading",'pots_garden', 'soil_type',"irrigation",'fertilizer','agfield', 'duration','temp_rep','lifestage',"fitness_type","fitness",'n_indiv','n_block','fit_var',"fit_var_type", 'continent_site','lat_site','long_site','elevation_m_site','elev_diff','geo_dist','thermal_dist')] summary(tp0) dim(tp0) #down to a mere 42 columns #fix names for variance tp0$var <- tp0$fit_var_type #create dummy variable tp0$fit_var_type <- as.factor(ifelse(tp0$var=='se' | tp0$var=='SE', 'se', ifelse(tp0$var=='sd' | tp0$var=='SD', 'sd', as.character(tp0$var)))) table(tp0$var, tp0$fit_var_type) tp0$var <- NULL #get rid of dummy variable summary(tp0) #fix taxonomic issues that will matter to random effect levels(tp0$taxon) tp0$g <- tp0$genus tp0$genus <- as.factor(ifelse(tp0$taxon=='Lychnis_flos_cuculi', 'Silene', ifelse(tp0$taxon=='Campanulastrum_americanum', 'Campanula', as.character(tp0$g)))) levels(tp0$genus) tp0$g <- NULL #__designate local foreign (local Y N) sources and get rid of data without both a local and foreign source ------------------------ #local.orig: whatever the authors said was local summary(tp0$geo_dist[tp0$native=='yes']); summary(tp0$elev_diff[tp0$native=='yes']) #but can be >1000 km or >800m away - not very local tp0$local.orig <- tp0$native #keep column for orginal local designation tp0$local.orig <- factor(tp0$local.orig, levels=c('yes','no')) tp0$native <- NULL #drop to avoid confusion unique(tp0$study[tp0$local.orig=='yes']) #159 studies #local100: criteria 1 (used in main analyses): local sources must be <100 km & <100 m elevation from transplant site --- tp0$local100 <- as.factor(ifelse(tp0$local.orig=='yes' & tp0$geo_dist <100 & abs(tp0$elev_diff) <100, 'yes', 'no')) #authors' 'local' & <100 km & <100m elev tp0$local100[tp0$study=='Thompsonetal_NewPhytol_1991'] <- tp0$local.orig[tp0$study=='Thompsonetal_NewPhytol_1991'] #override (missing locations but were all within 50km) tp0[is.na(tp0$geo_dist),-c(2,3,8:27)] #missing location info for Viveros only affects foreign sources summary(tp0) summary(tp0$geo_dist[tp0$local100=='yes']); #make sure Max <100km tp0$local100 <- factor(tp0$local100, levels=c('yes','no')) #how many studies have local source by this definition? unique(tp0$study[tp0$local100=='yes']) #148 studies have local source within 100km & 100m (loose 11 studies) #adjust Garrido 2012 NewPhytol mannually tp0[tp0$study=='Garridoetal_2012_NewPhytol',c('site','source','local.orig','local100','geo_dist')] #misclassified as elev_diff is NA tp0$local100[tp0$study=='Garridoetal_2012_NewPhytol' & tp0$source=='Pedregal'] <- 'no' #common garden >100masl from nearest source #local50: criteria 2 (used in SupInfo): local sources must be <50 km & <100 m elevation from transplant site --- tp0$local50 <- as.factor(ifelse(tp0$local.orig=='yes' & tp0$geo_dist <50 & abs(tp0$elev_diff) <100, 'yes', 'no')) tp0$local50[tp0$study=='Thompsonetal_NewPhytol_1991'] <- tp0$local.orig[tp0$study=='Thompsonetal_NewPhytol_1991'] #override as above tp0$local50 <- factor(tp0$local50, levels=c('yes','no')) summary(tp0$geo_dist[tp0$local50=='yes']) #make sure Max <50km #how many studies have local source by this definition? unique(tp0$study[tp0$local50=='yes']) #145 studies have local source within 50km & 100m (lose 3 more studies) #adjust Garrido 2012 NewPhytol mannually (as above) tp0[tp0$study=='Garridoetal_2012_NewPhytol',c('site','source','local.orig','local100','local50','geo_dist')] tp0$local50[tp0$study=='Garridoetal_2012_NewPhytol' & tp0$source=='Pedregal'] <- 'no' #make column denoting which site-lifestage-temprep-treat combos have both local & foreign source for each definition of local ------ --- summary(tp0) tp1 <- tp0 %>% group_by(study, taxon, site, fitness_type, lifestage, temp_rep, treatment_notes) %>% mutate(localforeignorig = ifelse((any(local.orig %in% "yes") & any(local.orig %in% "no")), "yes", "no"), localforeign100 = ifelse((any(local100 %in% "yes") & any(local100 %in% "no")), "yes", "no"), localforeign50 = ifelse((any(local50 %in% "yes") & any(local50 %in% "no")), "yes", "no")) tp1$localforeign100 = factor(tp1$localforeign100, levels=c('yes','no')) tp1$localforeign50 = factor(tp1$localforeign50, levels=c('yes','no')) tp1$localforeignorig = factor(tp1$localforeignorig, levels=c('yes','no')) summary(tp1[,c('local.orig','local100','local50','localforeignorig','localforeign100','localforeign50')]) #get rid of rows that don't have a local & foreign source by any definition (ie localforeignorig=='no'). These are in the database as the Bontrager study included 1-source-to-many-site designs but we will never use for local adaptation unique(tp1$study[tp1$localforeignorig=='yes']) #159 studies. so 1 study had no local sources at any site... # table(tp1$study, tp1$localforeignorig) #Welketal_PLoSOne_2014 is the culprit # Welk <- tp1[tp1$study=='Welketal_PLoSOne_2014',]; Welk <- droplevels(Welk) # table(Welk$source, Welk$site) #9 sites (caps), 9 sources, mostly fully crossed # table(Welk$source, Welk$local.orig, Welk$site) #looks like should have local & foreign sources at BS & HAL & JEN sites # summary(Welk) #but used multiple taxa... # table(Welk$source, Welk$local.orig, Welk$site, Welk$taxon) #at taxa level there are never local and foriegn genotypes at a site tp2 <- tp1[tp1$localforeignorig=='yes',]; tp2 <- droplevels(tp2) dim(tp2) #8722 #__designate 'best' fitness measure (closest to lifetime fitness) so can filter down---------- tp3 <- tp2 %>% group_by(study, taxon, site, lifestage, temp_rep, treatment_notes) %>% mutate(bestfitness = derivedFactor ( comp = (any(fitness_type %in% "comp")), repro = (any(fitness_type %in% "repro")), germ_surv = (any(fitness_type %in% "germ_surv")), surv = (any(fitness_type %in% "surv")), germ = (any(fitness_type %in% "germ")), .method = "first", .default = NA)) %>% ungroup() names(tp3) tp3$bestfitness <- factor(tp3$bestfitness, levels=c('germ', 'surv', 'germ_surv', 'repro', 'comp')) #compare fitness_type vs bestfitness summary(tp3$fitness_type); #full data set summary(tp3$bestfitness) #should have fewer germ (NB the best metric is in every column for that study, so there are more lines with 'comp' even though there are obviously not more datapoints with 'comp') levels(tp3$study) #159 studies dim(tp3) #8722 #__which variance estimates do we have? ---------- summary(tp3) summary(tp3[,c('fit_var', 'fit_var_type')]) #bvar = binomial variance. 3070 NAs dim(tp3) #of 8722 lines of data 3070/8722 #35% tp3$havevar <- as.factor(ifelse(is.na(tp3$fit_var), 'no','yes')) table(tp3$fitness_type, tp3$havevar) #missing variance for some of every fitness type, >half for germ_surv & comp #formal meta-analyses weight effect sizes by variance. we have some direct estimates, some estimates of SD from which var can be directly calculated (var = sd^2). Have lots of estimates of SE & can derive SD from SE if know n (SE = SD/sqrt(n)), BUT if authors used a blocked design n to calculate SE should be the # blocks, whereas we collected n individuals. If use wrong n will underestimate the variance. eg if SE=10, n blocks=5, nindiv=50: true SD = 10/sqrt(5) = 4.47, but we would calculate SD = 10/sqrt(50) = 1.4. went back to studies but n used to calculate se not always clear, so don't think we can calculate variance #estimate how many missing values are due to authors vs how many from our calculations test <- tp3[tp3$localforeign100=='yes',]; test <- droplevels(test) summary(test$fit_var) #2643 NAs dim(test) #7675 rows data 2643/7675 #34% #__explore taxonomic nestedness ---------- bysp <- aggregate(data=tp3, source ~ taxon + genus +family +order+class +phylum_division +kingdom, function(x) length(unique(x))) bysp <- droplevels(bysp) head(bysp) summary(bysp) #ideally want 5+ levels for every random factor. ##final overall data set before split into those with exp manipulations (Dataset1) and most natural treatment for all studies (Dataset2) summary(tp3) levels(tp3$study) tp <- tp3 #Dataset1(EXP) data prep for studies with controlled manipulation (generate df tpexp) ############################# tpx <- tp #tpx= tp w eXperimental manips. here is where Dataset1 departs from overall data #identify rows (abiotic_cat) that experimentally manipulated ABIOTIC env. #ct=not part of abioic exp manipulation, treat = part of experimental abiotic manipulation tpx$abiotic_cat <-as.character(tpx$treat_abiotic) tpx$abiotic_cat <-ifelse(tpx$abiotic_cat=="natural","ct",tpx$abiotic_cat) tpx$abiotic_cat <-ifelse(tpx$abiotic_cat=="no","ct",tpx$abiotic_cat) tpx$abiotic_cat <-ifelse(tpx$abiotic_cat=="none","ct",tpx$abiotic_cat) tpx$abiotic_cat <-ifelse(tpx$abiotic_cat=="yes","treat",tpx$abiotic_cat) tpx$abiotic_cat <-as.factor(tpx$abiotic_cat) #identify studies (exp.abiotic) that experimentally manipulated ABIOTIC env. abiotic_bystudy <- as.data.frame(ifelse(table(tpx$study,tpx$abiotic_cat)==0,0,1)) abiotic_bystudy$exp.abiotic <- ifelse(apply(abiotic_bystudy,MARGIN=1,FUN=sum)==1,"no","yes") abiotic_bystudy$study <- rownames(abiotic_bystudy) abiotic_bystudy <- abiotic_bystudy[,3:4] abiotic_bystudy #identify rows (biotic_cat) that manipulated BIOTIC env. #ct=not part of bioic exp manipulation, treat = part of experimental biotic manipulation summary(tpx) tpx$biotic_cat <-ifelse(tpx$treat_biotic=="natural" | tpx$treat_biotic=="no" | tpx$treat_biotic=="none", "ct", ifelse(tpx$treat_biotic=="yes", "treat", tpx$treat_biotic)) tpx$biotic_cat <-as.factor(tpx$biotic_cat) #identify studies (exp.biotic) that manipulated BIOTIC env. biotic_bystudy <-as.data.frame(ifelse(table(tpx$study, tpx$biotic_cat)==0,0,1)) biotic_bystudy$exp.biotic <-ifelse(apply(biotic_bystudy,MARGIN=1,FUN=sum)==1,"no","yes") biotic_bystudy$study <-rownames(biotic_bystudy) biotic_bystudy <- biotic_bystudy[,3:4] biotic_bystudy #add new columns to data frame exp_bystudy <- merge(abiotic_bystudy, biotic_bystudy, by="study") #join biotic & abiotic olumns head(exp_bystudy) head(tpx) tpx1 <- merge(tpx, exp_bystudy, by="study") #adds the 'exp.abiotic' & 'exp.biotic' columns to data sheet tpx1$exp.abiotic <-(as.factor(tpx1$exp.abiotic)) tpx1$exp.biotic <-(as.factor(tpx1$exp.biotic)) head(tpx1) summary(tpx1) dim(tpx1) #8722 but some will be eliminated once we apply 'local' cutoff (8567 rows for initial BioRXiv submission) length(levels(tpx1$study)) #159 studies using any def of local (149 studies for BioRXiv) #do any studies exp manipulate both biotic and abiotic environment? only 1 study and doesnt meet stricter definition of local table(tpx1$exp.biotic, tpx1$exp.abiotic) #yes tpx1$exp.both <- as.factor(ifelse(tpx1$exp.biotic=="yes" & tpx1$exp.abiotic=="yes", "both", "no")) unique(tpx1$study[tpx1$exp.both=='both']) #only 1 study does this: PutnamandReich_EcolMonograph_2017 table(tpx1$treat_effect[tpx1$exp.both=='both']) #manipulations both to improve transplant fitness #are any rows part of both biotic & abiotic experimental manipulations? table(tpx1$biotic_cat, tpx1$abiotic_cat) #yes tpx1[tpx1$biotic_cat=='treat' & tpx1$abiotic_cat=='treat',c(1:3)] #all from Putnam which lose when restrict def of local #how many studies do experimental manipulations (biotic OR abiotic env? summary(tpx1) tpx1$exp.any <- as.factor(ifelse(tpx1$exp.biotic=="yes"|tpx1$exp.abiotic=="yes", "any","none")) table(tpx1$exp.any, tpx1$treat_effect) table(tpx1$exp.any, tpx1$exp.both) #if local source can be any distance away test <- tpx1[tpx1$exp.any=='any' & tpx1$localforeignorig=='yes',]; test <- droplevels(test); nlevels(test$study) #25 studies #if local must be within 100km & 100masl test <- tpx1[tpx1$exp.any=='any' & tpx1$localforeign100=='yes',]; test <- droplevels(test); nlevels(test$study) #23 studies (lose Kreyling & Putnam) dim(tp); dim(tpx1) #should be same length #get rid of studies without a controlled experiment (calculate std fitness differently for analyses of 'most natural treatment') tpexp2 <- subset(tpx1, exp.any=="any") tpexp2 <- droplevels(tpexp2) summary(tpexp2) #check that treatments properly assigned test <- tpexp2[tpexp2$treat_biotic=='yes' & tpexp2$treat_effect=='positive',]; test <- droplevels(test); levels(test$treatment_notes) test <- tpexp2[tpexp2$treat_biotic=='yes' & tpexp2$treat_effect=='negative',]; test <- droplevels(test); levels(test$treatment_notes) test <- tpexp2[tpexp2$treat_abiotic=='yes' & tpexp2$treat_effect=='positive',]; test <- droplevels(test); levels(test$treatment_notes) test <- tpexp2[tpexp2$treat_abiotic=='yes' & tpexp2$treat_effect=='negative',]; test <- droplevels(test); levels(test$treatment_notes) #__calculate standardized fitness---------- #calculate av fitness of all sources & sites (NB: includes exp treatments in mean fitness estimate to keep relative fitness btwn 0-1) # first calculate average fitness of all sources at a given site (no longer used) tp_sitemeans <- tpexp2 %>% group_by(study, taxon, site, fitness_type, lifestage, temp_rep) %>% summarize(avg_fit_site = mean(fitness)) tpexp3 <- full_join(tpexp2, tp_sitemeans) summary(tpexp3) # calculate the fitness of each population relative to the best pop at that site # max_fit_site = maximum fitness at site (by any source in any treatment - incl exp treatments to keep relative fitness btwn 0-1) # std_fit_site = fitness / max_fit_site (ie standardized to range from 0 to 1) tpexp4 <- tpexp3 %>% group_by(study, taxon, fitness_type, lifestage, site, temp_rep) %>% mutate(max_fit_site = max(fitness)) %>% mutate(best_pop_at_site = ifelse((fitness == max(fitness, na.rm=TRUE)), "best", "not")) %>% mutate(best_pop_at_site = ifelse(max_fit_site == 0, "not", best_pop_at_site)) %>% ungroup() %>% mutate(std_fit_site = (fitness/max_fit_site)) %>% mutate(std_fit_site = ifelse(max_fit_site == 0, 0, std_fit_site)) # correct if denominator is 0 summary(tpexp4) #To test effect of whether the manipulation was biotic or abiotic, need a column designating ct, biotic, abiotic, or both. tpexp4$treat_type <- as.factor(ifelse(tpexp4$abiotic_cat=='ct' & tpexp4$biotic_cat=='ct', 'ct', ifelse(tpexp4$abiotic_cat=='treat' & tpexp4$biotic_cat=='ct', 'abiotic', ifelse(tpexp4$abiotic_cat=='ct' & tpexp4$biotic_cat=='treat', 'biotic', ifelse(tpexp4$abiotic_cat=='treat' & tpexp4$biotic_cat=='treat', 'both', 'wierd'))))) summary(tpexp4) ##nothing should be 'weird' in treat_type table(tpexp4$abiotic_cat, tpexp4$treat_effect, tpexp4$treat_type) #some abiotic treatments expected to improve performance (treat_effect=pos) table(tpexp4$abiotic_cat, tpexp4$biotic_cat) #some data points have both biotic and abiotic manipulations test <- (tpexp4[tpexp4$treat_effect=='negative',]); test <- droplevels(test) levels(test$treatment_notes) #treatments that negatively affect fitness are both abiotic and biotic #make new variable that designates unique combos of treat_effect and treat_type so can get fully crossed model tpexp4$treat_full <- as.factor(ifelse(tpexp4$treat_effect=='none', 'ct', paste(tpexp4$treat_effect, tpexp4$treat_type, sep='.'))) summary(tpexp4) tpexp4$treat_full <- factor(tpexp4$treat_full, levels=c('ct','positive.biotic','positive.abiotic','positive.both', 'negative.biotic','negative.abiotic')) #order for FigA1 levels(tpexp4$study) #25 studies (but some get excluded when restrict definition of local) dim(tpexp4) #834 #data set for analyses of standardized fitness and for calculating LA parameters below tpexp <- tpexp4 summary(tpexp) #Dataset1(EXP) Create LA data set with paired local & foreign fitness (generate df tpexpLA) ################## #__local_std_fit_site = mean standardized fitness of local sources at site relative to best source at the site ---------- summary(tpexp) #cant use kingdom or phylum or class as random effects as have <5 levels summary(tpexp$fitness) #restricting local sources to be <50 km and <100 masl away from site (used in Appendix A.2) tp_local50 <- tpexp %>% filter(local50=="yes") %>% #only sources within 50km to be considered local (same def as anomaly paper) group_by(study, taxon, functional_group, class, order, family, genus, site, fitness_type, bestfitness, lifestage, temp_rep, treat_effect, treat_type, treat_full, exp.abiotic, exp.biotic) %>% summarize(local50_fit = mean(fitness), local50_std_fit_site = mean(std_fit_site), local50_nindivsum = sum(n_indiv), local50_nblocksum = sum(n_block), local50_nsource = length(source), local50_nindivmean = mean(n_indiv), local50geo_dist = mean(geo_dist)) #should only be 1 value of n summary(tp_local50) #treat_type = both = 0 bc if say sources cant be local if >50km from site have no local sources dim(tpexp) #834 dim(tp_local50) tp_foreign50 <- tpexp %>% filter(local50=="no") %>% group_by(study, taxon, functional_group, class, order, family, genus, site, fitness_type, bestfitness, lifestage, temp_rep, treat_effect, treat_type, treat_full, exp.abiotic, exp.biotic) %>% summarize(foreign50_fit = mean(fitness), foreign50_fit_site = mean(std_fit_site), foreign50_nindivsum = sum(n_indiv), foreign50_nblocksum = sum(n_block), foreign50_nsource = length(source), foreign50_nindivmean = mean(n_indiv)) #should only be 1 value of n dim(tp_foreign50) #253 foreign sources summary(tp_foreign50) tpexp_50 <- full_join(tp_local50, tp_foreign50) summary(tpexp_50) #should have no NAs for any foreign fitness measure (might have NAs for sample size) #restricting local sources to be <100 km and <100 masl from site (used in main analyses) summary(tpexp) test <- tpexp[tpexp$local100=='yes' & tpexp$local50=='no',]; test <- droplevels(test) test #all 6 rows are from Hughesetal_EcolEvol_2017 #View(tpexp[tpexp$study=='Hughesetal_EcolEvol_2017',]) #but each site has 2 local sources so still have same # rows where localforiegn=yes tp_local100 <- tpexp %>% filter(local100=="yes") %>% #only sources within 100km considered local group_by(study, taxon, functional_group, class, order, family, genus, site, fitness_type, bestfitness, lifestage, temp_rep, treat_effect, treat_type, treat_full, exp.abiotic, exp.biotic) %>% summarize(local100_fit = mean(fitness), local100_std_fit_site = mean(std_fit_site), local100_nindivsum = sum(n_indiv), local100_nblocksum = sum(n_block), local100_nsource = length(source), local100_nindivmean = mean(n_indiv), local100geo_dist = mean(geo_dist), local100elev_diff = mean(elev_diff)) summary(tp_local100) dim(tp_local100) #227 - should be same or bigger than tp_local50 as allow more sources to be 'local' tp_local100[is.na(tp_local100$local100geo_dist),] #all NAs from Thompson 1991 so fine # #do we ever have >1 local source at a site? YES # test <- tpexp %>% # filter(local100=="yes") %>% #only sources within 100km to be considered local (same def as anomaly paper) # group_by(study, taxon, functional_group, site, source, fitness_type, bestfitness, lifestage, temp_rep, # treat_effect, treat_type, treat_full, exp.abiotic, exp.biotic) %>% # summarize(local100_fit = mean(fitness)) # dim(tp_local100) #227 # dim(test) #246 - so there ARE some sites with >1 local source tp_foreign100 <- tpexp %>% filter(local100=="no") %>% group_by(study, taxon, functional_group, class, order, family, genus, site, fitness_type, bestfitness, lifestage, temp_rep, treat_effect, treat_type, treat_full, exp.abiotic, exp.biotic) %>% summarize(foreign100_fit = mean(fitness), foreign100_fit_site = mean(std_fit_site), foreign100_nindivsum = sum(n_indiv), foreign100_nblocksum = sum(n_block), foreign100_nsource = length(source), foreign100_nindivmean = mean(n_indiv)) dim(tp_foreign100) #253 - smaller than when local <50 as more sources count as local summary(tp_foreign100) tpexp_100 <- full_join(tp_local100, tp_foreign100) dim(tpexp_100); dim(tp_local100); dim(tp_foreign100) #tpexp = tpforeign if no sites were missing local sources summary(tpexp_100) #26 NAs are rows that have local source by orig definition but not within 100km & 100m #join to dataframe so that all fitness estimates can be calculated relative to mean local fitness tpexp8 <- full_join(tpexp_50, tpexp_100) dim(tpexp8) #253 summary(tpexp8) #treat_full=right order levels(tpexp8$study) #25 studies (if local can be any distance away) #__calculate LA parameters ---------- #assumes that if had bigger n results would tail out to small numbers but not 0. So makes zeros 1% of mean fitness for other category (local or foreign) #used to calculate parameters with more restricted definition of local (within 50km) as well #how many times do we have 0 fitness for local source? summary(tpexp8) summary(tpexp8[tpexp8$local100_fit==0 & !is.na(tpexp8$local100_fit),]) #cases where mean fitness of local source = 0 summary(tpexp8[tpexp8$foreign100_fit==0 & !is.na(tpexp8$foreign100_fit),]) #cases where mean foreign fitness = 0 #get rid of rows where no local source within 100km and 100masl - neve use these tpexp8 <- tpexp8[!is.na(tpexp8$local100_fit),]; tpexp8 <- droplevels(tpexp8) summary(tpexp8) tpexp9 <- tpexp8 %>% #for <50km <100masl def of local mutate(local50_fit_corr = ifelse(local50_fit == 0, foreign50_fit*0.01, local50_fit)) %>% #new correction method mutate(foreign50_fit_corr = ifelse(foreign50_fit== 0, local50_fit*0.01, foreign50_fit)) %>% #new correction method mutate(corr_log_ratio50 = log((local50_fit_corr)/(foreign50_fit_corr))) %>% #if local=foreign, fraction is 1, log=0 #for <100km <100masl def of local mutate(local100_fit_corr = ifelse(local100_fit == 0, foreign100_fit*0.01, local100_fit)) %>% #new correction method mutate(foreign100_fit_corr = ifelse(foreign100_fit== 0, local100_fit*0.01, foreign100_fit)) %>% #new correction method mutate(corr_log_ratio100 = log((local100_fit_corr)/(foreign100_fit_corr))) #if local=foreign, fraction is 1, log=0 hist(tpexp9$corr_log_ratio50) hist(tpexp9$corr_log_ratio100) #get rid of correction columns - dont need tpexp9$local50_fit_corr <- tpexp9$foreign50_fit_corr <- tpexp9$local100_fit_corr <- tpexp9$foreign100_fit_corr <- NULL #rename for analyses tpexpLA <- tpexp9 #create categorical variable (is there LA, yes or no?) tpexpLA$LA50_cat <- as.factor(ifelse(tpexpLA$local50_fit > tpexpLA$foreign50_fit, 'yes', 'no')) tpexpLA$LA100_cat <- as.factor(ifelse(tpexpLA$local100_fit > tpexpLA$foreign100_fit, 'yes', 'no')) #check summary(tpexpLA) #NAs in corr_log_ratio100 should be rows where local_fit_corr is NA or both local & foreign =0 #calculate some options for n (to be used as weighting factor since can't weight by variance) tpexpLA$LA50_nindivsum <- tpexpLA$local50_nindivsum + tpexpLA$foreign50_nindivsum tpexpLA$LA50_nblocksum <- tpexpLA$local50_nblocksum + tpexpLA$foreign50_nblocksum tpexpLA$LA100_nindivsum <- tpexpLA$local100_nindivsum + tpexpLA$foreign100_nindivsum tpexpLA$LA100_nblocksum <- tpexpLA$local100_nblocksum + tpexpLA$foreign100_nblocksum #get rid of some of component n's as wont use them in this data tpexpLA$local50_nindivsum <- tpexpLA$local50_nblocksum <- tpexpLA$local50_nsource <- tpexpLA$local50_nindivmean <- NULL tpexpLA$foreign50_nindivsum <- tpexpLA$foreign50_nblocksum <- tpexpLA$foreign50_nsource <- tpexpLA$foreign50_nindivmean <- NULL tpexpLA$local100_nindivsum <- tpexpLA$local100_nblocksum <- tpexpLA$local100_nsource <- tpexpLA$local100_nindivmean <- NULL tpexpLA$foreign100_nindivsum <- tpexpLA$foreign100_nblocksum <- tpexpLA$foreign100_nsource <- tpexpLA$foreign100_nindivmean <- NULL summary(tpexpLA) #make sure catLA and continuousLA variables match. cases where LA=yes should have corr_log_ratio100 > 0 test <- tpexpLA[!is.na(tpexpLA$LA100_cat) & tpexpLA$LA100_cat=='yes' & tpexpLA$corr_log_ratio100<=0,]; dim(test) test <- tpexpLA[!is.na(tpexpLA$LA100_cat) & !is.na(tpexpLA$corr_log_ratio100) & tpexpLA$LA100_cat=='no' & tpexpLA$corr_log_ratio100 >=0,]; dim(test) #10 cases where corr-log-ratio=0 test <- droplevels(test); summary(test) #from 4 studies summary(tpexpLA) #Dataset2(ALL): data prep - Narrow to most natural treatment of all studies ############################### # last datframe before we subset to experiments was 'tp' - pick up from there #for Dataset2, get rid of rows that are part of an experimental manipulation tpCT1 <- tp[tp$treat_biotic!='yes' & tp$treat_abiotic!='yes',]; tpCT1 <- droplevels(tpCT1) summary(tpCT1[,c(1,9:15)]) #experiment can still be 'yes' (ie controls from experimental manips) but treat_effect should be all 'none' levels(tpCT1$treatment_notes) #all plausible most natural treatments? dim(tpCT1) #8309 levels(tpCT1$study) #still 159 studies will reduce when exclude those without local source from within 100km #__create factors saying what was altered and how it was supposed to affect transplant fitness. ---------- #check shading - normally plants are shaded by other plants and its bad, but in desert env maybe shading is unnatural amelioration of the abiotic env? Lost CarbajalNavarroetal_FrontiersEcolEvol_submitted shde-removed treatment as was part of experimental treatment summary(tpCT1$shading) test <- tpCT1[tpCT1$shading=='yes',]; test <- droplevels(test); levels(test$study) #4 studies used shading summary(test[test$study=='Andersenetal_ForestEcolManag_2008',])#seedlings of high elev rainforest tree protectd from sun. shading=natural summary(test[test$study=='Centeretal_AmJBot_2016',]) #shading to reduce photoinhibition of germination =positive abiotic ammelioration testV <-test[test$study=='Vizcaino-Palomaretal_PLoSOne_2014',]; testV <- droplevels(testV); summary(testV) #"Vizcaino-Palomaretal_PLoSOne_2014" - did experimental manipulation of shade, but deemed that shading was the most natural. summary(tpCT1) #pots_garden = biotic manipulation (eliminates above and below ground competition but also potentially mycorrhizal fungi) #think Hamann study on trees agfield should be coded as 'fertilizer' tpCT1$unCTabiotic <- as.factor( ifelse(tpCT1$irrigation=='yes' | tpCT1$fertilizer=='yes' | tpCT1$soil_type=='yes' | #will put yes for biotic too for this as changes soil biota as well tpCT1$agfield=='yes' | #will put yes for biotic too except for Hamann as changes everything tpCT1$shading=='yes' & (tpCT1$study=='Walteretal_Evolution_2016' | tpCT1$study=='Centeretal_AmJBot_2016'), 'pos.abiotic', 'nat.abiotic')) tpCT1$unCTbiotic <- as.factor( ifelse(tpCT1$fences_cages=='yes' | tpCT1$insecticide_traps=='yes' | tpCT1$herbicide_weeding=='yes' | tpCT1$pots_garden=='yes' | #unless soil is manipulated this primarily affects competition tpCT1$soil_type=='yes' | #this affects good and bad soil biota tpCT1$agfield=='yes' & (tpCT1$study!='Hamannetal_ForestEcolManag_2000'), #AA reckons this agfield ==fertilizer 'positive.biotic', 'nat.biotic')) summary(tpCT1) dim(tpCT1) #8309 #create variable denoting what was manipulated (for abio + bio manip) table(tpCT1$unCTabiotic, tpCT1$unCTbiotic) #only nat or pos - no alterations that might decrease transplant fitness (makes sense) tpCT1$unCTmanip.env <- as.factor(ifelse(tpCT1$unCTabiotic=='nat.abiotic' & tpCT1$unCTbiotic=='nat.biotic', 'nat', ifelse(tpCT1$unCTabiotic=='nat.abiotic' & tpCT1$unCTbiotic=='positive.biotic', 'biotic', ifelse(tpCT1$unCTabiotic=='pos.abiotic' & tpCT1$unCTbiotic=='nat.biotic', 'abiotic', 'both')))) tpCT1$unCTmanip.env <- factor(tpCT1$unCTmanip.env, levels=c('nat','biotic','abiotic','both')) table(tpCT1$unCTmanip.env) #create variable for intended effect of manipulation (natural or ameliorated - no one tried to universally decrease transplant fitness) tpCT1$unCTmanip.effect <- as.factor(ifelse(tpCT1$unCTmanip.env=='nat', 'nat', "ameliorated")) #__remove data points without local (from within 100km and 100 masl) and foreign source at the same site ---------- tpCT2 <- tpCT1[tpCT1$localforeign100=='yes',]; tpCT2 <- droplevels(tpCT2) dim(tpCT1); dim(tpCT2) #goes down to 7323 (7327 w Garrido) summary(tpCT2$fitness_type) # germ surv germ_surv repro comp # 380 3392 601 878 2072 nlevels(tpCT1$study); nlevels(tpCT2$study) #from 158 to 147 studies summary(tpCT2) #should be no NAs for fitness or nindiv, and no 'no's for localforeign100 # Dataset2(ALL): Calculate standardized fitness metrics (1 point per source) (generates df tpCTfit) ############### tpCT3 <- tpCT2 summary(tpCT3) levels(tpCT3$study) # calculate standardized fitness of each source relative to the maximum fitness achieved by any source at that site tpCT4 <- tpCT3 %>% group_by(study, taxon, fitness_type, lifestage, site, temp_rep, unCTmanip.env, unCTmanip.effect, lat_site) %>% mutate(max_fit_site = max(fitness)) %>% mutate(best_pop_at_site = ifelse((fitness == max(fitness, na.rm=TRUE)), "best", "not")) %>% mutate(best_pop_at_site = ifelse(max_fit_site == 0, "not", best_pop_at_site)) %>% ungroup() %>% mutate(std_fit_site = (fitness/max_fit_site)) %>% mutate(std_fit_site = ifelse(max_fit_site == 0, 0, std_fit_site)) # correct if denominator is 0 dim(tpCT4) #7323 56 summary(tpCT4) #check taxonomic info still there summary(tpCT4$std_fit_site) #med .74 mean .63 summary(tpCT4[,-c(3:8,10:27)]) #should be no NAs for fitness or std_fit_site # calculate # of ammeliorations: dont use in analyses (some have bigger effects) but do use in aggregating summary(tpCT4[,c(10:26)]) tpCT4$fences_cages01 <- ifelse(tpCT4$fences_cages=='yes', 1, 0) tpCT4$insecticide_traps01 <- ifelse(tpCT4$insecticide_traps=='yes', 1, 0) tpCT4$herbicide_weeding01 <- ifelse(tpCT4$herbicide_weeding=='yes', 1, 0) tpCT4$shading01 <- ifelse(tpCT4$shading=='yes', 1, 0) tpCT4$pots_garden01 <- ifelse(tpCT4$pots_garden=='yes', 1, 0) tpCT4$soil_type01 <- ifelse(tpCT4$soil_type=='yes', 1, 0) tpCT4$irrigation01 <- ifelse(tpCT4$irrigation=='yes', 1, 0) tpCT4$fertilizer01 <- ifelse(tpCT4$fertilizer=='yes', 1, 0) tpCT4$agfield01 <- ifelse(tpCT4$agfield=='yes', 1, 0) #calculate total sum of manipulations tpCT4$tot_biomanip <- tpCT4$fences_cages01 + tpCT4$insecticide_traps01 + tpCT4$herbicide_weeding01 + tpCT4$pots_garden01 tpCT4$tot_manip <- tpCT4$tot_biomanip + tpCT4$shading01 + tpCT4$soil_type01 + tpCT4$irrigation01 + tpCT4$fertilizer01 + tpCT4$agfield01 summary(tpCT4) #get rid of intermediary columns tpCT4$fences_cages01 <- tpCT4$insecticide_traps01 <- tpCT4$herbicide_weeding01 <- tpCT4$shading01 <- NULL tpCT4$pots_garden01 <- tpCT4$soil_type01 <- tpCT4$irrigation01 <- tpCT4$fertilizer01 <- tpCT4$agfield01 <- NULL #create variables to denote what is being altered tpCT4$herbiv <- as.factor(ifelse(tpCT4$fences_cages=='yes' | tpCT4$insecticide_traps=='yes', 'lessHerbiv', 'nat')) tpCT4$comp <- as.factor(ifelse(tpCT4$herbicide_weeding=='yes' | tpCT4$pots_garden=='yes', 'lessComp', 'nat')) table(tpCT4$herbiv[tpCT4$unCTmanip.env=='biotic'], tpCT4$comp[tpCT4$unCTmanip.env=='biotic']) tpCT4$biotic_manip <- as.factor( ifelse(tpCT4$herbiv=='lessHerbiv' & tpCT4$comp=='nat' & tpCT4$unCTabiotic=='nat.abiotic', 'herbiv', ifelse(tpCT4$herbiv=='nat' & tpCT4$comp=='lessComp' & tpCT4$unCTabiotic=='nat.abiotic', 'compet', ifelse(tpCT4$herbiv=='lessHerbiv' & tpCT4$comp=='lessComp' & tpCT4$unCTabiotic=='nat.abiotic', 'herbivcompet', ifelse(tpCT4$unCTmanip.env=='nat', 'nat', 'complicated'))))) summary(tpCT4) #localforeign100 should always be yes dim(tpCT4) #7323 data points levels(tpCT4$study) #147 studies levels(tpCT4$taxon) #157 taxa #rename dataset to use for analysis tpCTfit <- tpCT4 tpCTfit <- droplevels(tpCTfit) # Dataset2(All): Calculate LA metrics (need to create new df with 1 point per local-foreign comparison) (generates df tpCTLA) ########### summary(tpCTfit) #local50 tp_local50 <- tpCTfit %>% #for definition that local source must be from <100km & <100 masl away filter(local50=="yes") %>% group_by(study, kingdom, phylum_division, class, order, family, genus, taxon, functional_group, site, fitness_type, lifestage, temp_rep, lat_site, unCTmanip.env, unCTmanip.effect, biotic_manip, tot_biomanip, tot_manip, bestfitness) %>% #unCTmanip for tracking summarize(local50_fit = mean(fitness), local50_std_fit_site = mean(std_fit_site), local50_nindivsum = sum(n_indiv), local50_nblocksum = sum(n_block), local50_nsource = length(source), local50geo_dist = mean(geo_dist), local50elev_diff = mean(elev_diff)) summary(tp_local50) summary(tp_local50[tp_local50$study=='Garridoetal_2012_NewPhytol',]) #foreign50 tp_foreign50 <- tpCTfit %>% filter(local50=="no") %>% group_by(study, kingdom, phylum_division, class, order, family, genus, taxon, functional_group, site, fitness_type, lifestage, temp_rep, lat_site, unCTmanip.env, unCTmanip.effect, biotic_manip, tot_biomanip, tot_manip, bestfitness) %>% #unCTmanip for tracking summarize(foreign50_fit = mean(fitness), foreign50_std_fit_site = mean(std_fit_site), foreign50_nindivsum = sum(n_indiv), foreign50_nblocksum = sum(n_block), foreign50_nsource = length(source)) summary(tp_foreign50) #join tp_50 <- full_join(tp_local50, tp_foreign50) #check summary(tp_50) #should have no NAs for foreign50_fit. # #if do, look at them here... usually from GeberandEckhart_Evolution_2005 for no reason I can figure out # tp_50[is.na(tp_50$foreign50_fit),] #yup, central site # test <- tp_50[tp_50$study=='GeberandEckhart_Evolution_2005' & tp_50$taxon=='Clarkia_xantiana_ssp_parviflora' & # tp_50$site=='central',-c(1:7)]; test <- droplevels(test) # summary(test) # View(test[,-c(1:3,7:9,16:20)]) #first time happened was germ, this time surv & comp. so not about the actual data. wtf #local100 tp_local100 <- tpCTfit %>% #for definition that local source must be from <100km & <100 masl away filter(local100=="yes") %>% group_by(study, kingdom, phylum_division, class, order, family, genus, taxon, functional_group, site, fitness_type, lifestage, temp_rep, lat_site, unCTmanip.env, unCTmanip.effect, biotic_manip, tot_biomanip, tot_manip, bestfitness) %>% #unCTmanip for tracking summarize(local100_fit = mean(fitness), local100_std_fit_site = mean(std_fit_site), local100_nindivsum = sum(n_indiv), local100_nblocksum = sum(n_block), local100_nsource = length(source), local100geo_dist = mean(geo_dist), local100elev_diff = mean(elev_diff)) summary(tp_local100) tp_foreign100 <- tpCTfit %>% filter(local100=="no") %>% group_by(study, kingdom, phylum_division, class, order, family, genus, taxon, functional_group, site, fitness_type, lifestage, temp_rep, lat_site, unCTmanip.env, unCTmanip.effect, biotic_manip, tot_biomanip, tot_manip, bestfitness) %>% #unCTmanip for tracking summarize(foreign100_fit = mean(fitness), foreign100_std_fit_site = mean(std_fit_site), foreign100_nindivsum = sum(n_indiv), foreign100_nblocksum = sum(n_block), foreign100_nsource = length(source)) summary(tp_foreign100) tp_100 <- full_join(tp_local100, tp_foreign100) # join dataframes tpCT5 <- full_join(tp_50, tp_100) summary(tpCT5) #should be no NAs for local100_fit or foreign100_fit dim(tpCT5) #1098 44 (for BioRXiv was 1063 x 18) #changes to 1096 rows while sorting out merging problem - inore for now #create resp ratio metric for LA. corr_log_ratio100 makes zeros be 1% of fitness of other treatment (local or foreign) tpCT6 <- tpCT5 %>% mutate(local50_fit_corr = ifelse(local50_fit == 0, foreign50_fit*0.01, local50_fit)) %>% mutate(foreign50_fit_corr = ifelse(foreign50_fit== 0, local50_fit*0.01, foreign50_fit)) %>% mutate(corr_log_ratio50 = log((local50_fit_corr)/(foreign50_fit_corr))) %>% #for <100km & <100masl definition mutate(local100_fit_corr = ifelse(local100_fit == 0, foreign100_fit*0.01, local100_fit)) %>% mutate(foreign100_fit_corr = ifelse(foreign100_fit== 0, local100_fit*0.01, foreign100_fit)) %>% mutate(corr_log_ratio100 = log((local100_fit_corr)/(foreign100_fit_corr))) #get rid of uneeded columns tpCT6$local50_fit_corr <- tpCT6$foreign50_fit_corr <- tpCT6$local100_fit_corr <- tpCT6$foreign100_fit_corr <- NULL #create categorical variable tpCT6$LA50_cat <- as.factor(ifelse(tpCT6$local50_fit > tpCT6$foreign50_fit, 'yes', 'no')) tpCT6$LA100_cat <- as.factor(ifelse(tpCT6$local100_fit > tpCT6$foreign100_fit, 'yes', 'no')) summary(tpCT6) #name dataset used for analyses tpCTLA <- tpCT6 levels(tpCTLA$study) #calculate some options for n tpCTLA$LA50_nindivsum <- tpCTLA$local50_nindivsum + tpCTLA$foreign50_nindivsum tpCTLA$LA50_nblocksum <- tpCTLA$local50_nblocksum + tpCTLA$foreign50_nblocksum tpCTLA$LA100_nindivsum <- tpCTLA$local100_nindivsum + tpCTLA$foreign100_nindivsum tpCTLA$LA100_nblocksum <- tpCTLA$local100_nblocksum + tpCTLA$foreign100_nblocksum summary(tpCTLA) levels(tpCTLA$study) ## DATA SUMMARIZING ############################################ # #call in data so dont hve to go through all the processing # tpexp <- read.csv(file='data/tpexp 19 07 18.csv'); tpexp$X <- NULL # summary(tpexp) # #order variables # tpexp$local100 <- factor(tpexp$local100, levels=c('yes','no')) # tpexp$fitness_type <- factor(tpexp$fitness_type, levels=c('germ','surv','germ_surv','repro','comp')) # tpexp$treat_effect <- factor(tpexp$treat_effect, levels=c('none','positive','negative')) # tpexp$treat_full <- factor(tpexp$treat_full, # levels=c('ct','positive.biotic','positive.abiotic','positive.both','negative.biotic','negative.abiotic')) # # tpexpLA <- read.csv(file='data/tpexpLA 19 07 18.csv'); tpexpLA$X <- NULL # summary(tpexpLA) # #order variables # tpexpLA$fitness_type <- factor(tpexpLA$fitness_type, levels=c('germ','surv','germ_surv','repro','comp')) # tpexpLA$treat_effect <- factor(tpexpLA$treat_effect, levels=c('none','positive','negative')) # tpexpLA$treat_full <- factor(tpexpLA$treat_full, # levels=c('ct','positive.biotic','positive.abiotic','positive.both','negative.biotic','negative.abiotic')) # # tpCTfit <- read.csv(file='data/tpCTfit 19 07 18.csv'); tpCTfit$X <- NULL # summary(tpCTfit) # #order variables # tpCTfit$local100 <- factor(tpCTfit$local100, levels=c('yes','no')) # tpCTfit$fitness_type <- factor(tpCTfit$fitness_type, levels=c('germ','surv','germ_surv','repro','comp')) # tpCTfit$unCTmanip.env <- factor(tpCTfit$unCTmanip.env, levels=c('nat','biotic','abiotic','both')) # #need genus column # # tpCTfit$tx <- tpCTfit$taxon # # separate(tpCTfit, col=tx, into=c('genus','therest'), sep='_') # # summary(tpCTfit) # # tpCTfit$tx <- NULL #splitting didnt work, try merging # tax0 <- read.csv(file="data/tp (final before splitting dataset1v2) 19 07 18.csv"); #updated w taxonomy # summary(tax0) # tax1 <- tax0[,c('kingdom','phylum_division','class','order','family','genus','taxon', #taxonomy info we need # 'site','source','temp_rep','lifestage','fitness_type','n_indiv','n_block','lat_site')] #variables need for matching # summary(tax1) # tax4 <- aggregate(tax1['n_indiv'], by=tax1[c('genus','taxon')], mean) # summary(tax4) # tax4$n_indiv <- NULL # tpCTfit.orig <- tpCTfit # tpCTfit <- merge(tpCTfit, tax4, by=c('taxon'), all.x=T) # dim(tpCTfit.orig); dim(tpCTfit) # summary(tpCTfit) # # tpCTLA <- read.csv(file='data/tpCTLA 19 07 18.csv'); tpCTLA$X <- NULL # summary(tpCTLA) # #order variables # tpCTLA$unCTmanip.env <- factor(tpCTLA$unCTmanip.env, levels=c('nat','biotic','abiotic','both')) # tpCTLA$fitness_type <- factor(tpCTLA$fitness_type, levels=c('germ','surv','germ_surv','repro','comp')) # #add genus # tpCTLA.orig <- tpCTLA # tpCTLA <- merge(tpCTLA.orig, tax4, by=c('taxon'), all.x=T) # dim(tpCTLA.orig); dim(tpCTLA) # summary(tpCTLA) # # summary(tpexp[,c(1:14)]) # summary(tpexp$treat_full) # summary(tpexpLA[,c(1:12)]) #_Abstract-------------------- #__total studies & taxa ------- --- summary(tpCTfit) #all have local source from within 100km and 100masl nlevels(tpCTfit$study) #147 studies nlevels(tpCTfit$taxon) #157 taxa test <- tpCTfit[tpCTfit$unCTmanip.env=='biotic' | tpCTfit$unCTmanip.env=='nat',]; test <- droplevels(test) nlevels(test$study) #126 studies in main paper analyses nlevels(test$taxon) #130 taxa in main paper analyses #__tropical studies ------- --- test <- tpCTfit[!is.na(tpCTfit$lat_site) & abs(tpCTfit$lat_site)<=23.5 & tpCTfit$localforeign100=='yes',]; test <- droplevels(test) summary(test) levels(test$study) #15 studies (from 9 studies for BioRXiv) nlevels(test$study) / nlevels(tpCTfit$study) #7% -> now 10% of studies length(test$fitness) / length(tpCTfit$fitness) #only 2.4% of data rows #_Intro-------------------- #__fig1 (map) caption ---------- #X studies transplanted a local (within 100km and 100masl of planting site) and foreign source to a common field site nstudy <- tp[tp$localforeign100=='yes',]; nstudy <- droplevels(nstudy); nlevels(nstudy$study) #147 studies #of which Y studies also experimentally manipulated the biotic or abiotic environment of transplants with an appropriate control nstudyExp <- tpexp[tpexp$localforeign100=='yes',]; nstudyExp <- droplevels(nstudyExp); nlevels(nstudyExp$study) #23 studies #__table 1 ---------- summary(nstudyExp) table(nstudyExp$exp.abiotic, nstudyExp$treat_full) #Dataset 1 - how many studies do each kind of manipulation? tb <- unique(cbind(as.character(nstudyExp$treat_full), as.character(nstudyExp$study))) summary(tb) #Dataset 1 - how many taxa have each kind of manipulation? tb <- unique(cbind(as.character(nstudyExp$treat_full), as.character(nstudyExp$taxon))) summary(tb) #Dataset 2 - how many studies do each kind of manipulation? summary(tpCTfit) #all manipulations = positive (look at unCTmanip.effect) tb <- unique(cbind(as.character(tpCTfit$unCTmanip.env), as.character(tpCTfit$study))) summary(tb) #Dataset 2 - how many taxa have each kind of manipulation? tb <- unique(cbind(as.character(tpCTfit$unCTmanip.env), as.character(tpCTfit$taxon))) summary(tb) #what are the manipulations? (table footnotes) #biotic positive fences_cges insecticide_traps herbicide_weeding plus any exp manipulations that dont count as one of these? summary(nstudyExp[nstudyExp$treat_full=='positive.biotic', c('treatment_notes')]) #biotic pos #biotic negative summary(nstudyExp[nstudyExp$treat_full=='neg.biotic',c('treatment_notes')]) #biotic neg #abiotic positive - fertilizer + some cases of shading summary(nstudyExp[nstudyExp$treat_full=='pos.abiotic',c('treatment_notes')]) #abiotic pos from experiments summary(nstudy[nstudy$unCTmanip.effect=='ameliorated',c('shading','fertilizer')]) #abiotic pos from uncontrolled manipulations #abiotic negative summary(nstudyExp[nstudyExp$treat_full=='neg.abiotic',c('treatment_notes')]) #abiotic neg #both - pots, agricultural fields, combos #_Methods-------------------- #__1st pp ------- --- #median (mean) distance between source origin and transplant sites all studies #local sources summary(nstudy$geo_dist[nstudy$local100=='yes']) #median = 0 km, mean = 4.96 km #foreign sources summary(nstudy$geo_dist[nstudy$local100=='no']) #median = 230 km, mean = 584 km #dataset of x studies on Z taxa nlevels(nstudy$study) #147 nlevels(nstudy$taxon) #157 #__dataset1 ------- --- #dataset1 description ## #Controlled manipulations were done on... how many taxa per functional group? tb <- unique(cbind(as.character(nstudyExp$functional_group), as.character(nstudyExp$taxon))) summary(tb) #Controlled pos.biotic manipulations were done on... nstudyExp.pb <- nstudyExp[nstudyExp$treat_full=='positive.biotic',]; nstudyExp.pb <- droplevels(nstudyExp.pb) nstudyExp.pb #X studies nlevels(nstudyExp.pb$study) #15 studies nlevels(nstudyExp.pb$taxon) #22 taxa #X taxa per functional group? tb <- unique(cbind(as.character(nstudyExp.pb$functional_group), as.character(nstudyExp.pb$taxon))) summary(tb) #localness of local sources using local100 - full data (23 studies) summary(nstudyExp$geo_dist[nstudyExp$local100=='yes']) #0 to 87km, median = 0 km, mean = 5.1 km #localness of local sources using local100 - ct vs pos.biotic (15 studies) summary(nstudyExp.pb$geo_dist[nstudyExp.pb$local100=='yes']) #0 to 87km, median = 0 km, mean = 5.9 km #dataset2 description ## #total # taxa levels(tpCTfit$study) #147 taxa #taxa where env either natural or only biotic interactions altered summary(tpCTfit) tb <- tpCTfit[tpCTfit$unCTmanip.env=='biotic' | tpCTfit$unCTmanip.env=='nat',]; tb <- droplevels(tb) nlevels(tb$study) #126 studies of nlevels(tb$taxon) #130 taxa #taxa per functional group? tb2 <- unique(cbind(as.character(tb$functional_group), as.character(tb$taxon))) summary(tb2) levels(tb$functional_group) #localness of local sources using local100 #full dataset (147 studies) summary(tpCTfit$geo_dist[tpCTfit$local100=='yes']) #0 to 99km, median = 0 km, mean = 4.96 km #main analyses (126 studies) summary(tb$geo_dist[tb$local100=='yes']) #0 to 99km, median = 0 km, mean = 5.02 km #localness of local sources using local50 #full dataset (144 studies) nstudy50 <- tpCTfit[tpCTfit$localforeign50=='yes',]; nstudy50 <- droplevels(nstudy50); nlevels(nstudy50$study) #144 studies summary(nstudy50$geo_dist[nstudy50$local50=='yes']) #0 to 48.5km, median = 0 km, mean = 2.77 km #main analyses (123 studies) nstudy50 <- tpCTfit[tpCTfit$localforeign50=='yes' & (tpCTfit$unCTmanip.env=='biotic' | tpCTfit$unCTmanip.env=='nat'),]; nstudy50 <- droplevels(nstudy50); nlevels(nstudy50$study) #123 studies summary(nstudy50$geo_dist[nstudy50$local50=='yes']) #0 to 46km, median = 0 km, mean = 2.5 km #dataset1 manipulation ## #data points where adjusted fitness of local or foreign to avoid corr_log_ratio100 of +or- infinity summary(tpexpLA) #have some NAs in local50 fitness as dont always have local source by this definition summary(tpexpLA[tpexpLA$foreign100_fit==0 & tpexpLA$local100_fit!=0,]) #2 data points where all foriegn fitness = 0 summary(tpexpLA[tpexpLA$local100_fit==0 & tpexpLA$foreign100_fit!=0,]) #2 data points (different ones) where local fitness = 0 summary(tpexpLA[tpexpLA$local100_fit==0 & tpexpLA$foreign100_fit==0,]) #5 data points where fitness = 0 all sources #number of rows missing variance summary(tpexp$fit_var) #131 NAs dim(tpexp) #834 131/834 #16% of rows #__dataset2 ------- --- #dataset2 manipulation ## #cases where fitness zero for all sources summary(tpCTfit[tpCTfit$max_fit_site==0,]) #84 data points where max fitness at that site = 0 summary(tpCTLA[tpCTLA$local100_fit==0 & tpCTLA$foreign100_fit==0,]) #18 local-foreign combinations where fitness = 0 all sources #__tropical studies ------- --- test <- tpexp[!is.na(tpexp$lat_site) & abs(tpexp$lat_site)<=23.5 & tpexp$localforeign100=='yes',]; test <- droplevels(test) summary(test) #2 studies. treatments = water sup and wind reduction (Centeretal_AmJBot_2016; Fetcheretal_Biotropica_2000) #View(test[test$experiment=='yes',]) #_Results-------------------- summary(tpCTfit) #__1st pp ------ --- length(levels(tpCTfit$study)) #147 studies total #% studies that have at least some transplants in totally natural env tb <- tpCTfit[tpCTfit$unCTmanip.env=='nat',]; tb <- droplevels(tb); nlevels(tb$study) #63 studies 63/147 #43% tb <- tpCTfit[tpCTfit$unCTmanip.env!='nat',]; tb <- droplevels(tb); nlevels(tb$study) #88 studies altered env for at least some lifestage/site 88/147 #60% # number of studies that used plants levels(tpCTfit$functional_group) tb <- tpCTfit[tpCTfit$functional_group=='annual' | tpCTfit$functional_group=="herbaceous_perennial" | tpCTfit$functional_group=="woody_perennial",] tb <- droplevels(tb); nlevels(tb$study) #135 studies use plants 135/147 #91.8% of studies use plants nlevels(tb$taxon) #144 taxa 144/157 #92% of taxa are plants #which biotic interactions are being altered most often? --- summary(tpCTfit) #how many studies alter competition? bystudyComp <- aggregate(data=tpCTfit, study ~ herbicide_weeding+pots_garden, function(x) length(unique(x))); bystudyComp 27+26+10#63 studies #how many studies alter herbivory? bystudyHerbiv <- aggregate(data=tpCTfit, study ~ fences_cages+insecticide_traps, function(x) length(unique(x))); bystudyHerbiv 43+1+1#45 studies #bystudySoil <- aggregate(data=tpCTfit, study ~ soil_type, function(x) length(unique(x))); bystudySoil #11 studies #how many studies did exp manipulations? nstudyExp <- tpexp[tpexp$localforeign100=='yes',]; nstudyExp <- droplevels(nstudyExp); nlevels(nstudyExp$study) #23 studies #cite studies that do a manipulation AND have a totally natural control ct <- nstudyExp #shorter name cadillac <- ct[ct$fences_cages=='no' & ct$insecticide_traps=='no' & ct$herbicide_weeding=='no' & ct$shading=='no' & ct$pots_garden=='no' & ct$soil_type=='no' & ct$irrigation=='no' & ct$fertilizer=='no' & ct$agfield=='no', ]; cadillac <- droplevels(cadillac); summary(cadillac) levels(cadillac$study) #10 studies # (no longer any studies that manipulate both biotic and abiotic env - Putnam was only one) #_Discussion---------- #_eg of study that protects some transplants from negative BI with totally natural control treatment? bioX <- nstudyExp[nstudyExp$treat_full=='positive.biotic',]; bioX <- droplevels(bioX) levels(bioX$study) #15 studies, included all cadillac but Abdala (neg.biotic) #up the criteria to those that used a truly local source and measure composite fitness? bioX <- nstudyExp[nstudyExp$treat_full=='positive.biotic' & nstudyExp$geo_dist<50 & nstudyExp$fitness_type=='comp',]; bioX <- droplevels(bioX) levels(bioX$study) #7 studies, included all cadillac but Abdala (neg.biotic) #_Appendices ----- #A1 - how local is local? #using author definitions of local (or closest source if authors dont define local) summary(tp$geo_dist[tp$local.orig=='yes']) #can be >1600 km away summary(tp$elev_diff[tp$local.orig=='yes']) #canbe >800 km away #_Fig S1 plot distnace between local sources and transplant site after applying 100 km & 100m cutoff ---------- summary(tp$geo_dist[tp$local100=='yes']) quartz(width=10, height=4) bgdat2 <- adjustcolor('#482677FF', alpha.f=0.65); #colour for dataset 2 (full data) xlim=c(0,100); ylim=c(25,1100); cxlb=1.1 par(mfrow=c(1,2), bty='l', oma=c(2.5,3.5,1,0), mar=c(1,1,1,1)) hist(tp$geo_dist[tp$local100=='yes'], col=bgdat2, xlim=xlim, ylim=ylim, las=1, bty='l', ylab=NA, main=NA, breaks=20) mtext(side=2, line=3.2, text='# local data points per distance category', cex=cxlb) mtext(side=1, line=2.5, text='Geographic distance (km)', cex=cxlb) mtext(side=3, at=100, line=.5, text='Distance from transplant site to origin site of local sources', cex=cxlb+.1) hist(abs(tp$elev_diff[tp$local100=='yes']), col=bgdat2, xlim=xlim, ylim=ylim, las=1, ylab=NA, main=NA, breaks=20) mtext(side=1, line=2.5, text='Elevational distance (masl)', cex=cxlb) ## ANALYSES ####################################################################################################### #__general statistical notes ---------- # "If you're ever comparing models with different fixed effects you should always use ML and never use REML. Otherwise the results are likely to be garbage. – Ben Bolker Apr 7 '14 at 18:19 stackoverflow.com/questions/22892063". "Zuur et al. (2009; PAGE 122) suggest that "To compare models with nested fixed effects (but with the same random structure), ML estimation must be used and not REML." Faraway's (2006) Extending the linear model with R (p. 156): "The reason is that REML estimates the random effects by considering linear combinations of the data that remove the fixed effects. If these fixed effects are changed, the likelihoods of the two models will not be directly comparable." stats.stackexchange.com/questions/116770/ #glmer (used for binary LA (LA100_cat) and standardized fitness (std_fit_site)) automatically fit by maximum likelihood (Laplace Approximation) #lmer (used for LA log ratio (corr_log_ratio100)) autmoatically = fit by REML *BUT* #anova command to compare two mixed models automatically reftis them by ML ## QUESTION 1: is LA detected more often when biotic interactions intact (binary LA)? ############ summary(tpexpLA) #LA100_cat has 0 NAs but corr_log_ratio100 has 5 NAs, bc corr_log_ratio100 excludes 5 cases where fitness = 0 in for local & foreign #include nested random effect to control for phylogeny #cant use weight by n in binomial models unless n = number of bernoulli trials #__Q1 MAIN MS (Fig 2A): Dataset1 studies that experimentally ammeliorate biotic interactions, local <100km (CT vs pos.bio) ----- #2 studies (Abdala & Carbajal) manipulate biotic interactions but with negative fitness effect: make sure ct treatments from these not included in main analyses (will be exp.biotic==yes) table(tpexpLA$study[tpexpLA$exp.biotic=='yes'], tpexpLA$treat_full[tpexpLA$exp.biotic=='yes']) summary(tpexp[,c(1:10)]) #exclude kingdom, phylum & class as have <5 levels LAcatpB100.tx <- glmer(LA100_cat ~ treat_full + (1|study) + (1|taxon) + (1|fitness_type) + (1|order/family/genus), family=binomial, data=tpexpLA, subset=exp.biotic=='yes' & (treat_full=='ct' | treat_full=='positive.biotic') & study!='Abdala-RobertsandMarquis_Oecologia_2007' & #exclude ct treat from neg.biotic studies study!='CarbajalNavarroetal_FrontiersEcolEvol_submitted' & #exclude ct treat from neg.biotic studies !is.na(tpexpLA$corr_log_ratio100)) #exclude 5 cases where local and foreign fit = 0 LAcatpB100null.tx <- glmer(LA100_cat ~ (1|study) + (1|taxon) + (1|fitness_type) + (1|order/family/genus), family=binomial, data=tpexpLA, subset=exp.biotic=='yes' & (treat_full=='ct' | treat_full=='positive.biotic') & study!='Abdala-RobertsandMarquis_Oecologia_2007' & study!='CarbajalNavarroetal_FrontiersEcolEvol_submitted' & !is.na(tpexpLA$corr_log_ratio100)) #no warnings anova(LAcatpB100.tx, LAcatpB100null.tx, test='Chisq') #with taxonomic random effect - identical if don't include it lsmeans(LAcatpB100.tx, pairwise ~ treat_full) #neither significantly diff from 0 #get sample size for Fig caption n <- tpexpLA[!is.na(tpexpLA$LA100_cat) & tpexpLA$exp.biotic=='yes' & (tpexpLA$treat_full=='ct' | tpexpLA$treat_full=='positive.biotic') & tpexpLA$study!='Abdala-RobertsandMarquis_Oecologia_2007' & tpexpLA$study!='CarbajalNavarroetal_FrontiersEcolEvol_submitted' & !is.na(tpexpLA$corr_log_ratio100),]; n <- droplevels(n) length(levels(n$study)) #n = 15 studies summary(n$LA100_cat) #70 no + 85 yes = 155 data points #__Q1 MAIN MS (Fig 2B): Dataset2 most natural treatment from all studies. truly natural vs ammeliorate biotic interactions ----- table(tpCTLA$unCTmanip.env, tpCTLA$LA100_cat) #find LA >50% of time in all cases summary(tpCTLA) #can go all the way up to phylum (models take forever to run). LActcatB100.tx <- glmer(LA100_cat ~ unCTmanip.env + (1|study) + (1|taxon) + (1|fitness_type) + (1|phylum_division/class/order/family/genus), family=binomial, #no warning! data=tpCTLA, subset=(unCTmanip.env=='biotic' | unCTmanip.env=='nat') & !is.na(tpCTLA$corr_log_ratio100)) #make sure logratio & binary consider same data LActcatB100.nulltx <- glmer(LA100_cat ~ (1|study) + (1|taxon) + (1|fitness_type) + (1|phylum_division/class/order/family/genus), family=binomial, data=tpCTLA, subset=(unCTmanip.env=='biotic' | unCTmanip.env=='nat') & !is.na(tpCTLA$corr_log_ratio100)) #no warnings anova(LActcatB100.tx, LActcatB100.nulltx) #Table 3: P=.47 visreg(LActcatB100.tx, xvar='unCTmanip.env') #LA more common biotic > nat (but NS) lsmeans(LActcatB100.tx, ~ unCTmanip.env) #Table 3: both signif >0 so both treatments have signif overall LA lsmeans(LActcatB100.tx, ~ 1) #Table 3: overall measure of local adaptation since interaction is NS: >0 #get sample size for Fig cation n <- tpCTLA[(tpCTLA$unCTmanip.env=='biotic' | tpCTLA$unCTmanip.env=='nat') & !is.na(tpCTLA$LA100_cat) & !is.na(tpCTLA$corr_log_ratio100),] n <- droplevels(n) dim(n) # 958 data points length(levels(n$study)) #n = 126 studies #__Q1 Appendix A.2: Effect of localness of the 'local' sources? (NO qualitative, results slightly more NS) ----- #Ben Bolker: https://stats.stackexchange.com/questions/130707/can-i-include-covariate-as-random-effect-in-glmer #experimental design must support the estimation. if the covariate didn't vary within subjects (e.g. body mass for short-term experiment) across observations, then there's no way to estimate variation in slopes among individuals. I sometimes see confused people trying to fit models containing random effects of the form (1|covariate). This doesn't make sense: while the effect of covariates can vary across levels ((covariate|subject)), the grouping variable (the variable on the right side of the bar) must be categorical, or reasonably interpretable as a categorical variable (e.g. year, for an experiment with observations distributed across years). #model dist to local source as random covariate (Table S.1) #dataset 1 --- --- LAcatpB100gd <- glmer(LA100_cat ~ treat_full + (1|study) + (1|taxon) + (1|fitness_type) + (1|order/family/genus) #mx gr .0064 + (treat_full|local100geo_dist), family=binomial, data=tpexpLA, subset=exp.biotic=='yes' & (treat_full=='ct' | treat_full=='positive.biotic') & study!='Abdala-RobertsandMarquis_Oecologia_2007' & study!='CarbajalNavarroetal_FrontiersEcolEvol_submitted' &!is.na(tpexpLA$corr_log_ratio100)) LAcatpB100gdnull <- glmer(LA100_cat ~ (1|study) + (1|taxon) + (1|fitness_type) + (1|order/family/genus) + (treat_full|local100geo_dist),family=binomial, data=tpexpLA, subset=exp.biotic=='yes' & (treat_full=='ct' | treat_full=='positive.biotic') & study!='Abdala-RobertsandMarquis_Oecologia_2007' & study!='CarbajalNavarroetal_FrontiersEcolEvol_submitted' & !is.na(tpexpLA$corr_log_ratio100)) #no warnings anova(LAcatpB100gd, LAcatpB100gdnull, test='Chisq') #P was 0.11, now P=.295 - so more NS if control for geo_dist summary(LAcatpB100gdnull) #accounts for more variation than taxon lsmeans(LAcatpB100gd, ~ 1) #Table 3: overall measure of local adaptation: NS (no change) #n might be lower - not sure what happens to Thompson 1991 #dataset 2 --- --- summary(tpCTLA) LActcatB100gd <- glmer(LA100_cat ~ unCTmanip.env + (1|study) + (1|taxon) + (1|fitness_type) + (1|order/family/genus) + (unCTmanip.env|local100geo_dist), family=binomial, data=tpCTLA, subset=(unCTmanip.env=='biotic' | unCTmanip.env=='nat') & !is.na(tpCTLA$corr_log_ratio100)) #mxgrd .00138 LActcatB100gd.null <- glmer(LA100_cat ~ (1|study) + (1|taxon) + (1|fitness_type) + (1|order/family/genus) #mxgrd .0033 + (unCTmanip.env|local100geo_dist), family=binomial, data=tpCTLA, subset=(unCTmanip.env=='biotic' | unCTmanip.env=='nat') & !is.na(tpCTLA$corr_log_ratio100)) anova(LActcatB100gd, LActcatB100gd.null) lsmeans(LActcatB100gd, ~ 1) #overall measure of local adaptation: yes (no change) #use stricter def of local: <50km & <100masl (Table S1) #dataset 1 --- --- LAcatpB50 <- glmer(LA50_cat ~ treat_full + (1|study) + (1|taxon) + (1|fitness_type) + (1|order/family/genus), family=binomial, data=tpexpLA, subset=exp.biotic=='yes' & (treat_full=='ct' | treat_full=='positive.biotic') & study!='Abdala-RobertsandMarquis_Oecologia_2007' & study!='CarbajalNavarroetal_FrontiersEcolEvol_submitted' &!is.na(tpexpLA$corr_log_ratio50)) LAcatpB50null <- glmer(LA50_cat ~ (1|study) + (1|taxon) + (1|fitness_type) + (1|order/family/genus), family=binomial, data=tpexpLA, subset=exp.biotic=='yes' & (treat_full=='ct' | treat_full=='positive.biotic') & study!='Abdala-RobertsandMarquis_Oecologia_2007' & study!='CarbajalNavarroetal_FrontiersEcolEvol_submitted' & !is.na(tpexpLA$corr_log_ratio50)) #no warnings anova(LAcatpB50, LAcatpB50null, test='Chisq') #P was 0.11, now P=.22 - so more NS if control for geo_dist lsmeans(LAcatpB50, ~ 1) #overall measure of local adaptation: NS (no change) #dataset 2 --- --- LActcatB50 <- glmer(LA50_cat ~ unCTmanip.env + (1|study) + (1|taxon) + (1|fitness_type) + (1|phylum_division/class/order/family/genus), family=binomial, data=tpCTLA, subset=(unCTmanip.env=='biotic' | unCTmanip.env=='nat') & !is.na(tpCTLA$corr_log_ratio50)) LActcatB50.null <- glmer(LA50_cat ~ (1|study) + (1|taxon) + (1|fitness_type) + (1|phylum_division/class/order/family/genus), family=binomial, data=tpCTLA, subset=(unCTmanip.env=='biotic' | unCTmanip.env=='nat') & !is.na(tpCTLA$corr_log_ratio50)) anova(LActcatB50, LActcatB50.null) lsmeans(LActcatB50, ~ 1) #overall measure of local adaptation: yes (no change) #__Q1 TableS2 - use only best fitness metric (SUPP MAT) ---------- #Dataset 1 (no weights, w taxonomic correction) --- --- LAcatpB100.bestfit <- glmer(LA100_cat ~ treat_full + (1|study) + (1|taxon) + (1|fitness_type) + (1|order/family/genus), family=binomial, data=tpexpLA, subset=exp.biotic=='yes' & (treat_full=='ct' | treat_full=='positive.biotic') & study!='Abdala-RobertsandMarquis_Oecologia_2007' & study!='CarbajalNavarroetal_FrontiersEcolEvol_submitted' &!is.na(tpexpLA$corr_log_ratio100) & as.character(fitness_type)==as.character(bestfitness)) lsmeans(LAcatpB100.bestfit, ~ treat_full) #Table A1: no LA both overlap 0 LAcatpB100.bestfitnull <- glmer(LA100_cat ~ (1|study) + (1|taxon) + (1|fitness_type) + (1|order/family/genus), family=binomial, data=tpexpLA, subset=exp.biotic=='yes' & (treat_full=='ct' | treat_full=='positive.biotic') & study!='Abdala-RobertsandMarquis_Oecologia_2007' & study!='CarbajalNavarroetal_FrontiersEcolEvol_submitted' & !is.na(tpexpLA$corr_log_ratio100) & as.character(fitness_type)==as.character(bestfitness)) anova(LAcatpB100.bestfit, LAcatpB100.bestfitnull, test='Chisq') #Table A1 Q1 Dataset 1 p=.26 #Dataset 2 --- --- LActcatB100.bestf <- glmer(LA100_cat ~ unCTmanip.env + (1|study) + (1|taxon) + (1|fitness_type) + (1|phylum_division/class/order/family/genus), family=binomial, data=tpCTLA, subset=unCTmanip.env!='abiotic' & unCTmanip.env!='both' & !is.na(tpCTLA$corr_log_ratio100) & fitness_type==bestfitness) LActcatB100.bestfnull <- glmer(LA100_cat ~ (1|study) + (1|taxon) + (1|fitness_type) + (1|phylum_division/class/order/family/genus), family=binomial, data=tpCTLA, subset=unCTmanip.env!='abiotic' & unCTmanip.env!='both' & !is.na(tpCTLA$corr_log_ratio100) & fitness_type==bestfitness) #no warnings anova(LActcatB100.bestf, LActcatB100.bestfnull) #Table A1 Q1 Dataset 2 p=.40 lsmeans(LActcatB100.bestf, ~ unCTmanip.env) #both > 0 lsmeans(LActcatB100.bestf, ~ 1) #Table A1: overall measure of local adaptation: >0 #__Q1 SI FIG A1.A&B: all controlled manipulations (CT vs neg.bio vs pos.bio), local <100km ---------- #dataset 1 Fig A1 A------- -- LAcatall100.exp <- glmer(LA100_cat ~ treat_full + (1|study) + (1|taxon) + (1|fitness_type) + (1|order/family/genus), family=binomial, data=tpexpLA, subset=!is.na(tpexpLA$corr_log_ratio100)) LAcatall100.expnull <- glmer(LA100_cat ~ (1|study) + (1|taxon) + (1|fitness_type) + (1|order/family/genus), family=binomial, data=tpexpLA, subset=!is.na(tpexpLA$corr_log_ratio100)) anova(LAcatall100.exp, LAcatall100.expnull) #Chisq=4.3, P=.37 lsmeans(LAcatall100.exp, ~ treat_full) #no LA lsmeans(LAcatall100.exp, ~ 1) #no LA #n for Fig caption n <- tpexpLA[!is.na(tpexpLA$LA100_cat),]; n <- droplevels(n); nlevels(n$study) #23 studies dim(n) # 227 data points #dataset 2 Fig A1 B --- --- #lsmeans & visreg giving different orders of factors. look at raw data tpCTLA$LAcatnum <- ifelse(tpCTLA$LA100_cat=='yes',1,0) bwplot(LAcatnum ~ unCTmanip.env, data=tpCTLA, subset=!is.na(tpexpLA$corr_log_ratio100)) #abio > both > bio > nat test <- aggregate(tpCTLA['LAcatnum'], by=tpCTLA['unCTmanip.env'], mean); test #abio > both > bio > nat LAcatall100.ct <- glmer(LA100_cat ~ unCTmanip.env + (1|study) + (1|taxon) + (1|fitness_type) + (1|phylum_division/class/order/family/genus), family=binomial, data=tpCTLA) #0 warnings LAcatall100.ctnull <- glmer(LA100_cat ~ (1|study) + (1|taxon) + (1|fitness_type) + (1|phylum_division/class/order/family/genus), family=binomial, data=tpCTLA) #no warnings anova(LAcatall100.ct, LAcatall100.ctnull) #Chisq=1.8, P=0.62 visreg(LAcatall100.ct, xvar='unCTmanip.env', scale='response', ylim=c(0,1), rug=F) #'abiotic' poorly replicated lsmeans(LAcatall100.ct, pairwise ~ unCTmanip.env, transform='response') #n for Fig caption n <- tpCTLA$study[!is.na(tpCTLA$LA100_cat)]; n <- droplevels(n); nlevels(n) #147 studies ## QUESTION 2: is LA stronger (effect size) when biotic interactions are left intact? ########### summary(tpexpLA) #__Q2 MAIN MS (Fig 3A): Dataset1 effect size (log ratio), ct vs pos.bio (ignore neg.biotic which is poorly replicated) ---------- #weighting by sqrt(LA100_nindivsum) & with taxonomic random effect summary(tpexpLA) LAposB100.wsqtx <- lmer(corr_log_ratio100 ~ treat_full + (1|study) + (1|taxon) + (1|fitness_type) + (1|order/family/genus), weights=sqrt(LA100_nindivsum), data=tpexpLA, subset=exp.biotic=='yes' & (treat_full=='ct' | treat_full=='positive.biotic') & study!='Abdala-RobertsandMarquis_Oecologia_2007' & study!='CarbajalNavarroetal_FrontiersEcolEvol_submitted') #exclude ct treat from neg.biotic studies LAposB100.null.wsqtx <- lmer(corr_log_ratio100 ~ (1|study) + (1|taxon) + (1|fitness_type) + (1|order/family/genus), weights=sqrt(LA100_nindivsum), data=tpexpLA, subset=exp.biotic=='yes' & (treat_full=='ct' | treat_full=='positive.biotic') & study!='Abdala-RobertsandMarquis_Oecologia_2007'& study!='CarbajalNavarroetal_FrontiersEcolEvol_submitted') anova(LAposB100.wsqtx, LAposB100.null.wsqtx, test='Chisq') #P becomes .378 - slightly less NS lsmeans(LAposB100.wsqtx, ~ treat_full) #LA NS in either treatment lsmeans(LAposB100.wsqtx, ~ 1) #LA NS overall #check sample size (should match Q1 dataset 1) n <- tpexpLA[!is.na(tpexpLA$corr_log_ratio100) & tpexpLA$exp.biotic=='yes' & (tpexpLA$treat_full=='ct' | tpexpLA$treat_full=='positive.biotic') & tpexpLA$study!='Abdala-RobertsandMarquis_Oecologia_2007' & tpexpLA$study!='CarbajalNavarroetal_FrontiersEcolEvol_submitted',] n <- droplevels(n); length(levels(n$study)) #n=15 studies dim(n) #155 data points. #__Q2 MAIN MS (Fig 3B): Dataset2 effect size (log ratio), most nat treatment from all studies, only nat or +biotic -------- #weighting by sqrt(n) and with taxonomic random effect summary(tpCTLA) LAB100ctc.wsqtx <- lmer(corr_log_ratio100 ~ unCTmanip.env + (1|study) + (1|taxon) + (1|fitness_type) + (1|phylum_division/class/order/family/genus), weights=sqrt(LA100_nindivsum), data=tpCTLA, subset=(unCTmanip.env=='biotic' | unCTmanip.env=='nat')) LAB100ctc.null.wsqtx <- lmer(corr_log_ratio100 ~ (1|study) + (1|taxon) + (1|fitness_type) + (1|phylum_division/class/order/family/genus), weights=sqrt(LA100_nindivsum), data=tpCTLA, subset=(unCTmanip.env=='biotic' | unCTmanip.env=='nat')) anova(LAB100ctc.wsqtx, LAB100ctc.null.wsqtx) #Table 3 Q2 Dataset2 effect size p=.043 visreg(LAB100ctc.wsqtx, xvar='unCTmanip.env') lsmeans(LAB100ctc.wsqtx, ~ unCTmanip.env) #neither has signif LA, but LA biotic > LA nat lsmeans(LAB100ctc.wsqtx, ~ 1) #overall 0.067 #get n n <- tpCTLA[!is.na(tpCTLA$corr_log_ratio100) & (tpCTLA$unCTmanip.env=='biotic' | tpCTLA$unCTmanip.env=='nat'),] n <- droplevels(n); length(levels(n$study)) #126 studies dim(n) #958 data points #__Q2 Appendix A.2 Effect of localness of the local source? effect size (log ratio) ----- #using stricter definition of local (<50km and <100masl) #dataset 1 --- --- LAposB50 <- lmer(corr_log_ratio50 ~ treat_full + (1|study) + (1|taxon) + (1|fitness_type) + (1|order/family/genus), weights=sqrt(LA100_nindivsum), data=tpexpLA, subset=exp.biotic=='yes' & (treat_full=='ct' | treat_full=='positive.biotic') & study!='Abdala-RobertsandMarquis_Oecologia_2007' & study!='CarbajalNavarroetal_FrontiersEcolEvol_submitted') LAposB50.null <- lmer(corr_log_ratio50 ~ (1|study) + (1|taxon) + (1|fitness_type)+ (1|order/family/genus), weights=sqrt(LA100_nindivsum), data=tpexpLA, subset=exp.biotic=='yes' & (treat_full=='ct' | treat_full=='positive.biotic') & study!='Abdala-RobertsandMarquis_Oecologia_2007' & study!='CarbajalNavarroetal_FrontiersEcolEvol_submitted') anova(LAposB50, LAposB50.null, test='Chisq') #Table 3 Q2 Dataset1 effect size p was .59, now .58 lsmeans(LAposB50, ~ 1) #no overall LA (so no change) #effect on n main analyses test100 <- tpexpLA[!is.na(tpexpLA$corr_log_ratio100) & tpexpLA$exp.biotic=='yes' & (tpexpLA$treat_full=='ct' | tpexpLA$treat_full=='positive.biotic') & tpexpLA$study!='Abdala-RobertsandMarquis_Oecologia_2007' & tpexpLA$study!='CarbajalNavarroetal_FrontiersEcolEvol_submitted',]; test100 <- droplevels(test100) test50 <- tpexpLA[!is.na(tpexpLA$corr_log_ratio50) & tpexpLA$exp.biotic=='yes' & (tpexpLA$treat_full=='ct' | tpexpLA$treat_full=='positive.biotic') & tpexpLA$study!='Abdala-RobertsandMarquis_Oecologia_2007' & tpexpLA$study!='CarbajalNavarroetal_FrontiersEcolEvol_submitted',]; test50 <- droplevels(test50) nlevels(test100$study); nlevels(test50$study) #14 - lose 1 study dim(test100); dim(test50) #dataset 2 --- --- summary(tpCTLA$phylum_division[!is.na(tpCTLA$corr_log_ratio50)]) #still have 7 phyla summary(tpCTLA$class[!is.na(tpCTLA$corr_log_ratio50)]) #still have 10 classes LAB50ctc <- lmer(corr_log_ratio50 ~ unCTmanip.env + (1|study) + (1|taxon) + (1|fitness_type) + (1|phylum_division/class/order/family/genus), weights=sqrt(LA100_nindivsum), data=tpCTLA, subset=(unCTmanip.env=='biotic' | unCTmanip.env=='nat')) LAB50ctc.null <- lmer(corr_log_ratio50 ~ (1|study) + (1|taxon) + (1|fitness_type) + (1|phylum_division/class/order/family/genus), weights=sqrt(LA100_nindivsum), data=tpCTLA, subset=(unCTmanip.env=='biotic' | unCTmanip.env=='nat')) anova(LAB50ctc, LAB50ctc.null) #DOES CHANGE - from p=.043 to p=.10 even though dont exclude that much data lsmeans(LAB100ctc.wsqtx, ~ unCTmanip.env) #nat=-0.082, bio.amel=0.214 (takes a while to run for some reason) lsmeans(LAB50ctc, ~ unCTmanip.env) #nat=-0.051, bio.amel=0.185 lsmeans(LAB100ctcgd, ~1) #no overall LA (so no change) #effect on n summary(tpCTLA) test100 <- tpCTLA[!is.na(tpCTLA$corr_log_ratio100) & (tpCTLA$unCTmanip.env=='biotic' | tpCTLA$unCTmanip.env=='nat'),]; test100 <- droplevels(test100) test50 <- tpCTLA[!is.na(tpCTLA$corr_log_ratio50) & (tpCTLA$unCTmanip.env=='biotic' | tpCTLA$unCTmanip.env=='nat'),]; test50 <- droplevels(test50) nlevels(test100$study); nlevels(test50$study) #lose 4 of 126 studies (<3%) nlevels(test100$taxon); nlevels(test50$taxon) #lose 3 of 130 taxa (<3%) nlevels(test100$class); nlevels(test50$class) #dont lose any classes nlevels(test100$order); nlevels(test50$order) #dont lose any orders nlevels(test100$family); nlevels(test50$family) #lose 1 family nlevels(test100$genus); nlevels(test50$genus) #lose 2 genera summary(test100$corr_log_ratio100); summary(test100$corr_log_ratio50) #same min and max, mean and median higher. ie LA no longer stronger in biotically-ameliorated studies. So including not-super-local sources might actually be obscuring local adaptation test100$excl50 <- ifelse(is.na(test100$corr_log_ratio50), 'excluded', 'kept') bwplot(corr_log_ratio100 ~ excl50, data=test100) #least local sources have slightly lower corr_log_ratio bwplot(LA100_nindivsum ~ excl50, data=test100) #least local sources have slightly higher n summary(test100[is.na(test100$corr_log_ratio50),]) summary(test100$LA100_nindivsum); summary(test100$LA50_nindivsum) #these studies are also bigger than average so are having relatively strong effects on the results #using random covariate #dataset 1 --- --- LAposB100gd <- lmer(corr_log_ratio100 ~ treat_full + (1|study) + (1|taxon) + (1|fitness_type) + (1|order/family/genus) + #2 warnings (treat_full|local100geo_dist), weights=sqrt(LA100_nindivsum), data=tpexpLA, subset=exp.biotic=='yes' & (treat_full=='ct' | treat_full=='positive.biotic') & study!='Abdala-RobertsandMarquis_Oecologia_2007' & study!='CarbajalNavarroetal_FrontiersEcolEvol_submitted') LAposB100gd.null <- lmer(corr_log_ratio100 ~ (1|study) + (1|taxon) + (1|fitness_type) + (1|order/family/genus) + (treat_full|local100geo_dist), weights=sqrt(LA100_nindivsum), data=tpexpLA, subset=exp.biotic=='yes' & (treat_full=='ct' | treat_full=='positive.biotic') & study!='Abdala-RobertsandMarquis_Oecologia_2007' & study!='CarbajalNavarroetal_FrontiersEcolEvol_submitted') anova(LAposB100gd, LAposB100gd.null, test='Chisq') #doesnt lsmeans(LAposB100gd, ~ 1) #no overall LA (so no change) #dataset 2 --- --- LAB100ctcgd <- lmer(corr_log_ratio100 ~ unCTmanip.env + (1|study) + (1|taxon) + (1|fitness_type) + (1|phylum_division/class/order/family/genus) + (unCTmanip.env|local100geo_dist), weights=sqrt(LA100_nindivsum), data=tpCTLA, subset=(unCTmanip.env=='biotic' | unCTmanip.env=='nat')) #introduces 2 warnings LAB100ctcgd.null <- lmer(corr_log_ratio100 ~ (1|study) + (1|taxon) + (1|fitness_type) + (1|phylum_division/class/order/family/genus) + (unCTmanip.env|local100geo_dist), weights=sqrt(LA100_nindivsum), data=tpCTLA, subset=(unCTmanip.env=='biotic' | unCTmanip.env=='nat')) anova(LAB100ctcgd, LAB100ctcgd.null) #still barely signif lsmeans(LAB100ctcgd, ~1) #no overall LA (so no change) #__Q2 SI Table A1: effect size (logratio) using only best fitness metric ###### #datset 1 --- --- LAposB100.bestfit <- lmer(corr_log_ratio100 ~ treat_full + (1|study) + (1|taxon) + (1|fitness_type) + (1|order/family/genus), weights=sqrt(LA100_nindivsum), data=tpexpLA, subset=exp.biotic=='yes' & (treat_full=='ct' | treat_full=='positive.biotic') & study!='Abdala-RobertsandMarquis_Oecologia_2007' & study!='CarbajalNavarroetal_FrontiersEcolEvol_submitted' & as.character(fitness_type)==as.character(bestfitness)) #no warnings LAposB100.bestfitnull <- lmer(corr_log_ratio100 ~ (1|study) + (1|taxon) + (1|fitness_type) + (1|order/family/genus), weights=sqrt(LA100_nindivsum), data=tpexpLA, subset=exp.biotic=='yes' & (treat_full=='ct' | treat_full=='positive.biotic') & study!='Abdala-RobertsandMarquis_Oecologia_2007' & study!='CarbajalNavarroetal_FrontiersEcolEvol_submitted'& as.character(fitness_type)==as.character(bestfitness)) #2 warnings anova(LAposB100.bestfit, LAposB100.bestfitnull, test='Chisq') #Table A1 Q2 Dataset1 p=.74 lsmeans(LAposB100.bestfit, ~ treat_full) #no signal LA #dataset 2 --- --- LAB100ctc.bf <- lmer(corr_log_ratio100 ~ unCTmanip.env + (1|study) + (1|taxon) + (1|fitness_type) + (1|phylum_division/class/order/family/genus), weights=sqrt(LA100_nindivsum), data=tpCTLA, subset=unCTmanip.env!='abiotic' & unCTmanip.env!='both' & fitness_type==bestfitness) LAB100ctc.bfnull <- lmer(corr_log_ratio100 ~ (1|study) + (1|taxon) + (1|fitness_type)+ #2 warnings (1|phylum_division/class/order/family/genus), weights=sqrt(LA100_nindivsum), data=tpCTLA, subset=unCTmanip.env!='abiotic' & unCTmanip.env!='both' & fitness_type==bestfitness) anova(LAB100ctc.bf, LAB100ctc.bfnull) #now signif (so is main ms analysis so fine) lsmeans(LAB100ctc.bf, ~ unCTmanip.env) #if reporting contrasts for each manip even though manip was NS lsmeans(LAB100ctc.bf, ~ 1) #Table A1 LA: if want 1 overall value #__Q2 SI FIG A1 C&D: effect size all treatments ---------- #dataset 1 Fig A1 C ------- -- LAall100exp <- lmer(corr_log_ratio100 ~ treat_full + (1|study) + (1|taxon) + (1|fitness_type) + (1|order/family/genus), weights=sqrt(LA100_nindivsum), data=tpexpLA) LAall100exp.null<- lmer(corr_log_ratio100 ~ (1|study) + (1|taxon) + (1|fitness_type) + (1|order/family/genus), weights=sqrt(LA100_nindivsum), data=tpexpLA) anova(LAall100exp, LAall100exp.null) #Chisq=1.5, P=.83 visreg(LAall100exp, xvar='treat_full') #totally no difference lsmeans(LAall100exp, ~ treat_full) #CI overlap 0 in all cases - 'both' gone as lost Putnam #dataset 2 Fig A1 D ------- -- LAall100ctc <- lmer(corr_log_ratio100 ~ unCTmanip.env + (1|study) + (1|taxon) + (1|fitness_type) + (1|phylum_division/class/order/family/genus), weights=sqrt(LA100_nindivsum), data=tpCTLA) LAall100ctc.null<- lmer(corr_log_ratio100 ~ (1|study) + (1|taxon) + (1|fitness_type) + (1|phylum_division/class/order/family/genus), weights=sqrt(LA100_nindivsum), data=tpCTLA) anova(LAall100ctc, LAall100ctc.null) #Chisq=4.8, P=0.19 visreg(LAall100ctc, xvar='unCTmanip.env') lsmeans(LAall100ctc, ~ unCTmanip.env) #NS ## QUESTION 2: is LA stronger (std fitness) when biotic interactions are left intact? ########### #after some sleauthing I *think* glm binomial does not automatically round proportional data to integers 0 or 1 (unlike poisson and neg.bin which do round data automatically). quasi binomial does not throw 'non-integer #successes in a binomial glm!' warning but cant do with glmer #__Q2 MAIN MS (Fig 3C): Dataset1 std fitness from EXP: is LA stronger & fitness lower when neg bio interactions left intact? -------- #full model w interaction)-- summary(tpexp) #with taxonomic radom effect (doesnt change results at all) summary(tpexp[,c(1:10)]) #can do up to order fit.B100kmtx <- glmer(std_fit_site ~ local100*treat_full + (1|study) + (1|taxon) + (1|fitness_type) + (1|order/family/genus), data=tpexp, subset=localforeign100=='yes' & (treat_full=='ct' | treat_full=='positive.biotic') & exp.biotic=='yes' & study!='Abdala-RobertsandMarquis_Oecologia_2007' & study!='CarbajalNavarroetal_FrontiersEcolEvol_submitted', family="binomial", control=glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5))) fit.BnoX100kmtx <- glmer(std_fit_site ~ local100+treat_full + (1|study) + (1|taxon) + (1|fitness_type) + (1|order/family/genus), data=tpexp, subset=localforeign100=='yes' & (treat_full=='ct' | treat_full=='positive.biotic') & exp.biotic=='yes' & study!='Abdala-RobertsandMarquis_Oecologia_2007'& study!='CarbajalNavarroetal_FrontiersEcolEvol_submitted', family="binomial", control=glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5))) anova(fit.B100kmtx, fit.BnoX100kmtx) #Table 3 Q2 dataset1 std fitness: interaction p=.18 fit.BnoN100kmtx <- glmer(std_fit_site ~ treat_full + (1|study) + (1|taxon) + (1|fitness_type) + (1|order/family/genus), data=tpexp, subset=localforeign100=='yes' & (treat_full=='ct' | treat_full=='positive.biotic') & exp.biotic=='yes' & study!='Abdala-RobertsandMarquis_Oecologia_2007' & study!='CarbajalNavarroetal_FrontiersEcolEvol_submitted', family="binomial", control=glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5))) #anova(fit.B100kmtx, fit.BnoN100kmtx) #Table 3 Q2 dataset1 std fitness: LA p=.033 if testing against full model anova(fit.BnoX100kmtx, fit.BnoN100kmtx) #Table 3 Q2 dataset1 std fitness: LA p=.025 if teting against noX model #get n for fig3 caption (C) n <- tpexp[!is.na(tpexp$std_fit_site) & tpexp$localforeign100=='yes' & (tpexp$treat_full=='ct' | tpexp$treat_full=='positive.biotic') & tpexp$exp.biotic=='yes' & tpexp$study!='Abdala-RobertsandMarquis_Oecologia_2007' & tpexp$study!='CarbajalNavarroetal_FrontiersEcolEvol_submitted',] n <- droplevels(n); length(levels(n$study)) #n=15 studies dim(n) #456 data points. #__Q2 MAIN MS (Fig 3D): Dataset2 std fitness ---------- #with taxonomic random effect summary(tpCTfit[,c(1:10)]) Byn100b.tx <- glmer(std_fit_site ~ local100*unCTmanip.env + (1|study) + (1|taxon) + (1|fitness_type) + (1|phylum_division/class/order/family/genus), family="binomial", control=glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5)), data=tpCTfit, subset=(unCTmanip.env=='biotic' | unCTmanip.env=='nat') ) visreg(Byn100b.tx, xvar="local100", by="unCTmanip.env") #random effect has no effect on estimates Byn100b.noXtx <- glmer(std_fit_site ~ local100+unCTmanip.env + (1|study) + (1|taxon) + (1|fitness_type)+ (1|phylum_division/class/order/family/genus), family="binomial", control=glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5)), data=tpCTfit, subset=(unCTmanip.env=='biotic' | unCTmanip.env=='nat') ) anova(Byn100b.tx, Byn100b.noXtx, test='Chisq') #Table 3 Q2 dataset2 stdfit Chsq=13.7 lsmeans(Byn100b.tx, pairwise ~ local100|unCTmanip.env) #local > foreign in both categories #get n for Fig3 caption (D) n <- tpCTfit[!is.na(tpCTfit$std_fit_site) & (tpCTfit$unCTmanip.env=='biotic' | tpCTfit$unCTmanip.env=='nat'),] n <- droplevels(n); summary(n) nlevels(n$study) #126 studies dim(n) #6688 data points #define which model to report in results winFit.all <- Byn100b.tx #__Q2 Apendix A.2 Effect of localness of the local source? --------- #for std fitness doesnt make sense to add geo_dist as random effect - have to add for local and foreign sources so just sucks up most of the variation of interest between local and foreign sources (foreign sources are always from farther away). # #with local50 - more explicitly controls for distance of local source without confounding definition of local with fact that foreign sources always come from farther away than local ones #dataset1 --- --- fit.B50km <- glmer(std_fit_site ~ local50*treat_full + (1|study) + (1|taxon) + (1|fitness_type) + (1|order/family/genus), data=tpexp, subset=localforeign50=='yes' & (treat_full=='ct' | treat_full=='positive.biotic') & exp.biotic=='yes' & study!='Abdala-RobertsandMarquis_Oecologia_2007'& study!='CarbajalNavarroetal_FrontiersEcolEvol_submitted', family="binomial", control=glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5))) fit.BnoX50km <- glmer(std_fit_site ~ local50+treat_full + (1|study) + (1|taxon) + (1|fitness_type) + (1|order/family/genus), data=tpexp, subset=localforeign50=='yes' & (treat_full=='ct' | treat_full=='positive.biotic') & exp.biotic=='yes' & study!='Abdala-RobertsandMarquis_Oecologia_2007'& study!='CarbajalNavarroetal_FrontiersEcolEvol_submitted', family="binomial", control=glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5))) anova(fit.B50km, fit.BnoX50km) # even less signif #what is n? doesnt change n <- tpexp[!is.na(tpexp$std_fit_site) & tpexp$localforeign50=='yes' & (tpexp$treat_full=='ct' | tpexp$treat_full=='positive.biotic') & tpexp$exp.biotic=='yes' & tpexp$study!='Abdala-RobertsandMarquis_Oecologia_2007' & tpexp$study!='CarbajalNavarroetal_FrontiersEcolEvol_submitted',] n <- droplevels(n); length(levels(n$study)) #still 15 studies! just lose 3 rows #official test for local adaptation fit.BnoN50km <- glmer(std_fit_site ~ treat_full + (1|study) + (1|taxon) + (1|fitness_type) + (1|order/family/genus), data=tpexp, subset=localforeign50=='yes' & (treat_full=='ct' | treat_full=='positive.biotic') & exp.biotic=='yes' & study!='Abdala-RobertsandMarquis_Oecologia_2007' & study!='CarbajalNavarroetal_FrontiersEcolEvol_submitted', family="binomial", control=glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5))) anova(fit.BnoX50km, fit.BnoN50km) #Table 3 Q2 dataset1 std fitness: #dataset2 --- --- Byn50b <- glmer(std_fit_site ~ local50*unCTmanip.env + (1|study) + (1|taxon) + (1|fitness_type) + (1|phylum_division/class/order/family/genus), family="binomial", control=glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5)), data=tpCTfit, subset=(unCTmanip.env=='biotic' | unCTmanip.env=='nat')) #get official stats for interaction (effect of bio amel on LA) Byn50b.noX <- glmer(std_fit_site ~ local50+unCTmanip.env + (1|study) + (1|taxon) + (1|fitness_type) + (1|phylum_division/class/order/family/genus), family="binomial", control=glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5)), data=tpCTfit, subset=(unCTmanip.env=='biotic' | unCTmanip.env=='nat')) anova(Byn50b, Byn50b.noX, test='Chisq') #Table 3 Q2 dataset2 stdfit X: p was .0002, now .0011 #get official test for local adaptation lsmeans(Byn50b, pairwise ~ local50|unCTmanip.env) #local > foreign in both categories (so no change) #__Q2 SI Table A1 std fitness: results change if use only best fitness metric? ---------- #dataset 1 --- --- fit.B100kmbestf <- glmer(std_fit_site ~ local100*treat_full + (1|study) + (1|taxon) + (1|fitness_type) + (1|order/family/genus), data=tpexp, subset=localforeign100=='yes' & (treat_full=='ct' | treat_full=='positive.biotic') & as.character(fitness_type)==as.character(bestfitness), #fitness type closest to lifetime fitness family="binomial", control=glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5))) fit.B100kmbestf.noX <- glmer(std_fit_site ~ local100+treat_full + (1|study) + (1|taxon) + (1|fitness_type), data=tpexp, subset=localforeign100=='yes' & (treat_full=='ct' | treat_full=='positive.biotic') & as.character(fitness_type)==as.character(bestfitness), #only fitness type closest to lifetime fitness family="binomial", control=glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5))) anova(fit.B100kmbestf, fit.B100kmbestf.noX) #still NS #dataset 2 --- --- Byn100b.bestfit <- glmer(std_fit_site ~ local100*unCTmanip.env + (1|study) + (1|taxon) + (1|fitness_type) + (1|phylum_division/class/order/family/genus), family="binomial", data=tpCTfit, subset=unCTmanip.env!='abiotic' & unCTmanip.env!='both' & fitness_type==bestfitness, control=glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5))) Byn100b.bestfitnoX <- glmer(std_fit_site ~ local100+unCTmanip.env + (1|study) + (1|taxon) + (1|fitness_type) + (1|phylum_division/class/order/family/genus), family="binomial", data=tpCTfit, subset=unCTmanip.env!='abiotic' & unCTmanip.env!='both' & fitness_type==bestfitness, control=glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5))) anova(Byn100b.bestfit, Byn100b.bestfitnoX, test='Chisq') #Table A1 Q2 dataset2 (interaction) stdfit: p=.378 lsmeans(Byn100b.bestfit, pairwise ~ local100 | unCTmanip.env) #__Q2 SI FIG A1 E&F all categories ---------- #dataset 1 Fig A1 E ------- -- table(tpexp$genus, tpexp$treat_full) #neg.abio has 1 genus, neg.abio has 2, but genera always shared among at least 2 treatments so ok table(tpexp$family, tpexp$treat_full) table(tpexp$order, tpexp$treat_full) fitall100km.exp <- glmer(std_fit_site ~ local100*treat_full + (1|study) + (1|taxon) + (1|fitness_type) + (1|order/family/genus), data=tpexp, subset=localforeign100=='yes', family="binomial", control=glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5))) fitall100km.exp.noX <- glmer(std_fit_site ~ local100+treat_full + (1|study) + (1|taxon) + (1|fitness_type) + (1|order/family/genus), data=tpexp, subset=localforeign100=='yes', family="binomial", control=glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5))) anova(fitall100km.exp, fitall100km.exp.noX, test='Chisq') #P= 0.026 - effect of BI significant lsmeans(fitall100km.exp, pairwise ~ local100 | treat_full) #signif for neg.biotic but MALadaptation and pos.biotic (LA) #dataset 2 Fig A1 F ------- -- table(tpexp$class, tpexp$treat_full) table(tpexp$order, tpexp$treat_full) table(tpexp$family, tpexp$treat_full) fitall100km.ct <- glmer(std_fit_site ~ local100*unCTmanip.env + (1|study) + (1|taxon) + (1|fitness_type) + #runs forever but no warnings (1|phylum_division/class/order/family/genus), family="binomial", control=glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5)), data=tpCTfit) fitall100km.ctnoX <- glmer(std_fit_site ~ local100+unCTmanip.env + (1|study) + (1|taxon) + (1|fitness_type) + (1|phylum_division/class/order/family/genus), family="binomial", control=glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5)), data=tpCTfit) anova(fitall100km.ct, fitall100km.ctnoX, test='Chisq') #Fig A1.F interaction p=.00363 lsmeans(fitall100km.ct, pairwise ~ local100 | unCTmanip.env) #signif for nat (LA), biotic (LA), & both (LA), ## QUESTION 3: Dataset1 ONLY std fitness do biotic interactions affect fitness? ################## #can only do with standardized fitness (not LA metrics) and experimental manipulations (otherwise biotic treatment confounded with study) #this = reduced model from Q2 std fitness without interaction) #__Q3 MAIN MS: do biotic interactions affect fitness? (ie is treatment signif?) ###### #with taxonomic random effect summary(tpexp[,c(1:10)]) #can do up to order fit.BnoT100km.tx <- glmer(std_fit_site ~ local100 + (1|study) + (1|taxon) + (1|fitness_type) + (1|order/family/genus), data=tpexp, subset=localforeign100=='yes' & (treat_full=='ct' | treat_full=='positive.biotic') & exp.biotic=='yes' & study!='Abdala-RobertsandMarquis_Oecologia_2007' & study != "CarbajalNavarroetal_FrontiersEcolEvol_submitted", family="binomial", control=glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5))) anova(fit.BnoT100km.tx, fit.BnoX100kmtx) #Table3 Q3 dat1 effect of biotic interactions on fitness (Chsq=70.4) lsmeans(fit.BnoX100kmtx, ~ treat_full, transform='response') #Q3 text: #ct=0.498, pos.biotic=0.867 #is there an overall signal of local adaptation? fit.BnoN100km.tx <- glmer(std_fit_site ~ treat_full + (1|study) + (1|taxon) + (1|fitness_type) + (1|order/family/genus), data=tpexp, subset=localforeign100=='yes' & (treat_full=='ct' | treat_full=='positive.biotic') & exp.biotic=='yes' & study!='Abdala-RobertsandMarquis_Oecologia_2007' & study != "CarbajalNavarroetal_FrontiersEcolEvol_submitted", family="binomial", control=glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5))) anova(fit.BnoN100km.tx, fit.BnoX100kmtx) #Table3 Q3 for LA vis.mod<-visreg(fit.BnoX100km, xvar="local") vis.mod<-visreg(fit.BnoX100km, xvar="treat_full") #Effect by localness of local source? (skip adding random covariate as doesnt make sense for std fitness, see Q2) #use stricter definition of local #dataset 1 fit.BnoX50km #local + treat fit.BnoT50km <- glmer(std_fit_site ~ local50 + (1|study) + (1|taxon) + (1|fitness_type) + (1|order/family/genus), data=tpexp, subset=localforeign50=='yes' & (treat_full=='ct' | treat_full=='positive.biotic') & exp.biotic=='yes' & study!='Abdala-RobertsandMarquis_Oecologia_2007' & study != "CarbajalNavarroetal_FrontiersEcolEvol_submitted", family="binomial", control=glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5))) anova(fit.BnoT50km, fit.BnoX50km) #Table3 Q3 for effect of biotic interactions on fitness #is there an overall signal of local adaptation? fit.BnoN50km <- glmer(std_fit_site ~ treat_full + (1|study) + (1|taxon) + (1|fitness_type) + (1|order/family/genus), data=tpexp, subset=localforeign50=='yes' & (treat_full=='ct' | treat_full=='positive.biotic') & exp.biotic=='yes' & study!='Abdala-RobertsandMarquis_Oecologia_2007' & study != "CarbajalNavarroetal_FrontiersEcolEvol_submitted", family="binomial", control=glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5))) anova(fit.BnoN50km, fit.BnoX50km) #Table3 Q3 for LA #__Q3 SI TABLEA1: ct vs pos.biotic, only best fitness metric (still >1 data points per site if have >1 temporal replicates) ----- summary(tpexp) fit.B100kmbestf.noX <- glmer(std_fit_site ~ local100 + treat_full + (1|study) + (1|taxon) + (1|fitness_type) + (1|order/family/genus), data=tpexp, subset=localforeign100=='yes' & (treat_full=='ct' | treat_full=='positive.biotic') & as.character(fitness_type)==as.character(bestfitness), family="binomial", control=glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5))) fit.B100kmbestf.noT <- glmer(std_fit_site ~ local100 + (1|study) + (1|taxon) + (1|fitness_type) + (1|order/family/genus), data=tpexp, subset=localforeign100=='yes' & (treat_full=='ct' | treat_full=='positive.biotic') & as.character(fitness_type)==as.character(bestfitness), family="binomial", control=glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5))) anova(fit.B100kmbestf.noT, fit.B100kmbestf.noX) #effect of treatmnet fit.B100kmbestf.noN <- glmer(std_fit_site ~ treat_full + (1|study) + (1|taxon) + (1|fitness_type) + (1|order/family/genus), data=tpexp, subset=localforeign100=='yes' & #excl if missing local or foreign source (treat_full=='ct' | treat_full=='positive.biotic') & as.character(fitness_type)==as.character(bestfitness), family="binomial", control=glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5))) anova(fit.B100kmbestf.noN, fit.B100kmbestf.noX) #overall LA) ## QUESTION 4: cases where biotic+ treatment detects maladaptation and ct treatment detects LA? ##### #code redone (simplified) Nov 2018. #categorize results from control treatments (just do all - will filter to correct comparisons in next step) ct.full <- tpexpLA[!is.na(tpexpLA$corr_log_ratio100) & tpexpLA$treat_full=='ct',]; ct.full <- droplevels(ct.full) ct.full$LA100_3cat <- as.factor(ifelse(ct.full$LA100_cat=='yes' & ct.full$corr_log_ratio100>0, 'LA', #corr_log_ratio>0 should be redundant ifelse(ct.full$corr_log_ratio100==0, 'NoA', #no adaptation ifelse(ct.full$LA100_cat=='no' & ct.full$corr_log_ratio100<0, 'MA', NA)))) #maladaptation ct.full$LA50_3cat <- as.factor( ifelse(!is.na(ct.full$LA50_cat) & ct.full$LA50_cat=='yes' & ct.full$corr_log_ratio50>0, 'LA', ifelse(!is.na(ct.full$LA50_cat) & ct.full$corr_log_ratio50==0, 'NoA', #no adaptation ifelse(!is.na(ct.full$LA50_cat) & ct.full$LA50_cat=='no' & ct.full$corr_log_ratio50<0, 'MA', NA)))) #maladaptation summary(ct.full) #should be no NAs in LA100_3cat & 1 NA in LA50 #make streamlined dataframe and rename variables ct <- ct.full[,c('study','taxon','site','fitness_type','lifestage','temp_rep','treat_full', 'corr_log_ratio100','LA100_3cat','corr_log_ratio50','LA50_3cat')]; head(ct) colnames(ct) <- c('study','taxon','site','fitness_type','lifestage', 'temp_rep', 'cttreat_full', 'ctcorr_log_ratio100','ctLA100_3cat','ctcorr_log_ratio50','ctLA50_3cat') #categorize results from treatments that ammeliorate biotic env pb.full <- tpexpLA[!is.na(tpexpLA$corr_log_ratio100) & tpexpLA$treat_full=='positive.biotic',]; pb.full <- droplevels(pb.full) pb.full$LA100_3cat <- as.factor(ifelse(pb.full$LA100_cat=='yes' & pb.full$corr_log_ratio100>0, 'LA', #local adaptation ifelse(pb.full$corr_log_ratio100==0, 'NoA', #no adaptation ifelse(pb.full$LA100_cat=='no' & pb.full$corr_log_ratio100<0, 'MA', NA)))) #maladaptation pb.full$LA50_3cat <- as.factor( ifelse(!is.na(pb.full$LA50_cat) & pb.full$LA50_cat=='yes' & pb.full$corr_log_ratio50>0, 'LA', ifelse(!is.na(pb.full$LA50_cat) & pb.full$corr_log_ratio50==0, 'NoA', #no adaptation ifelse(!is.na(pb.full$LA50_cat) & pb.full$LA50_cat=='no' & pb.full$corr_log_ratio50<0, 'MA', NA)))) #maladaptation pb <- pb.full[,c('study', 'taxon', 'site', 'fitness_type', 'lifestage', 'temp_rep', 'treat_full', 'corr_log_ratio100','LA100_3cat','corr_log_ratio50','LA50_3cat')]; head(pb) colnames(pb) <- c('study', 'taxon', 'site', 'fitness_type', 'lifestage', 'temp_rep', 'pbtreat_full', 'pbcorr_log_ratio100','pbLA100_3cat','pbcorr_log_ratio50','pbLA50_3cat') summary(pb) #merge together so that ct and pb treatments line up dim(ct); dim(pb) #118 vs 76 rows LAvMAvNA <- merge(ct,pb, all.y=T) dim(LAvMAvNA) LAvMAvNA$cttreat_full <- LAvMAvNA$pbtreat_full <- NULL #no longer applies as both treats represented #check whether every lifestage/metric in pb had a corresponding ct treatment summary(LAvMAvNA[is.na(LAvMAvNA$ctLA100_3cat),]) # 2 pb from SambattiandRice_Evolution_2006 at missing a ct value #make column designating whether categorical conclusion about LA differed between ct and pos.biotic treatment LAvMAvNA$LA100comparison <- as.factor(ifelse(LAvMAvNA$ctLA100_3cat==LAvMAvNA$pbLA100_3cat, 'same', ifelse(LAvMAvNA$ctLA100_3cat!=LAvMAvNA$pbLA100_3cat, 'different', NA))) LAvMAvNA$LA50comparison <- as.factor(ifelse(LAvMAvNA$ctLA50_3cat==LAvMAvNA$pbLA50_3cat, 'same', ifelse(LAvMAvNA$ctLA50_3cat!=LAvMAvNA$pbLA50_3cat, 'different', NA))) summary(LAvMAvNA$LA100comparison) #22 diff, 52 same, 2 NAs summary(LAvMAvNA$LA50comparison) #20 diff, 53 same, 3 NAs #__methods ----- --- summary(LAvMAvNA) #LAcomparison: 22 diff, 52 same, 2 NAs #total of 74 cases with ct and pos.biotic treatment #__Q4: 1st: how often do ct vs bitoic+ treatments reach same vs different conclusions re LA? ----- --- summary(LAvMAvNA$LA100comparison) #changes the result 22 / (22+52) times, ie 22/(22+52) #29.7% of the time. summary(LAvMAvNA$LA50comparison) #changes the result 20 / (20+53) times, ie 20/(20+53) #27.4% of the time. #__Q4: 2nd: how often does changing BI lead to 'false maladaptation' (ie LA detected in ct, MA detected in pos.biotic treat)? vs the reverse? #original definition of local (<100km & <100 masl) LAvMAvNA$falseLAMA <- as.factor(ifelse(LAvMAvNA$ctLA100_3cat=='LA' & LAvMAvNA$pbLA100_3cat=='MA', 'falseMA', ifelse(LAvMAvNA$ctLA100_3cat=='MA' & LAvMAvNA$pbLA100_3cat=='LA', 'falseLA', NA))) summary(LAvMAvNA) #falseLA=6, falseMA=15 LAvMA <- LAvMAvNA[!is.na(LAvMAvNA$falseLAMA),]; LAvMA$pbcorr_log_ratio<-LAvMA$ctcorr_log_ratio100 <-LAvMA$LAcomparison <-LAvMA$temp_rep <-NULL; LAvMA <-droplevels(LAvMA) dim(LAvMA) #only 21 - ie lose case where one treatment was neutral LAvMA #two cases (Bromus_fasciculatus & Elymus_glaucus) where same taxon x site xtemprep are contributing 2 data points - one for comp and one for surv. if take out surv in each case have falseLA=6, falseMA=13 #false maladpataion ('success') vs false local-adaptation ('failure') #two sided binomial test: does biotic amelioration result in prob of false maladaptation diffferent from 0.5? binom.test(x=c(13,6), p=0.5, alterlocal='two.sided') #NS: P=0.167 (1 data point per taxon x site). #one sided binomial test: does biotic amelioration result in prob false maladaptation >0.5 (alternate)? binom.test(x=c(13,6), p=0.5, alterlocal='greater') #NS: P=0.0835 (1 data point per taxon x site). one-sided #which studies? 11 studies table(LAvMA$study, LAvMA$falseLAMA) #stricter definition of local (<50km & <100 masl) LAvMAvNA$falseLAMA <- as.factor(ifelse(LAvMAvNA$ctLA50_3cat=='LA' & LAvMAvNA$pbLA50_3cat=='MA', 'falseMA', ifelse(LAvMAvNA$ctLA50_3cat=='MA' & LAvMAvNA$pbLA50_3cat=='LA', 'falseLA', NA))) summary(LAvMAvNA) #falseLA=6, falseMA=15 (NO CHANGE) ### QUESTION 5: Effect of life stage (both datasets) ################################################ #only use the catLA and logratioLA. If use standardized fitness need a 3-way interaction and we dont have the statistical power #do detect more / stronger LA as get closer to LTF? #dont inlucde composite fitness metrics as confound different lifestages summary(tpexpLA) bwplot(corr_log_ratio100 ~ treat_full|fitness_type, data=tpexpLA, subset=(treat_full=='ct' | treat_full=='positive.biotic')) #__Q5: Dataset1 binary LA (Table4 row1) ---------- #with taxonomic random effect (no weights as dont work for binary data) summary(tpexpLA[,c(1:10)]) #leave out order as models stop converging (results are exactly the same as if leave order in) LAcatvsf.exptx <- glmer(LA100_cat ~ treat_full*fitness_type + (1|study) + (1|taxon) + (1|family/genus), data=tpexpLA, subset=(treat_full=='ct' | treat_full=='positive.biotic') & study!='Abdala-RobertsandMarquis_Oecologia_2007' & study!='CarbajalNavarroetal_FrontiersEcolEvol_submitted' & fitness_type!='comp' & fitness_type!='germ_surv', family=binomial, control=glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5))) LAcatvsf.expnoXtx <- glmer(LA100_cat ~ treat_full+fitness_type + (1|study) + (1|taxon) + (1|family/genus), data=tpexpLA, subset=(treat_full=='ct' | treat_full=='positive.biotic') & study!='Abdala-RobertsandMarquis_Oecologia_2007' & study!='CarbajalNavarroetal_FrontiersEcolEvol_submitted' & fitness_type!='comp' & fitness_type!='germ_surv', family=binomial, control=glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5))) #same warning anova(LAcatvsf.expnoXtx, LAcatvsf.exptx, test='Chisq') #Table 4 row1: interaction NS #treat only (treat = signif) LAcatvsf.expTtx <- glmer(LA100_cat ~ treat_full + (1|study) + (1|taxon) + (1|family/genus), data=tpexpLA, subset=(treat_full=='ct' | treat_full=='positive.biotic') & study!='Abdala-RobertsandMarquis_Oecologia_2007' & study!='CarbajalNavarroetal_FrontiersEcolEvol_submitted' & fitness_type!='comp' & fitness_type!='germ_surv', family=binomial, control=glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5))) #2 warnings #fitness component only LAcatvsf.expFtx <- glmer(LA100_cat ~ fitness_type + (1|study) + (1|taxon) + (1|family/genus), data=tpexpLA, subset=(treat_full=='ct' | treat_full=='positive.biotic') & study!='Abdala-RobertsandMarquis_Oecologia_2007' & study!='CarbajalNavarroetal_FrontiersEcolEvol_submitted' & fitness_type!='comp' & fitness_type!='germ_surv', family=binomial, control=glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5))) #1 warning anova(LAcatvsf.expnoXtx, LAcatvsf.expFtx) #Table 4 row1: effect treat compared to no interaction model lsmeans(LAcatvsf.expnoXtx, ~ treat_full) #ct > pos.biotic anova(LAcatvsf.expnoXtx, LAcatvsf.expTtx) #Table 4 row1: effect fitness comp compared to no interaction model #overall signal of LA? since treat is signif divide by treat lsmeans(LAcatvsf.expnoXtx, ~ treat_full) #no (CI onverlap 0) lsmeans(LAcatvsf.expnoX, ~ treat_full | fitness_type) #no (CI onverlap 0) #__Q5: Dataset1 effect size (log ratio) (Table 4) ---------- #with taxonomic correction random effect and weighted by sqrt(n) summary(tpexpLA[,c(1:10)]) #leave out order - makes model almost unidentifiable names(tpexpLA) table(tpexpLA$fitness_type, tpexpLA$class) #2/3 orders have >1 fitcomp and none have all 5 table(tpexpLA$fitness_type, tpexpLA$order) #6/14 orders have >1 fitcomp and none have all 5 table(tpexpLA$fitness_type, tpexpLA$family) #6/16 families have >1 fitcomp and none have all 5 nlevels(tpexpLA$family) logLAvsf.expwtx <- lmer(corr_log_ratio100 ~ treat_full*fitness_type + (1|study) + (1|taxon) + (1|family/genus), weights=sqrt(LA100_nindivsum), data=tpexpLA, subset=(treat_full=='ct' | treat_full=='positive.biotic') & study!='Abdala-RobertsandMarquis_Oecologia_2007' & study!='CarbajalNavarroetal_FrontiersEcolEvol_submitted' & fitness_type!='comp' & fitness_type!='germ_surv') logLAvsf.expnoXwtx <- lmer(corr_log_ratio100 ~ treat_full+fitness_type + (1|study) + (1|taxon) + (1|family/genus), weights=sqrt(LA100_nindivsum), data=tpexpLA, subset=(treat_full=='ct' | treat_full=='positive.biotic') & study!='Abdala-RobertsandMarquis_Oecologia_2007' & study!='CarbajalNavarroetal_FrontiersEcolEvol_submitted' & fitness_type!='comp' & fitness_type!='germ_surv') anova(logLAvsf.expnoXwtx, logLAvsf.expwtx) #Table 4 dat1 effectsize: int NS (P=.86) #look at effect of treatment (model fitness type only) logLAvsf.expFwtx <- lmer(corr_log_ratio100 ~ fitness_type + (1|study) + (1|taxon) + (1|family/genus), weights=sqrt(LA100_nindivsum), data=tpexpLA, subset=(treat_full=='ct' | treat_full=='positive.biotic') & study!='Abdala-RobertsandMarquis_Oecologia_2007' & study!='CarbajalNavarroetal_FrontiersEcolEvol_submitted' & fitness_type!='comp' & fitness_type!='germ_surv') anova(logLAvsf.expnoXwtx, logLAvsf.expFwtx) #Table 4 dat1 effectsize: treatment lsmeans(logLAvsf.expnoXtx, pairwise ~ treat_full) #look at effect of fitness type (model treatment only) logLAvsf.expTwtx <- lmer(corr_log_ratio100 ~ treat_full + (1|study) + (1|taxon) + (1|family/genus), weights=sqrt(LA100_nindivsum), data=tpexpLA, subset=(treat_full=='ct' | treat_full=='positive.biotic') & study!='Abdala-RobertsandMarquis_Oecologia_2007' & study!='CarbajalNavarroetal_FrontiersEcolEvol_submitted' & fitness_type!='comp' & fitness_type!='germ_surv') anova(logLAvsf.expnoXtx, logLAvsf.expTtx) #Table 4 dat1 effectsize: fitness component p=.158 #overall signal of LA? lsmeans(logLAvsf.expFtx, ~ fitness_type) #no (CI overlap 0) lsmeans(logLAvsf.expTtx, ~ treat_full) #no (CI overlap 0) lsmeans(logLAvsf.exptx, ~ treat_full | fitness_type) #no #__Q5: Dataset2 binary LA (Table4 row3) ---------- #with taxonomic random effect, no weighting summary(tpCTLA[,c(1:10)]) LAcatvsf.cttx <- glmer(LA100_cat ~ unCTmanip.env*fitness_type + (1|study) + (1|taxon) + (1|phylum_division/class/order/family/genus), data=tpCTLA, subset=(unCTmanip.env=='nat' | unCTmanip.env=='biotic') & #fig A1 shows all possible manipulations (fitness_type!='comp' & fitness_type!='germ_surv'), #composite metrics confound lifestages family=binomial, control=glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5))) LAcatvsf.ctnoXtx <- glmer(LA100_cat ~ unCTmanip.env + fitness_type + (1|study) + (1|taxon) + (1|phylum_division/class/order/family/genus), data=tpCTLA, subset=(unCTmanip.env=='nat' | unCTmanip.env=='biotic') & (fitness_type!='comp' & fitness_type!='germ_surv'), family=binomial, control=glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5))) anova(LAcatvsf.ctnoXtx, LAcatvsf.cttx) #Table 4 dat2 binaryLA: interaction bck to significant! P=0.037 lsmeans(LAcatvsf.cttx, pairwise ~ unCTmanip.env | fitness_type) #signif for repro (LA bio > nat) lsmeans(LAcatvsf.cttx, pairwise ~ fitness_type | unCTmanip.env) #signif for repro (LA bio > nat) #__Q5: Dataset2 LA effect size (log ratio) (Table 4 row4) ---------- #original definition of local , with taxonomic random effect & weighted by sqrt(n) summary(tpCTLA) logLAvsf.ctwtx <- lmer(corr_log_ratio100 ~ unCTmanip.env*fitness_type + (1|study) + (1|taxon)+ (1|phylum_division/class/order/family/genus), weights=sqrt(LA100_nindivsum), data=tpCTLA, subset=(unCTmanip.env=='nat' | unCTmanip.env=='biotic') & (fitness_type!='comp' & fitness_type!='germ_surv')) logLAvsf.ctnoXwtx <- lmer(corr_log_ratio100 ~ unCTmanip.env + fitness_type + (1|study) + (1|taxon) + (1|phylum_division/class/order/family/genus), weights=sqrt(LA100_nindivsum), data=tpCTLA, subset=(unCTmanip.env=='nat' | unCTmanip.env=='biotic') & (fitness_type!='comp' & fitness_type!='germ_surv')) anova(logLAvsf.ctnoXwtx, logLAvsf.ctwtx) #Table 4 dat2 last row: int NS #explore effect of manip (model fitness comp only) logLAvsf.ctFwtx <- lmer(corr_log_ratio100 ~ fitness_type + (1|study) + (1|taxon) + (1|phylum_division/class/order/family/genus), weights=sqrt(LA100_nindivsum), data=tpCTLA, subset=(unCTmanip.env=='nat' | unCTmanip.env=='biotic') & (fitness_type!='comp' & fitness_type!='germ_surv')) anova(logLAvsf.ctnoXwtx, logLAvsf.ctFwtx) #Table 4 dat2 last row: treatment NS #effect of lifestage (model manip only) logLAvsf.ctMwtx <- lmer(corr_log_ratio100 ~ unCTmanip.env + (1|study) + (1|taxon) + (1|phylum_division/class/order/family/genus), weights=sqrt(LA100_nindivsum), data=tpCTLA, subset=(unCTmanip.env=='nat' | unCTmanip.env=='biotic') & (fitness_type!='comp' & fitness_type!='germ_surv')) anova(logLAvsf.ctnoXwtx, logLAvsf.ctMwtx) #Table 4 fitness component: NS #overall LA? lsmeans(logLAvsf.ctFwtx, ~ fitness_type) #nope lsmeans(logLAvsf.ctMwtx, ~ unCTmanip.env) #nope lsmeans(logLAvsf.ctwtx, ~ unCTmanip.env | fitness_type) #nope ### QUESTION 6: Effect of latitude (dataset 2; dataset 1 not enough power) ################################################# summary(tpCTLA) summary(tpCTLA[is.na(tpCTLA$lat_site),]) #NAs from Thompsonetal_NewPhytol_1991 - done in Dee estuary UK so force to be temperate tpCTLA$abs_lat_site <- ifelse(tpCTLA$study=='Thompsonetal_NewPhytol_1991', 53.3, abs(tpCTLA$lat_site)) #absolute lat summary(tpCTLA$abs_lat_site) tpCTLA$bin_lat_site <- as.factor(ifelse(tpCTLA$abs_lat_site < 23.5, 'tropics', 'temperate')) #try categorical comparison summary(tpCTLA[,c('abs_lat_site','bin_lat_site')]) levels(tpCTLA$study) #__QLat Dataset2: Probability of LA, ---------- #again with taxonomic correction - only include levels with at least 5 levels of overlap table(tpCTLA$bin_lat_site, tpCTLA$phylum_division) #5 temp, 3trop, only 1 (vasc plants) overlap table(tpCTLA$bin_lat_site, tpCTLA$class) #2 overlap (both plants) table(tpCTLA$bin_lat_site, tpCTLA$order) #10 overlap table(tpCTLA$bin_lat_site, tpCTLA$family) #5 overlap table(tpCTLA$bin_lat_site, tpCTLA$genus) #only 1 of overlap out of 155 nlevels(tpCTLA$genus) #155 levels LAcatvslat.cttx <- glmer(LA100_cat ~ unCTmanip.env*bin_lat_site + (1|study) + (1|taxon) + (1|fitness_type) + (1|order/family), data=tpCTLA, subset=(unCTmanip.env=='nat' | unCTmanip.env=='biotic'), family=binomial, control=glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5))) LAcatvslat.ctnoXtx <- glmer(LA100_cat ~ unCTmanip.env + bin_lat_site + (1|study) + (1|taxon) + (1|fitness_type)+ (1|order/family), data=tpCTLA, subset=(unCTmanip.env=='nat' | unCTmanip.env=='biotic'), family=binomial, control=glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5))) anova(LAcatvslat.cttx, LAcatvslat.ctnoXtx, test='Chisq') #p=.0381 lsmeans(LAcatvslat.cttx, ~ unCTmanip.env | bin_lat_site) #signif LA in temperate zone (both treatments), not in tropics (neither treatment) lsmeans(LAcatvslat.cttx, pairwise ~ unCTmanip.env | bin_lat_site) #in tropics LA stronger in nat but NS p=.086 #get n n <- tpCTLA[tpCTLA$unCTmanip.env=='nat' | tpCTLA$unCTmanip.env=='biotic',]; n <- droplevels(n) table(n$unCTmanip.env, n$bin_lat_site) #full breakdown ntemp <- n[n$bin_lat_site=='temperate',]; ntemp <- droplevels(ntemp) tb <- unique(cbind(as.character(ntemp$unCTmanip.env), as.character(ntemp$study))); summary(tb) ntrop <- n[n$bin_lat_site=='tropics',]; ntrop <- droplevels(ntrop) tb <- unique(cbind(as.character(ntrop$unCTmanip.env), as.character(ntrop$study))); summary(tb) tpCTLA[tpCTLA$bin_lat_site=='tropics' & tpCTLA$unCTmanip.env=='nat',] #4 studies! woot woot #__QLat Dataset2: strength (effect size) LA --- --- #weighting by sqrt(n) with taxonomic correction logLAvslat.ctwtx <- lmer(corr_log_ratio100 ~ unCTmanip.env*bin_lat_site + (1|study) + (1|taxon) + (1|fitness_type) + (1|order/family), weights=sqrt(LA100_nindivsum), data=tpCTLA, subset=(unCTmanip.env=='nat' | unCTmanip.env=='biotic')) logLAvslat.ctnoXwtx <- lmer(corr_log_ratio100 ~ unCTmanip.env+bin_lat_site + (1|study) + (1|taxon) + (1|fitness_type) + (1|order/family), weights=sqrt(LA100_nindivsum), data=tpCTLA, subset=(unCTmanip.env=='nat' | unCTmanip.env=='biotic')) anova(logLAvslat.ctwtx, logLAvslat.ctnoXwtx) #X NS p=.51 (p=.49 w Garrido) #__QLat Dataset2: strength (std fitness) --- --- tpCTfit$abs_lat_site <- ifelse(tpCTfit$study=='Thompsonetal_NewPhytol_1991', 53.3, abs(tpCTfit$lat_site)) tpCTfit$bin_lat_site <- as.factor(ifelse(tpCTfit$abs_lat_site < 23.5, 'tropics', 'temperate')) #try categorical comparison summary(tpCTfit) #with taxonomic correction fitLAvsblat.cttx <- lmer(std_fit_site ~ unCTmanip.env*local100*bin_lat_site + (1|study) + (1|taxon) + (1|fitness_type)+ (1|order/family), data=tpCTfit, subset=(unCTmanip.env=='nat' | unCTmanip.env=='biotic')) fitLAvsblat.ctno3Xtx <- lmer(std_fit_site ~ unCTmanip.env*bin_lat_site + unCTmanip.env*local100 + bin_lat_site*local100 + (1|study) + (1|taxon) + (1|fitness_type) + (1|order/family), data=tpCTfit, subset=(unCTmanip.env=='nat' | unCTmanip.env=='biotic')) anova(fitLAvsblat.cttx, fitLAvsblat.ctno3Xtx, test='Chisq') #3way interaction NS P=.088 (P=.069 w Garrido) lsmeans(fitLAvsblat.cttx, pairwise ~ local100 | unCTmanip.env*bin_lat_site) # #Q3: is effect of ammeliorating bioX bigger in tropics? dont do as should be using dataset1 ## FIGURES ############################################################################################ #__funnel plots --------------- #where should reference line be? >0 if local adaptation is more common than maladaptation levels(tpexpLA$fitness_type) plot(1 ~ 1, pch=21, cex=1.5, bg='#440154FF', xlim=c(0,1.2), ylim=c(0,1.2)) points(.9, .9, pch=21, cex=1.5, bg='#482677FF') points(.8, .8, pch=21, cex=1.5, bg='#404788FF') points(.7, .7, pch=21, cex=1.5, bg='#33638DFF') points(.6, .6, pch=21, cex=1.5, bg='#287D8EFF') points(.5, .5, pch=21, cex=1.5, bg='#1F968BFF') points(.4, .4, pch=21, cex=1.5, bg='#29AF7FFF') points(.3, .3, pch=21, cex=1.5, bg='#3CBB75FF') points(.2, .2, pch=21, cex=1.5, bg='#55C667FF') points(.1, .1, pch=21, cex=1.5, bg='#B8DE29FF') points(.0, .0, pch=21, cex=1.5, bg='#DCE319FF') #__funnel plot to match map: sqrt n by mainvs full data, col=datset ----- quartz(width=5, height=7.5) layout(matrix(c(1,2), nrow=2, ncol=1)) par(mar=c(4,1,3,.5), oma=c(0,3,0,1.5)); xlim=c(0,150); ylim=c(-5,5) cxp=.9; cxtitl=1; ltit=-1.5 #coldat1='#DCE319FF'; coldat2='#482677FF'; bgdat1=NA; bgdat2=NA coldat1=NA; coldat2=NA; bgdat1 <- adjustcolor('#DCE319FF', alpha.f=0.8); bgdat2 <- adjustcolor('#482677FF', alpha.f=0.65); tpexpLAin <- tpexpLA[tpexpLA$treat_type=='ct' | tpexpLA$treat_type=='biotic',] plot(corr_log_ratio100 ~ sqrt(LA100_nindivsum), cex=cxp, pch=21, col=coldat2, bg=bgdat2, las=1, xaxt='n', xlab=NA, xlim=xlim, ylim=ylim, data=tpCTLA) abline(h=0) points(corr_log_ratio100 ~ sqrt(LA100_nindivsum), cex=cxp, pch=21, col=coldat1, bg=bgdat1, data=tpexpLA) #mtext(side=3, line=ltit, 'full data', cex=cxtitl) axis(side=1, at=seq(0,150,by=50), labels=NA) xd =rep(118,2); yd=c(3,4); cxlg=.9 points(x=xd, y=yd, cex=2, pch=22, col=NA, bg=c(bgdat2,bgdat1)) text(x=xd, y=yd, labels=c('147 studies', '23 studies'), cex=cxlg, pos=4) plot(corr_log_ratio100 ~ (sqrt(LA100_nindivsum)), cex=cxp, pch=21, col=coldat2, bg=bgdat2, las=1, xlab=NA, xlim=xlim, ylim=ylim, data=tpCTLAin) abline(h=0) abline(h=.091, col='#DCE319FF') #estimated corrlog ratio for dataset1 abline(h=0.067, col=bgdat2) points(corr_log_ratio100 ~ sqrt(LA100_nindivsum), cex=cxp, pch=21, col=coldat1, bg=bgdat1, data=tpexpLAin) points(x=xd, y=yd, cex=2, pch=22, col=NA, bg=c(bgdat2,bgdat1)) text(x=xd, y=yd, labels=c('126 studies', '15 studies'), cex=cxlg, pos=4) mtext(side=2, line=2.5, at=5, 'response ratio (ln(meanLocal/meanForeign))', cex=cxtitl) mtext(side=1, line=2.8, 'square root of sample size', cex=cxtitl) #__Fig 1) map - where do studies come from? ########################################################## library (maptools) library (raster) library (rgeos) #not in JLY original script but crop command failed without it #Full data set: first need dataframe with 1 point per study summary(tpexp) #has lats & longs but >1 point per site summary(tpCTfit) #has lats & longs but >1 point per site #reduce to get 1 row per study x taxon x site) map1 <- tpCTfit[tpCTfit$local100=='yes' & (tpCTfit$fitness_type==tpCTfit$bestfitness) & tpCTfit$treat_effect=='none',]; map1 <- droplevels(map1); summary(map1) #missing location info for 2 data map1[is.na(map1$lat_site),] #both are from Thompson 1991 - use rought coordinates of Dee Estuary so get a point dim(tpCTfit); dim(map1) #gets from >7000 to 765 data points nlevels(map1$study) #147 studies #boil down to 1 row per study x taxon x site (experiment, functional group, continent are just for tracking) map2 <- aggregate(map1[c('lat_site','long_site')], by=map1[c('study','functional_group','taxon','site','experiment','continent_site')], mean, na.rm=T) map2 <- droplevels(map2) colnames(map2) <- c('study','functional_group','taxon','site','experiment','continent','lat_site','long_site') summary(map2) levels(map2$study) #147 studies (kept Thompson oddly) dim(map2) #516 sites #get mean lat and mean long for each study, while counting the number of rows (site x taxa combinations) map3 <- do.call(data.frame, aggregate(map2[c('lat_site','long_site')], #by=map2[c('study','functional_group','experiment','continent')], #Hancocketal_ has 2 function groups by=map2[c('study','experiment','continent')], function(x) c(av=mean(x), nsite=length(x)))) map3$lat_site.nsite <- NULL summary(map3) colnames(map3) <- c('study','experiment','continent','lat_siteav','long_siteav','nsitetax') dim(map3) #147 data points nlevels(map3$study) #147 studies #add in rough coordinates for Thompson study (Dee Estuary according to wikipedia) map3$lat_siteav[map3$study=='Thompsonetal_NewPhytol_1991'] <- 53.3 map3$long_siteav[map3$study=='Thompsonetal_NewPhytol_1991'] <- -3.16 head(map3) summary(map3) mapdatfull <- map3 #reduce to only points included in main analyses - Dataset 2 map1b <- map1[map1$unCTmanip.env!='abiotic' & map1$unCTmanip.env!='both',]; map1b <- droplevels(map1b) nlevels(map1b$study) #126 studies #boil down to 1 row per study x taxon x site (experiment, functional group, continent are just for tracking) map2b <- aggregate(map1b[c('lat_site','long_site')], by=map1b[c('study','functional_group','taxon','site','experiment','continent_site')], mean, na.rm=T) colnames(map2b) <- c('study','functional_group','taxon','site','experiment','continent','lat_site','long_site') nlevels(map2b$study) #126 studies (kept Thompson oddly) dim(map2b) #472 sites #get mean lat and mean long for each study, while counting the number of rows (site x taxa combinations) map3b <- do.call(data.frame, aggregate(map2b[c('lat_site','long_site')], by=map2b[c('study','experiment','continent')], function(x) c(av=mean(x), nsite=length(x)))) map3b$lat_site.nsite <- NULL summary(map3b) colnames(map3b) <- c('study','experiment','continent','lat_siteav','long_siteav','nsitetax') dim(map3b) #126 data points #add in rough coordinates for Thompson study (Dee Estuary according to wikipedia) map3b$lat_siteav[map3b$study=='Thompsonetal_NewPhytol_1991'] <- 53.3 map3b$long_siteav[map3b$study=='Thompsonetal_NewPhytol_1991'] <- -3.16 head(map3b) summary(map3b) mapdat2main <- map3b #again for dataset1 summary(tpexp) map1c <- tpexp[tpexp$local100=='yes' & tpexp$treat_full=='positive.biotic',]; map1c <- droplevels(map1c) nlevels(map1c$study) #15 studies #boil down to 1 row per study x taxon x site (experiment, functional group, continent are just for tracking) map2c <- aggregate(map1c[c('lat_site','long_site')], by=map1c[c('study','functional_group','taxon','site','experiment','continent_site')], mean, na.rm=T) colnames(map2c) <- c('study','functional_group','taxon','site','experiment','continent','lat_site','long_site') nlevels(map2c$study) #15 studies (kept Thompson oddly) dim(map2c) #57 sites #get mean lat and mean long for each study, while counting the number of rows (site x taxa combinations) map3c <- do.call(data.frame, aggregate(map2c[c('lat_site','long_site')], by=map2c[c('study','experiment','continent')], function(x) c(av=mean(x), nsite=length(x)))) map3c$lat_site.nsite <- NULL summary(map3c) summary(map3c[is.na(map3c$lat_siteav),]) colnames(map3c) <- c('study','experiment','continent','lat_siteav','long_siteav','nsitetax') dim(map3c) #15 data points #add in rough coordinates for Thompson study (Dee Estuary according to wikipedia) map3c$lat_siteav[map3c$study=='Thompsonetal_NewPhytol_1991'] <- 53.3 map3c$long_siteav[map3c$study=='Thompsonetal_NewPhytol_1991'] <- -3.16 head(map3c) summary(map3c) mapdat1main <- map3c summary(mapdatfull) summary(mapdat1main) summary(mapdat2main) # Get map shapefile world <- readShapePoly("ne_110m_land/ne_110m_land.shp") # Crop to study area (do later in ppt or word) dev.off() plot(world, border="white",col="grey") abline(h=0) #equator # get transparent colours coldat2 <- adjustcolor('#482677FF', alpha.f=0.7) coldat1 <- adjustcolor('white', alpha.f=0.9) #higher 3=less transparent coldat1 <- adjustcolor('#DCE319FF', alpha.f=0.9) coltrop <- adjustcolor('grey', alpha.f=0.3) #Map BioRXiv full dataset: 1 point per study with size indicating n taxa and sites quartz(width=9, height=4) #shape does not adjust to quartz dimensions, just get lots fo white space par(mar=c(0,0,0,0)) plot(world,border="white",col="grey"); #plot equator & tropics abline(h=0, lwd=.8); #equator rect(xleft=-180, xright=180, ytop=23.5, ybot=-23.5, border=NA, col=coltrop) #tropics #plot 1 point per study, size representing 3 of sites + taxa in study #studies in dataset2 only points(mapdatfull[mapdatfull$experiment=='no', c("long_siteav","lat_siteav")], cex=mapdatfull$nsitetax[mapdatfull$experiment=='no']*.15+.4, pch=21, bg=coldat2) #studies in dataset1(&2) points(mapdatfull[mapdatfull$experiment=='yes', c("long_siteav","lat_siteav")], cex=mapdatfull$nsitetax[mapdatfull$experiment=='yes']*.15+.4, pch=21, bg=coldat1) #legend legx=-148; xsz=rep(legx,3); ysz=y=c(0,-8.5,-19)-33 points(x=xsz, y=c(ysz[1],ysz[2]+.5,ysz[3]), cex=c(1,8,16)*.15+.4, pch=21, col="black") text(x=xsz+10, y=ysz, labels=c('1','8','16'), cex=.8) text(x=legx, y=ysz[1]+7, labels='sites x taxa', cex=.8) xd = rep(legx,2); yd=c(ysz[2],ysz[3])-1; xd=rep(-118,2); yd=ysz[2:3] points(x=xd, y=yd, cex=2, pch=22, col=NA, bg=c(coldat1,coldat2)) text(x=xd+2, y=yd-1, labels=c('dataset 1', 'dataset 2'), cex=.8, pos=4) #Map added for revised version showing only studies in main analyses (ie nat/ct vs biotically ammerliorated) quartz(width=9, height=4) par(mar=c(0,0,0,0)) plot(world,border="white",col="grey"); #plot equator & tropics abline(h=0, lwd=.8); #equator rect(xleft=-180, xright=180, ytop=23.5, ybot=-23.5, border=NA, col=coltrop) #tropics #plot 1 point per study, size representing # of sites x taxa in study #studies in dataset2 only points(mapdat2main[mapdat2main$experiment=='no', c("long_siteav","lat_siteav")], cex=mapdat2main$nsitetax[mapdat2main$experiment=='no']*.15+.4, pch=21, bg=coldat2) points(mapdat1main[,c("long_siteav","lat_siteav")], cex=mapdat1main$nsitetax*.15+.4, pch=21, bg=coldat1) mapdat1main[,c("long_siteav","lat_siteav")] mapdat2main[mapdat2main$experiment=='no', c("long_siteav","lat_siteav")] #legend legx=-148; xsz=rep(legx,3); ysz=y=c(0,-8.5,-19)-33 points(x=xsz, y=c(ysz[1],ysz[2]+.5,ysz[3]), cex=c(1,8,16)*.15+.4, pch=21, col="black") text(x=xsz+10, y=ysz, labels=c('1','8','16'), cex=.8) text(x=legx, y=ysz[1]+7, labels='sites x taxa', cex=.8) xd = rep(legx,2); yd=c(ysz[2],ysz[3])-1; xd=rep(-118,2); yd=ysz[2:3] points(x=xd, y=yd, cex=2, pch=22, col=NA, bg=c(coldat1,coldat2)) text(x=xd+2, y=yd-1, labels=c('dataset 1', 'dataset 2'), cex=.8, pos=4) #__Fig 2) Is LA more common w nat biotic interactions? ################################################# #panel 1 data - experiments (define which model want to use for plotting) winLAcat.exp <- LAcatpB100.tx #get backtransformed means & confidence intervals winLAcat.explsm <- data.frame(summary(lsmeans(winLAcat.exp, ~ treat_full))); winLAcat.explsm winLAcat.explsm$btmean <- exp(winLAcat.explsm$lsmean)/(1+exp(winLAcat.explsm$lsmean)) winLAcat.explsm$btLCL <- exp(winLAcat.explsm$asymp.LCL)/(1+exp(winLAcat.explsm$asymp.LCL)) winLAcat.explsm$btUCL <- exp(winLAcat.explsm$asymp.UCL)/(1+exp(winLAcat.explsm$asymp.UCL)) winLAcat.explsm #get backtransformed residuals vcatLAexp <- visreg(winLAcat.exp, xvar='treat_full', partial=T, scale='response'); vcatLAexp #backtransformed scale head(vcatLAexp$fit) #fitted means head(vcatLAexp$res) #residuals summary(vcatLAexp$res) vcatLAexp$res$treat_fulln <- ifelse(vcatLAexp$res$treat_full=='ct',1,2) #hack to get numeric xaxis so can jitter #panel 2 data - experiments (define which model want to use for plotting) winLAcat.all <- LActcatB100.tx Anova(winLAcat.all) #get backtransformed means & confidence intervals winLAcat.alllsm <- data.frame(summary(lsmeans(winLAcat.all, ~ unCTmanip.env))); winLAcat.alllsm winLAcat.alllsm$btmean <- exp(winLAcat.alllsm$lsmean)/(1+exp(winLAcat.alllsm$lsmean)) winLAcat.alllsm$btLCL <- exp(winLAcat.alllsm$asymp.LCL)/(1+exp(winLAcat.alllsm$asymp.LCL)) winLAcat.alllsm$btUCL <- exp(winLAcat.alllsm$asymp.UCL)/(1+exp(winLAcat.alllsm$asymp.UCL)) #get backtransformed residuals #vcatLAall <- visreg(winLAcat.all, xvar='unCTmanip', partial=T) #transformed scale, asking for residuals vcatLAall <- visreg(winLAcat.all, xvar='unCTmanip.env', partial=T, scale='response') #backtransformed w residuals head(vcatLAall$fit) #fitted means vcatLAall$fit$unCTmanipn <- factor(vcatLAall$fit$unCTmanip.env, levels=levels(vcatLAall$fit$unCTmanip.env)[c(2,1)]) #hack to get numeric xaxis so can re-order with nat first head(vcatLAall$res) #residuals vcatLAall$res$unCTmanipn <- ifelse(vcatLAall$res$unCTmanip.env=='nat',1,2) #hack to get numeric xaxis so can jitter #Fig 2. 2 col x 1 rows quartz(width=5, height=2.9) layout(matrix(c(1,2), nrow=1, ncol=2, byrow=F)) mar=c(1.2,1,1,.2); #xlim=c(0,.6); c(seq(0,.6,by=.1)) ymin=0; ymax=1; ref=.5; ystat=0; yline=2.7; ytck=seq(ymin,ymax,length.out=6); #backtransformed scale ylab='Probability of local adaptation'; colline=c(1,.3)+.3 ylimLA=c(ymin,ymax); ytck par(oma=c(1,3,1.3,0), bty='l') #mfcol draws plots by column cxp=.55; abw=1.5; xtck=c(seq(0,1)); #ytck=c(seq(0,60,by=10)) #if want to plot points have to extend x axis xl=2.9; xtl=.55 #line for xax label & xax tick labels cxaxt=.85; cxsig=.8 #size axis tick labels & test giving model significance colCT='grey'; colpB='grey'; colCT='chartreuse4'; colpB='darkorange2'; colCTCI <- adjustcolor(colCT, alpha.f=.6); colpBCI <- adjustcolor(colpB, alpha.f=.6) #colors for ct/nat vs biotic amelioration (pos.bio) # par(mar=mar) #Fig 2 panel A (dataset1) --- --- plot(vcatLAexp$fit$visregFit ~ vcatLAexp$fit$treat_full, ylim=ylimLA, las=1, yaxt='n', xaxt='n') axis(side=2, at=ytck, labels=NA, las=1) mtext(side=2, at=ytck, text=ytck, cex=.9, las=1, line=1) abline(h=ref, col='grey', lwd=abw) xct=1; r=.4; rect(xleft=xct-r, xright=xct+r, ytop=winLAcat.explsm$btUCL[1], ybottom =winLAcat.explsm$btLCL[1], col=colCTCI, border=NA) #backtrnsfrm rect(xleft=2-r, xright=2+r, ytop=winLAcat.explsm$btUCL[2], ybottom =winLAcat.explsm$btLCL[2], col=colpBCI, border=NA) #backtransfrmed points(jitter(x=vcatLAexp$res$treat_fulln[vcatLAexp$res$treat_full=='ct'], factor=18), y=vcatLAexp$res$visregRes[vcatLAexp$res$treat_full=='ct'], cex=cxp, pch=16, col=colCTCI) points(jitter(x=vcatLAexp$res$treat_fulln[vcatLAexp$res$treat_full=='positive.biotic'], factor=9), y=vcatLAexp$res$visregRes[vcatLAexp$res$treat_full=='positive.biotic'], cex=cxp, pch=16, col=colpBCI) par(new=T); plot(vcatLAexp$fit$visregFit ~ vcatLAexp$fit$treat_full, ylim=ylimLA, xaxt='n', yaxt='n') #replots lines incase dots obscured them mtext(side=3, line=colline, text=c('controlled manipulations','within studies'), cex=.9) axis(side=1, at=c(1,2), labels=NA) mtext(side=1, at=c(1,2), line=xtl, text=c('control', 'biotic'), cex=.9) mtext(side=1, at=c(1,2), line=xtl+.7, text=c('', 'amelioration'), cex=.9) mtext(side=2, line=yline, text=ylab, cex=.9) text(x=.4, y=ystat, labels='biotic manipulation NS', pos=4, cex=cxsig) mtext(at=.6, side=3, line=-.7, text='A', cex=1.2) # #Fig 2 panel B (dataset2) --- --- plot(vcatLAall$fit$visregFit ~ vcatLAall$fit$unCTmanipn, ylim=ylimLA, las=1, yaxt='n', xaxt='n') axis(side=2, at=ytck, labels=NA) abline(h=ref, col='grey', lwd=abw) rect(xleft=1-r, xright=1+r, ytop=winLAcat.alllsm$btUCL[2], ybottom=winLAcat.alllsm$btLCL[2], col=colCTCI, border=NA) rect(xleft=2-r, xright=2+r, ytop=winLAcat.alllsm$btUCL[1], ybottom=winLAcat.alllsm$btLCL[1], col=colpBCI, border=NA) points(jitter(x=vcatLAall$res$unCTmanipn[vcatLAall$res$unCTmanip.env=='nat'], factor=18), y=vcatLAall$res$visregRes[vcatLAall$res$unCTmanip.env=='nat'], cex=cxp, pch=16, col=colCTCI) points(jitter(x=vcatLAall$res$unCTmanipn[vcatLAall$res$unCTmanip.env=='biotic'], factor=9), y=vcatLAall$res$visregRes[vcatLAall$res$unCTmanip.env=='biotic'], cex=cxp, pch=16, col=colpBCI) par(new=T); plot(vcatLAall$fit$visregFit ~ vcatLAall$fit$unCTmanipn, ylim=ylimLA, yaxt='n', xaxt='n') mtext(side=3, line=colline, text=c('uncontrolled manipulations','among studies'), cex=.9) axis(side=1, at=c(1,2), labels=NA) mtext(side=1, at=c(1,2), line=xtl, text=c('natural', 'biotic'), cex=.9) mtext(side=1, at=c(1,2), line=xtl+.7, text=c('', 'amelioration'), cex=.9) text(x=.4, y=ystat, labels='biotic manipulation NS', pos=4, cex=cxsig) text(x=1.5, y=.94, labels='*', cex=2) text(x=c(1,2), y=.94, labels=c('*','*'), cex=2, col=c(colCT,colpB)) mtext(at=.6, side=3, line=-.7, text='B', cex=1.2) #__Fig 3) Is LA stronger w nat biotic interactions? ##################################### #data prep --- --- # (A) top row: LA (log diff) / left col: experiments winLA.exp <- LAposB100.wsqtx winLA.explsm <- data.frame(summary(lsmeans(winLA.exp, ~ treat_full))); winLA.explsm visreg(winLA.exp, xvar='treat_full') #need y axis limits to 4.5 vLAexp <- visreg(winLA.exp, xvar='treat_full', partial=T, plot=F) #transformed scale, asking for residuals vLAexp$res$treat_fulln <- ifelse(vLAexp$res$treat_full=='ct',1,2) #hack to get numeric xaxis so can jitter # (B) top row: LA (log diff) / right col: all studies par(mfrow=c(1,1), mar=c(4,3,1,1)) winLA.all <- LAB100ctc.wsqtx; winLA.all winLA.alllsm <- data.frame(summary(lsmeans(winLA.all, ~ unCTmanip.env))); winLA.alllsm visreg(winLA.all, xvar='unCTmanip.env', partial=T) #need y axis to 6 vLAall <- visreg(winLA.all, xvar='unCTmanip.env', partial=T, plot=F) #transformed scale, askign for residuals head(vLAall$fit) #fitted means (nat should be 1st treatment) vLAall$fit$unCTmanipn <- as.factor(ifelse(vLAall$fit$unCTmanip.env=='nat','nat','xbiotic')) #hack to make sure order right head(vLAall$fit) #fitted means head(vLAall$res) #residuals vLAall$res$unCTmanipn <- ifelse(vLAall$res$unCTmanip.env=='nat',1,2) #hack to get numeric xaxis so can jitter # (C) bottom row: fitness / left col: experiments ### winFit.exp <- fit.B100kmtx; winFit.exp winFit.explsm <- data.frame(summary(lsmeans(winFit.exp, ~ local100| treat_full))); winFit.explsm winFit.explsm$btmean <- exp(winFit.explsm$lsmean)/(1+exp(winFit.explsm$lsmean)) winFit.explsm$btLCL <- exp(winFit.explsm$asymp.LCL)/(1+exp(winFit.explsm$asymp.LCL)) winFit.explsm$btUCL <- exp(winFit.explsm$asymp.UCL)/(1+exp(winFit.explsm$asymp.UCL)) winFit.explsm summary(tpexp$std_fit_site) #bounded between 0 and 1, median is .56 Anova(winFit.exp) vFitexp <- visreg(winFit.exp, xvar='local100', by="treat_full", partial=T, layout=c(2,1)) vFitexp <- visreg(winFit.exp, xvar='local100', by="treat_full", partial=T, scale='response', layout=c(2,1)) vFitexp$fit$localn <- factor(ifelse(vFitexp$fit$local100=='yes', 1, 2)); vFitexp$fit #hack to get numeric adj=.1 vFitexp$fit$order <- ifelse(vFitexp$fit$treat_full=='ct' & vFitexp$fit$local100=='yes', 1-adj, ifelse(vFitexp$fit$treat_full=='ct' & vFitexp$fit$local100=='no', 2-adj, ifelse(vFitexp$fit$treat_full=='positive.biotic' & vFitexp$fit$local100=='yes', 3+adj, 4+adj))) head(vFitexp$fit) #fitted means head(vFitexp$res) #residuals vFitexp$res$treat_fulln <- ifelse(vFitexp$res$treat_full=='ct',1,2) #hack to get numeric xaxis so can jitter vFitexp$res$localn <- ifelse(vFitexp$res$local100=='yes',1,2) vFitexp$res$order <- ifelse(vFitexp$res$treat_full=='ct' & vFitexp$res$local100=='yes', 1-adj, ifelse(vFitexp$res$treat_full=='ct' & vFitexp$res$local100=='no', 2-adj, ifelse(vFitexp$res$treat_full=='positive.biotic' & vFitexp$res$local=='yes', 3+adj, 4+adj))) # (D) bottom row: fitness / right col: ct all studies winFit.all <- Byn100b.tx winFit.alllsm <- data.frame(summary(lsmeans(winFit.all, ~ local100 | unCTmanip.env))); winFit.alllsm winFit.alllsm$btmean <- exp(winFit.alllsm$lsmean)/(1+exp(winFit.alllsm$lsmean)) winFit.alllsm$btLCL <- exp(winFit.alllsm$asymp.LCL)/(1+exp(winFit.alllsm$asymp.LCL)) winFit.alllsm$btUCL <- exp(winFit.alllsm$asymp.UCL)/(1+exp(winFit.alllsm$asymp.UCL)) winFit.alllsm visreg(winFit.all, xvar="local100",by="unCTmanip.env", layout=c(2,1)) vFitall <- visreg(winFit.all, xvar='local100', by="unCTmanip.env", partial=T) vFitall <- visreg(winFit.all, xvar='local100', by="unCTmanip.env", partial=T, scale='response') vFitall$fit$order <- ifelse(vFitall$fit$unCTmanip.env=='nat' & vFitall$fit$local100=='yes', 1-adj, ifelse(vFitall$fit$unCTmanip.env=='nat' & vFitall$fit$local100=='no', 2-adj, ifelse(vFitall$fit$unCTmanip=='biotic' & vFitall$fit$local100=='yes', 3+adj, 4+adj))) head(vFitall$fit) #fitted means head(vFitall$res) #residuals vFitall$res$order <- ifelse(vFitall$res$unCTmanip.env=='nat' & vFitall$res$local100=='yes', 1-adj, ifelse(vFitall$res$unCTmanip.env=='nat' & vFitall$res$local100=='no', 2-adj, ifelse(vFitall$res$unCTmanip.env=='biotic' & vFitall$res$local100=='yes', 3+adj, 4+adj))) #__plotting guidelines --- --- quartz(width=5, height=5.5) layout(matrix(c(1,2,3,4), nrow=2, ncol=2, byrow=T)) mar=c(1.5,1,1,.5); par(oma=c(2,3.5,1.5,0), bty='l') #mfcol draws plots by column cxp=.5; xtck=c(seq(0,1)); ytck=c(seq(0,60,by=10)) #if want to plot points have to extend x axis abw=1.5; xl=1.8; yl=2.6 cxaxt=.85; cxsig=.9 #size axis tick labels colCT='grey'; colpB='grey'; colCT='chartreuse4'; colpB='darkorange2'; colCTCI <- adjustcolor(colCT, alpha.f=.6); colpBCI <- adjustcolor(colpB, alpha.f=.6) xct=1; r=.4 #defines width of rectangle # (A) ---- - par(mar=mar); ymin=-5.3; ymax=6; ylimLA=c(ymin, ymax); plot(vLAexp$fit$visregFit ~ vLAexp$fit$treat_full, ylim=ylimLA, xaxt='n',las=1) abline(h=0, col='grey', lwd=abw) rect(xleft=xct-r, xright=xct+r, ytop=winLA.explsm$lower.CL[1], ybottom =winLA.explsm$upper.CL[1], col=colCTCI, border=NA) rect(xleft=2-r, xright=2+r, ytop=winLA.explsm$lower.CL[2], ybottom =winLA.explsm$upper.CL[2], col=colpBCI, border=NA) points(jitter(x=vLAexp$res$treat_fulln[vLAexp$res$treat_full=='ct'], factor=18), y=vLAexp$res$visregRes[vLAexp$res$treat_full=='ct'], cex=cxp, pch=16, col=colCTCI) points(jitter(x=vLAexp$res$treat_fulln[vLAexp$res$treat_full=='positive.biotic'], factor=9), y=vLAexp$res$visregRes[vLAexp$res$treat_full=='positive.biotic'], cex=cxp, pch=16, col=colpBCI) par(new=T); plot(vLAexp$fit$visregFit ~ vLAexp$fit$treat_full, ylim=ylimLA, xaxt='n', yaxt='n') axis(side=1, at=c(1,2), labels=NA) mtext(side=1, at=c(1,2), line=xtl, text=c('control', 'biotic ameliorat.'), cex=.9) mtext(side=3, line=colline, text=c('controlled manipulations','within studies'), cex=.9) mtext(side=2, line=yl, text='Effect size local adaptation', cex=.9) text(x=.4, y=ymin+.1, labels='biotic manipulation NS', pos=4, cex=cxsig) mtext(at=.6, side=3, line=-.8, text='A', cex=1.2) # (B) get more points w this code than original visreg - havent figured out why. doesnt happen with panel A --- #WATCH OUT! even if treatments are properly ordered (nat first) visreg plots them alphabetically, ie biotic then natrual plot(vLAall$fit$visregFit ~ vLAall$fit$unCTmanipn, ylim=ylimLA, las=1, yaxt='n', xaxt='n') axis(side=2, at=c(-4,-2,0,2,4,6), labels=NA) axis(side=1, at=c(1,2), labels=NA) mtext(side=1, at=c(1,2), line=xtl, text=c('natural', 'biotic ameliorat.'), cex=.9) abline(h=0, col='grey', lwd=abw) points(jitter(x=vLAall$res$unCTmanipn[vLAall$res$unCTmanip.env=='nat'], factor=18), y=vLAall$res$visregRes[vLAall$res$unCTmanip.env=='nat'], cex=cxp, pch=16, col=colCTCI) points(jitter(x=vLAall$res$unCTmanipn[vLAall$res$unCTmanip.env=='biotic'], factor=9), y=vLAall$res$visregRes[vLAall$res$unCTmanip.env=='biotic'], cex=cxp, pch=16, col=colpBCI) rect(xleft=xct-r, xright=xct+r, ytop=winLA.alllsm$lower.CL[winLA.alllsm$unCTmanip.env=='nat'], ybottom=winLA.alllsm$upper.CL[winLA.alllsm$unCTmanip.env=='nat'], col=colCTCI, border=NA) rect(xleft=2-r, xright=2+r, ytop=winLA.alllsm$lower.CL[winLA.alllsm$unCTmanip.env=='biotic'], ybottom=winLA.alllsm$upper.CL[winLA.alllsm$unCTmanip.env=='biotic'], col=colpBCI, border=NA) par(new=T); plot(vLAall$fit$visregFit ~ vLAall$fit$unCTmanipn, ylim=ylimLA, las=1, yaxt='n', xaxt='n') mtext(side=3, line=colline, text=c('uncontrolled manipulations','among studies'), cex=.9) text(x=.4, y=ymin+.1, labels='biotic manipulation P = 0.043', pos=4, cex=cxsig) mtext(at=.6, side=3, line=-.8, text='B', cex=1.2) # (C) --- - winFit.explsm cxp; cxp=.4; xlim=c(0.4,4.5); x=c(1-adj,2-adj,3+adj,4+adj); xtxt=c('local','foreign','local','foreign') ymin=-1.9; ylimFit=c(ymin,3.6); stdfytck=c(-2,-1,0,1,2,3,4) #for logit scale ymin=0; ylimFit=c(ymin,1); stdfytck=seq(0,1,length.out=6) #for transformed scale plot(vFitexp$fit$visregFit ~ vFitexp$fit$order, xlim=xlim, ylim=ylimFit, las=1, pch='—', cex=2.5, xaxt='n') axis(side=1, at=x, labels=NA) mtext(side=1, line=xtl-.2, at=x, text=xtxt, cex=.8) # ct local points(jitter(x=rep(1-adj,length=length(vFitexp$res$visregRes[vFitexp$res$treat_full=='ct' & vFitexp$res$local100=='yes'])), factor=20), y=vFitexp$res$visregRes[vFitexp$res$treat_full=='ct' & vFitexp$res$local100=='yes'], cex=cxp, pch=16, col=colCTCI) rect(xleft=1-adj-r, xright=1-adj+r, ytop=winFit.explsm$btUCL[1], ybottom =winFit.explsm$btLCL[1], col=colCTCI, border=NA) # ct foreign points(jitter(x=rep(2-adj, length=length(vFitexp$res$visregRes[vFitexp$res$treat_full=='ct' & vFitexp$res$local100=='no'])), factor=9), y=vFitexp$res$visregRes[vFitexp$res$treat_full=='ct' & vFitexp$res$local100=='no'], cex=cxp, pch=16, col=colCTCI) rect(xleft=2-adj-r, xright=2-adj+r, ytop=winFit.explsm$btUCL[2], ybottom =winFit.explsm$btLCL[2], col=colCTCI, border=NA) # biotic+ local points(jitter(x=rep(3+adj,length=length(vFitexp$res$visregRes[vFitexp$res$treat_full=='positive.biotic' & vFitexp$res$local100=='yes'])), factor=6), y=vFitexp$res$visregRes[vFitexp$res$treat_full=='positive.biotic' & vFitexp$res$local100=='yes'], cex=cxp, pch=16, col=colpBCI) rect(xleft=3+adj-r, xright=3+adj+r, ytop=winFit.explsm$btUCL[3], ybottom =winFit.explsm$btLCL[3], col=colpBCI, border=NA) # biotic+ foreign points(jitter(x=rep(4+adj,length=length(vFitexp$res$visregRes[vFitexp$res$treat_full=='positive.biotic' & vFitexp$res$local100=='no'])), factor=4.), y=vFitexp$res$visregRes[vFitexp$res$treat_full=='positive.biotic' & vFitexp$res$local100=='no'], cex=cxp, pch=16, col=colpBCI) rect(xleft=4+adj-r, xright=4+adj+r, ytop=winFit.explsm$btUCL[4], ybottom =winFit.explsm$btLCL[4], col=colpBCI, border=NA) par(new=T); plot(vFitexp$fit$visregFit ~ vFitexp$fit$order, xlim=xlim, ylim=ylimFit, xaxt='n', yaxt='n', pch='—', cex=2.5) #redo mean lines mtext(side=2, line=yl, text='Standardized fitness', cex=.9) mtext(at=.6, side=3, line=-.9, text='C', cex=1.2) xat=c(1.5-adj,3.5+adj+adj); mtext(at=xat, side=1, line=xl, text=c('control','biotic ameliorat.'), cex=.9) text(x=.2, y=ymin, labels='source x biotic manipulation = NS', pos=4, cex=cxsig) text(x=.2, y=ymin+.08, labels='biotic manipulation P < 0.0001', pos=4, cex=cxsig) #text(x=2.48, y=.99, labels='*', cex=2.4) #signif across treatments as treat x source NS text(x=c(1.4,3.65), y=1, labels=c('*','*'), cex=2.4, col=c(colCT,colpB)) #signif bt treatments # (D) winFit.alllsm #treatments now in same order as figure plot(vFitall$fit$visregFit ~ vFitall$fit$order, xlim=xlim, ylim=ylimFit, las=1, pch='—', cex=2.5, xaxt='n', yaxt='n') axis(side=1, at=x, labels=NA) axis(side=2, at=stdfytck, labels=NA) mtext(side=1, line=xtl-.2, at=x, text=xtxt, cex=.8) # ct local points(jitter(x=rep(1-adj,length=length(vFitall$res$visregRes[vFitall$res$unCTmanip.env=='nat' & vFitall$res$local100=='yes'])), factor=19), y=vFitall$res$visregRes[vFitall$res$unCTmanip.env=='nat' & vFitall$res$local100=='yes'], cex=cxp, pch=16, col=colCTCI) rect(xleft=1-adj-r, xright=1-adj+r, ytop=winFit.alllsm$btUCL[1], ybottom =winFit.alllsm$btLCL[1], col=colCTCI, border=NA) # ct foreign points(jitter(x=rep(2-adj, length=length(vFitall$res$visregRes[vFitall$res$unCTmanip.env=='nat' & vFitall$res$local100=='no'])), factor=9), y=vFitall$res$visregRes[vFitall$res$unCTmanip.env=='nat' & vFitall$res$local100=='no'], cex=cxp, pch=16, col=colCTCI) rect(xleft=2-adj-r, xright=2-adj+r, ytop=winFit.alllsm$btUCL[2], ybottom =winFit.alllsm$btLCL[2], col=colCTCI, border=NA) # biotic+ local points(jitter(x=rep(3+adj,length=length(vFitall$res$visregRes[vFitall$res$unCTmanip.env=='biotic' & vFitall$res$local=='yes'])), factor=5), y=vFitall$res$visregRes[vFitall$res$unCTmanip.env=='biotic' & vFitall$res$local=='yes'], cex=cxp, pch=16, col=colpBCI) rect(xleft=3+adj-r, xright=3+adj+r, ytop=winFit.alllsm$btUCL[3], ybottom =winFit.alllsm$btLCL[3], col=colpBCI, border=NA) # biotic+ foreign points(jitter(x=rep(4+adj,length=length(vFitall$res$visregRes[vFitall$res$unCTmanip.env=='biotic' & vFitall$res$local=='no'])), factor=4), y=vFitall$res$visregRes[vFitall$res$unCTmanip.env=='biotic' & vFitall$res$local=='no'], cex=cxp, pch=16, col=colpBCI) rect(xleft=4+adj-r, xright=4+adj+r, ytop=winFit.alllsm$btUCL[4], ybottom =winFit.alllsm$btLCL[4], col=colpBCI, border=NA) par(new=T); plot(vFitall$fit$visregFit ~ vFitall$fit$order, xlim=xlim, ylim=ylimFit, xaxt='n', yaxt='n', pch='—', cex=2.5) text(x=.2, y=ymin, labels='source x biotic manip. P < 0.001', pos=4, cex=cxsig) mtext(at=.6, side=3, line=-.9, text='D', cex=1.2) mtext(at=xat, side=1, line=xl, text=c('natural','biotic ameliorat.'), cex=.9) text(x=xat, y=c(rep(1,2)), labels=c('*','*'), col=c(colCT,colpB), cex=2.4) #__Fig 4) Effect of latitude on biotic interactions (standard fitness) ---------- #make new data frame for plotting raw data points - just to shorten code tropsf <- tpCTfit[tpCTfit$bin_lat_site=='tropics',]; tropsf <- droplevels(tropsf); levels(tropsf$study) tempsf <- tpCTfit[tpCTfit$bin_lat_site=='temperate',] #Prob LA: extract means from full model LAcatvslat.cttx winLAcat.lat <- LAcatvslat.cttx #get backtransformed means & confidence intervals winpLALatlsm <- data.frame(summary(lsmeans(LAcatvslat.cttx, ~ unCTmanip.env*bin_lat_site, emm_options(pbkrtest.limit = 6580)))); winpLALatlsm winpLALatlsm$btmean <- exp(winpLALatlsm$lsmean)/(1+exp(winpLALatlsm$lsmean)) winpLALatlsm$btLCL <- exp(winpLALatlsm$asymp.LCL)/(1+exp(winpLALatlsm$asymp.LCL)) winpLALatlsm$btUCL <- exp(winpLALatlsm$asymp.UCL)/(1+exp(winpLALatlsm$asymp.UCL)) winpLALatlsm #get backtransformed residuals vcatLA <- visreg(winLAcat.lat, xvar='unCTmanip.env', by="bin_lat_site", partial=T, scale='response'); vcatLA #backtransformed scale head(vcatLA$fit) #fitted means head(vcatLA$res) #residuals summary(vcatLA$res) adj=0.14 vcatLA$res$xn <- ifelse(vcatLA$res$bin_lat_site=='temperate' & vcatLA$res$unCTmanip.env=='nat', 1-adj, ifelse(vcatLA$res$bin_lat_site=='temperate' & vcatLA$res$unCTmanip.env=='biotic', 2-adj, ifelse(vcatLA$res$bin_lat_site=='tropics' & vcatLA$res$unCTmanip.env=='nat', 3+adj, 4+adj))) #hack to get numeric xaxis to jitter vcatLA$fit$order <- ifelse(vcatLA$fit$bin_lat_site=='temperate' & vcatLA$fit$unCTmanip.env=='nat', 1-adj, ifelse(vcatLA$fit$bin_lat_site=='temperate' & vcatLA$fit$unCTmanip.env=='biotic', 2-adj, ifelse(vcatLA$fit$bin_lat_site=='tropics' & vcatLA$fit$unCTmanip.env=='nat', 3+adj, 4+adj))) #Fig 4 effect of lat on frequency of LA temperate vs tropics --- --- quartz(width=4.6, height=3) mar=c(1.2,1,1,1.5); z=.3; xlim=c(0+z,5-z) ymin=0; ymax=1.1; ylimLA=c(ymin,ymax); ref=.5; ystat=0; yline=2.7; ytck=seq(ymin,ymax,length.out=6); ylab='Probability of local adaptation'; cxp=.55; abw=1.5; xtck=c(seq(0,1)); #ytck=c(seq(0,60,by=10)) #if want to plot points have to extend x axis xl=2.9; xtl=.55 #line for xax label & xax tick labels cxt=.95; cxsig=.8 #size axis tick labels & test giving model significance colCT='chartreuse4'; colpB='darkorange2'; colCTCI <- adjustcolor(colCT, alpha.f=.6); colpBCI <- adjustcolor(colpB, alpha.f=.6) # par(oma=c(1,3,.5,0), bty='l', mar=mar) plot(vcatLA$fit$visregFit ~ vcatLA$fit$order, xlim=xlim, ylim=ylimLA, las=1, pch='—', cex=2.5, xaxt='n') abline(h=ref, col='grey', lwd=abw) xct=1-adj; r=.37; rect(xleft=xct-r, xright=xct+r, ytop=winpLALatlsm$btUCL[1], ybottom =winpLALatlsm$btLCL[1], col=colCTCI, border=NA) #backtrnsfrm rect(xleft=xct+1-r, xright=xct+1+r, ytop=winpLALatlsm$btUCL[2], ybottom =winpLALatlsm$btLCL[2], col=colpBCI, border=NA) #backtransfrmed rect(xleft=3+adj-r, xright=3+adj+r, ytop=winpLALatlsm$btUCL[3], ybottom =winpLALatlsm$btLCL[3], col=colCTCI, border=NA) #backtrnsfrm rect(xleft=4+adj-r, xright=4+adj+r, ytop=winpLALatlsm$btUCL[4], ybottom =winpLALatlsm$btLCL[4], col=colpBCI, border=NA) #backtransfrmed f=.8; points(jitter(x=vcatLA$res$xn[vcatLA$res$unCTmanip.env=='nat'], factor=f), y=vcatLA$res$visregRes[vcatLA$res$unCTmanip.env=='nat'], cex=cxp, pch=16, col=colCTCI) points(jitter(x=vcatLA$res$xn[vcatLA$res$unCTmanip.env=='biotic'], factor=f), y=vcatLA$res$visregRes[vcatLA$res$unCTmanip.env=='biotic'], cex=cxp, pch=16, col=colpBCI) par(new=T); plot(vcatLA$fit$visregFit ~ vcatLA$fit$order, xlim=xlim, ylim=ylimLA, xaxt='n', yaxt='n', pch='—', cex=3.2) #replots lines text(c(1.5-adj,3.5+adj), y=rep(1.1,2), labels=c('temperate','tropical'), cex=1) x=c(1-adj,2-adj,3+adj,4+adj) axis(side=1, at=x, labels=NA) mtext(side=1, at=x, line=xtl, text=c('natural', 'biotic'), cex=cxt) mtext(side=1, at=x, line=xtl+.7, text=c('', 'amelioration'), cex=cxt) mtext(side=2, line=yline, text=ylab, cex=1) text(x=.2, y=ystat, labels='latitude x biotic amelioration P = 0.038', pos=4, cex=cxt-.1) text(x=x, y=c(rep(.98,4)), labels=c('*','*','',''), col=c(colCT,colpB), cex=2.4) #LA only signif in temperate zone #_SI Fig A1) Freq & strength of LA in multiple enviornments ################################# # could go back ito original std fit database creation and order treat_full -> would prob simplify code # NOTE - CODE IS TRICKY!!!! sometimes order of means gets mixed up between visreg and lsmeans (eg visreg$fit and lsmeans have same numebrs but assigned to diff treatments). has happened twice, no idea why, resolves by clearing cache and running whole code again, not sure whats happening. but double check comparing lsmeans and visreg means for each plot #data prep #PLOT w EVERYTHING 2 col x 1 rows quartz(width=6.2, height=6.6) layout(matrix(c(1,2,3,4,5,6), nrow=3, ncol=2, byrow=T), widths=c(1,.65), heights=c(1,1,1)) marL=c(1,1,.5,1); marR=c(1,.5,.5,1); par(oma=c(2.2,3,1.5,0), bty='l'); jit=1.5 cxp=.45; cxbar=3; cxylab=.8; cxxlab=cxylab; cxaxt=.85; cxsig=1; #size axis tick labels & test giving model significance abw=1.5; xtck=c(seq(0,1)); ytck=c(seq(0,60,by=10)) #if want to plot points have to extend x axis titl=.6; llet=-1; lylab=2.3 #line for xax label & xax tick labels & title colCT='chartreuse4'; colpB='darkorange2'; #colors for ct vs biotic amelioration colCTCI <- adjustcolor(colCT, alpha.f=.6); colpBCI <- adjustcolor(colpB, alpha.f=.6) colgreyCI <- adjustcolor('grey', alpha.f=.5); #other manip CI grey colgreyp <- adjustcolor('grey50', alpha.f=.7) #other manip point grey colgreystar='grey50' gap=.5; xct=1; xbp=2+gap; xap=3+gap; xbothp=4+gap; xbn=5+2*gap; xan=6+2*gap; atxLA=c(xct,xbp,xap,xbn,xan) #no longer have a 'both' treatment atxexpAll=c(xct,xbp,xap,xbothp,xbn,xan) #__A probability LA dataset1 (backtransform to 0-1 scale) ----- winLAcat.exp <- LAcatall100.exp lsmeans(winLAcat.exp, ~1) #overall no LA (overlaps 0) #get backtransformed means & confidence intervals (should already be in order) winLAcat.explsm <- data.frame(summary(lsmeans(winLAcat.exp, ~ treat_full))); winLAcat.explsm winLAcat.explsm$btmean <- exp(winLAcat.explsm$lsmean)/(1+exp(winLAcat.explsm$lsmean)) winLAcat.explsm$btLCL <- exp(winLAcat.explsm$asymp.LCL)/(1+exp(winLAcat.explsm$asymp.LCL)) winLAcat.explsm$btUCL <- exp(winLAcat.explsm$asymp.UCL)/(1+exp(winLAcat.explsm$asymp.UCL)) winLAcat.explsm #btmean all 0 to 1 (proportion). ct > posbio > neg.abio > pos.abio > neg.bio #get backtransformed residuals #quartz(width=7, height=4); visreg(winLAcat.exp, xvar='treat_full', partial=T, scale='response') #see what should look like vcatLAexp <- visreg(winLAcat.exp, xvar='treat_full', partial=T, scale='response', plot=F); #backtransformed scale head(vcatLAexp$fit) #fitted means head(vcatLAexp$res) #residuals summary(vcatLAexp$res) #visregRes between 0 and 1 vcatLAexp$res$treat_fulln <- ifelse(vcatLAexp$res$treat_full=='ct', xct, ifelse(vcatLAexp$res$treat_full=='positive.biotic', xbp, ifelse(vcatLAexp$res$treat_full=='positive.abiotic', xap, ifelse(vcatLAexp$res$treat_full=='negative.biotic',xbn, xan)))) #hack to get num xaxis to jittr summary(vcatLAexp$res$treat_fulln[vcatLAexp$res$treat_full=='positive.abiotic']) #A (top left panel) prob LA dataset1 ------ -- vcatLAexp$fit winLAcat.explsm par(mar=marL); gap=.5; ymin=0; ylimcLA=c(ymin,1); ab=.5; lsig=ymin+.01; xmin=.5; xlimL=c(xmin,7.5); #order goal: ct - pos.bio pos.abio pos.both - neg.bio neg.abio plot(vcatLAexp$fit$visregFit ~ atxLA, xlim=xlimL, ylim=ylimcLA, las=1, pch='—', cex=2.5, xaxt='n') abline(h=.5, col='grey', lwd=abw) #axis(side=1, at=atxL.AC, labels=NA) axis(side=1, at=c(xct,xbp,xap,xbothp,xbn,xan), labels=NA); xtl=-.3; #tick for 'both' r=.39; rect(xleft=atxLA-r, xright=atxLA+r, ytop=winLAcat.explsm$btUCL, ybottom=winLAcat.explsm$btLCL, col=c(colCTCI,colpBCI,colgreyCI,colgreyCI,colgreyCI), border=NA) #ct (green) points (from visreg) points(jitter(x=vcatLAexp$res$treat_fulln[vcatLAexp$res$treat_full=='ct'], factor=17), y=vcatLAexp$res$visregRes[vcatLAexp$res$treat_full=='ct'], cex=cxp, pch=16, col=colCTCI) #pos.biotic (orange) points points(jitter(x=vcatLAexp$res$treat_fulln[vcatLAexp$res$treat_full=='positive.biotic'], factor=7), y=vcatLAexp$res$visregRes[vcatLAexp$res$treat_full=='positive.biotic'], cex=cxp, pch=16, col=colpBCI) #pos.abiotic (grey) points points(jitter(x=vcatLAexp$res$treat_fulln[vcatLAexp$res$treat_full=='positive.abiotic'], factor=jit+.3), y=vcatLAexp$res$visregRes[vcatLAexp$res$treat_full=='positive.abiotic'], cex=cxp, pch=16, col=colgreyp) #neg.biotic (grey) points points(jitter(x=vcatLAexp$res$treat_fulln[vcatLAexp$res$treat_full=='negative.biotic'], factor=jit+.3), y=vcatLAexp$res$visregRes[vcatLAexp$res$treat_full=='negative.biotic'], cex=cxp, pch=16, col=colgreyp) #neg.abiotic (grey) points points(jitter(x=vcatLAexp$res$treat_fulln[vcatLAexp$res$treat_full=='negative.abiotic'], factor=jit+.3), y=vcatLAexp$res$visregRes[vcatLAexp$res$treat_full=='negative.abiotic'], cex=cxp, pch=16, col=colgreyp) par(new=T); plot(vcatLAexp$fit$visregFit ~ atxLA, xlim=xlimL, ylim=ylimcLA, pch='—', cex=cxbar, yaxt='n', xaxt='n') mtext(side=3, line=titl, text='controlled treatments', cex=cxylab) mtext(side=2, line=lylab, text='Probability of local adaptation', cex=cxylab) text(x=xmin-.2, y=lsig, labels='treatment P = 0.37', pos=4, cex=cxsig) mtext(at=xmin, side=3, line=llet, text='A', cex=1.2) abline(v=c( (xct+xbp)/2, (xbothp+xbn)/2), col='black', lwd=abw/2) #if doing by effect then type #__B probability LA dataset2 (backtransform to 0-1 scale) ---------- winLAcat.all <- LAcatall100.ct lsmeans(winLAcat.all, ~1) #overall LA (>0) #get backtransformed means & confidence intervals winLAcat.alllsm <- data.frame(summary(lsmeans(winLAcat.all, ~ unCTmanip.env))); winLAcat.alllsm winLAcat.alllsm$btmean <- exp(winLAcat.alllsm$lsmean)/(1+exp(winLAcat.alllsm$lsmean)) winLAcat.alllsm$btLCL <- exp(winLAcat.alllsm$asymp.LCL)/(1+exp(winLAcat.alllsm$asymp.LCL)) winLAcat.alllsm$btUCL <- exp(winLAcat.alllsm$asymp.UCL)/(1+exp(winLAcat.alllsm$asymp.UCL)) winLAcat.alllsm #get backtransformed residuals #quartz(width=6, height=4); visreg(winLAcat.all, xvar='unCTmanip.env', partial=T, scale='response') #abio > both > bio > nat vcatLAall <- visreg(winLAcat.all, xvar='unCTmanip.env', partial=T, scale='response', plot=F) #backtransformed w residuals vcatLAall$fit #fitted means head(vcatLAall$res) #residuals vcatLAall$res$unCTmanipn <- ifelse(vcatLAall$res$unCTmanip.env=='nat',1, ifelse(vcatLAall$res$unCTmanip.env=='biotic',2, #hack to get numeric xaxis so can jitter ifelse(vcatLAall$res$unCTmanip.env=='abiotic',3,4))) #B (top right panel) prob LA ----- --- vcatLAall$fit winLAcat.alllsm xmin=.5; xlimR=c(xmin,5); atxR=c(xct,xbp,xap,xbothp) par(mar=marR) plot(vcatLAall$fit$visregFit ~ atxR, xlim=xlimR, ylim=ylimcLA, las=1, pch='—', cex=cxbar, yaxt='n', xaxt='n') axis(side=2, at=seq(0,1,by=.2), labels=NA) abline(h=ab, col='grey', lwd=abw) rect(xleft=atxR-r, xright=atxR+r, ytop=winLAcat.alllsm$btUCL, ybottom =winLAcat.alllsm$btLCL, col=c(colCTCI, colpBCI, colgreyCI, colgreyCI), border=NA) #green (nat) points points(jitter(x=vcatLAall$res$unCTmanipn[vcatLAall$res$unCTmanip.env=='nat'], factor=17), y=vcatLAall$res$visregRes[vcatLAall$res$unCTmanip.env=='nat'], cex=cxp, pch=16, col=colCTCI) points(jitter(x=vcatLAall$res$unCTmanipn[vcatLAall$res$unCTmanip.env=='biotic']+gap, factor=7), y=vcatLAall$res$visregRes[vcatLAall$res$unCTmanip.env=='biotic'], cex=cxp, pch=16, col=colpBCI) points(jitter(x=vcatLAall$res$unCTmanipn[vcatLAall$res$unCTmanipn >2.5]+ gap, factor=jit+.3), y=vcatLAall$res$visregRes[vcatLAall$res$unCTmanipn >2.5], cex=cxp, pch=16, col=colgreyp) par(new=T); plot(vcatLAall$fit$visregFit ~ atxR, xlim=xlimR, ylim=ylimcLA, pch='—', cex=cxbar, yaxt='n', xaxt='n') mtext(side=3, line=titl, text='uncontrolled alterations', cex=cxylab) axis(side=1, at=atxR, labels=NA) text(x=xmin-.2, y=lsig, labels='alteration P = 0.62', pos=4, cex=cxsig) ##Chisq=1.8, P=0.62 mtext(at=xmin+.1, side=3, line=llet, text='B', cex=1.2) abline(v=(1+2+gap)/2, col='black', lwd=abw/2) text(x=((xbp+xap)/2)-.08, y=.99, labels='*', cex=2.4) #signif across manipulations as unCTmanip NS #__C (mid left) effect size LA (dont need to backtransform) -------- winLA.exp <- LAall100exp lsmeans(winLA.exp, ~1) #overall no LA winLA.explsm <- data.frame(summary(lsmeans(winLA.exp, ~ treat_full))) winLA.explsm #residuals from visreg #quartz(width=6, height=5); visreg(LAall100exp, xvar='treat_full', partial=T) vLAexp <- visreg(LAall100exp, xvar='treat_full', partial=T, plot=F) vLAexp$res$treat_fulln <- ifelse(vLAexp$res$treat_full=='ct',1,2) #hack to get numeric xaxis so can jitter residuals head(vLAexp$res) #residuals vLAexp$res$treat_fulln <- ifelse(vLAexp$res$treat_full=='ct',xct, #hack to get numeric xaxis so can jitter ifelse(vLAexp$res$treat_full=='positive.biotic',xbp, ifelse(vLAexp$res$treat_full=='negative.biotic',xbn, ifelse(vLAexp$res$treat_full=='positive.abiotic',xap, ifelse(vLAexp$res$treat_full=='negative.abiotic',xan, xbothp))))) summary(vLAexp$res) vLAexp$fit winLA.explsm #order pos.abio > neg.bio > ct> neg.abio > pos.bio. all >0 except pos.bio par(mar=marL); ymin=-5; ymax=6; ylimLA=c(ymin, ymax); lsig=ymin plot(vLAexp$fit$visregFit ~ atxLA, xlim=xlimL, ylim=ylimLA, las=1, pch='—', cex=2.5, xaxt='n') abline(h=0, col='grey', lwd=abw) #axis(side=1, at=atxLA, labels=NA) axis(side=1, at=c(xct,xbp,xap,xbothp,xbn,xan), labels=NA); xtl=-.3; #tick for 'both' rect(xleft=atxLA-r, xright=atxLA+r, ytop=winLA.explsm$upper.CL, ybottom =winLA.explsm$lower.CL, col=c(colCTCI,colpBCI,colgreyCI,colgreyCI,colgreyCI), border=NA) points(jitter(x=vLAexp$res$treat_fulln[vLAexp$res$treat_full=='ct'], factor=16), y=vLAexp$res$visregRes[vLAexp$res$treat_full=='ct'], cex=cxp, pch=16, col=colCTCI) points(jitter(x=vLAexp$res$treat_fulln[vLAexp$res$treat_full=='positive.biotic'], factor=6.5), y=vLAexp$res$visregRes[vLAexp$res$treat_full=='positive.biotic'], cex=cxp, pch=16, col=colpBCI) points(jitter(x=vLAexp$res$treat_fulln[vLAexp$res$treat_fulln>2.5], factor=jit+.3), y=vLAexp$res$visregRes[vLAexp$res$treat_fulln>2.5], cex=cxp, pch=16, col=colgreyp) par(new=T); plot(vLAexp$fit$visregFit ~ atxLA, xlim=xlimL, ylim=ylimLA, pch='—', cex=cxbar, xaxt='n', yaxt='n') mtext(side=2, line=lylab, text='Effect size local adaptation', cex=cxylab) text(x=xmin-.2, y=lsig, labels='treatment P = 0.83', pos=4, cex=cxsig) ##Chisq=1.5, P=.83 mtext(at=xmin, side=3, line=llet, text='C', cex=1.2) abline(v=c( (xct+xbp)/2, (xbothp+xbn)/2), col='black', lwd=abw/2) #if doing by effect then type #__D (mid right) effect size LA: all studies ---------- winLA.all <- LAall100ctc; winLA.all #lsmeans(winLA.all, ~1) #no overall LA (overlaps 0) - takes a while to run w taxonomic effect winLA.alllsm <- data.frame(summary(lsmeans(winLA.all, ~ unCTmanip.env))); winLA.alllsm #dont have to baktransform vLAall <- visreg(winLA.all, xvar='unCTmanip.env', partial=T, plot=F) head(vLAall$fit) #fitted means head(vLAall$res) #residuals vLAall$res$unCTmanipn <- ifelse(vLAall$res$unCTmanip.env=='nat',1, ifelse(vLAall$res$unCTmanip.env=='biotic',2+gap, ifelse(vLAall$res$unCTmanip.env=='abiotic',3+gap,4+gap))) #hack to get numeric xaxis so can jitter winLA.alllsm #only nat is <0. rest are v close, both > abio > bio, all CI overlap 0 par(mar=marR); plot(vLAall$fit$visregFit ~ atxR, xlim=xlimR, ylim=ylimLA, las=1, pch='—', cex=2.5, yaxt='n', xaxt='n') axis(side=2, at=seq(-4,ymax,by=2), labels=NA) abline(h=0, col='grey', lwd=abw) rect(xleft=atxR-r, xright=atxR+r, ytop=winLA.alllsm$upper.CL, ybottom =winLA.alllsm$lower.CL, col=c(colCTCI, colpBCI, colgreyCI, colgreyCI), border=NA) points(jitter(x=vLAall$res$unCTmanipn[vLAall$res$unCTmanip.env=='nat'], factor=16), y=vLAall$res$visregRes[vLAall$res$unCTmanip.env=='nat'], cex=cxp, pch=16, col=colCTCI) points(jitter(x=vLAall$res$unCTmanipn[vLAall$res$unCTmanip.env=='biotic'], factor=7), y=vLAall$res$visregRes[vLAall$res$unCTmanip.env=='biotic'], cex=cxp, pch=16, col=colpBCI) points(jitter(x=vLAall$res$unCTmanipn[vLAall$res$unCTmanipn > 2.5], factor=jit+.25), y=vLAall$res$visregRes[vLAall$res$unCTmanipn > 2.5], cex=cxp, pch=16, col=colgreyp) par(new=T); plot(vLAall$fit$visregFit ~ atxR, xlim=xlimR, ylim=ylimLA, pch='—', cex=cxbar, yaxt='n', xaxt='n') axis(side=1, at=atxR, labels=NA) text(x=xmin-.2, y=lsig, labels='alteration P = 0.19', pos=4, cex=cxsig) ##Chisq=4.8, P=0.19 mtext(at=xmin+.1, side=3, line=llet, text='D', cex=1.2) abline(v=(1+2+gap)/2, col='black', lwd=abw/2) #__E (bottom left) std_fit LA: experiments (backtransform) ---------- winFit.exp <- fitall100km.exp #using all possible data Anova(winFit.exp) #interaction signif so cant get overall LA. lsmeans(winFit.exp, pairwise ~ local100 | treat_full) #signif for pos.biotic, negative.biotic winFit.explsm <- data.frame(summary(lsmeans(winFit.exp, ~ local100 | treat_full))); winFit.explsm$btmean <- exp(winFit.explsm$lsmean)/(1+exp(winFit.explsm$lsmean)) winFit.explsm$btLCL <- exp(winFit.explsm$asymp.LCL)/(1+exp(winFit.explsm$asymp.LCL)) winFit.explsm$btUCL <- exp(winFit.explsm$asymp.UCL)/(1+exp(winFit.explsm$asymp.UCL)) winFit.explsm #residuals vFitexp <- visreg(winFit.exp, xvar='local100', by="treat_full", scale='response', partial=T, plot=F) #vFitexp$fit$localn <- factor(vFitexp$fit$local, levels=levels(vFitexp$fit$local)[c(2,1)]) #hack to order adj=.2; xctL=xct-adj; xctF=xct+adj; xbpL=xbp-adj; xbpF=xbp+adj; xapL=xap-adj; xapF=xap+adj; xbothpL=xbothp-adj; xbothpF=xbothp+adj; xbnL=xbn-adj; xbnF=xbn+adj; xanL=xan-adj; xanF=xan+adj vFitexp$res$localn <- ifelse(vFitexp$res$local100=='yes',1,2) vFitexp$res$treatn <- ifelse(vFitexp$res$treat_full=='ct' & vFitexp$res$local100=='yes', xctL, ifelse(vFitexp$res$treat_full=='ct' & vFitexp$res$local100=='no', xctF, ifelse(vFitexp$res$treat_full=='positive.biotic' & vFitexp$res$local100=='yes', xbpL, ifelse(vFitexp$res$treat_full=='positive.biotic' & vFitexp$res$local100=='no', xbpF, ifelse(vFitexp$res$treat_full=='positive.abiotic' & vFitexp$res$local100=='yes', xapL, ifelse(vFitexp$res$treat_full=='positive.abiotic' & vFitexp$res$local100=='no', xapF, ifelse(vFitexp$res$treat_full=='positive.both' & vFitexp$res$local100=='yes', xbothpL, ifelse(vFitexp$res$treat_full=='positive.both' & vFitexp$res$local100=='no', xbothpF, ifelse(vFitexp$res$treat_full=='negative.biotic' & vFitexp$res$local100=='yes', xbnL, ifelse(vFitexp$res$treat_full=='negative.biotic' & vFitexp$res$local100=='no', xbnF, ifelse(vFitexp$res$treat_full=='negative.abiotic' & vFitexp$res$local100=='yes', xanL, xanF))))))))))) vFitexp$fit winFit.explsm par(mar=marL) xl=2.9; xtl=.7; cxpf=cxp-.12; cxbar=1.3 #ymin=-2.5; ylimfit=c(ymin,4.6); lsig=ymin; #on transformed scale atxLf=c(xctL, xctF, xbpL, xbpF, xapL, xapF, xbnL, xbnF, xan-adj, xan+adj) plot(vFitexp$fit$visregFit ~ atxLf, xlim=xlimL, ylim=ylimcLA, las=1, pch='—', cex=cxbar, xaxt='n') # one by one: control local (do separate to ensure plotting local vs foreign correctly) points(jitter(x=vFitexp$res$treatn[vFitexp$res$treat_full=='ct'], factor=2), y=vFitexp$res$visregRes[vFitexp$res$treat_full=='ct'], cex=cxpf, pch=16, col=colCTCI) r=.17; rect(xleft=xctL-r, xright=xctL+r, ytop=winFit.explsm$btUCL[1], ybottom =winFit.explsm$btLCL[1], col=colCTCI, border=NA) # ct foreign rect(xleft=xctF-r, xright=xctF+r, ytop=winFit.explsm$btUCL[2], ybottom =winFit.explsm$btLCL[2], col=colCTCI, border=NA) # biotic+ local & foreign (think have sorting figured out now so plot together) points(jitter(x=vFitexp$res$treatn[vFitexp$res$treat_full=='positive.biotic'], factor=2), y=vFitexp$res$visregRes[vFitexp$res$treat_full=='positive.biotic'], cex=cxpf, pch=16, col=colpBCI) rect(xleft=c(xbpL,xbpF)-r, xright=c(xbpL,xbpF)+r, ytop=winFit.explsm$btUCL[3:4], ybottom=winFit.explsm$btLCL[3:4], col=colpBCI, border=NA) #rest (grey) points(jitter(x=vFitexp$res$treatn[vFitexp$res$treatn > xbpF], factor=1), y=vFitexp$res$visregRes[vFitexp$res$treatn > xbpF], cex=cxpf, pch=16, col=colgreyp) rect(xleft=atxLf[5:10]-r, xright=atxLf[5:10]+r, ytop=winFit.explsm$btUCL[c(5:10)], ybottom=winFit.explsm$btLCL[c(5:10)], col=colgreyCI, border=NA) par(new=T); plot(vFitexp$fit$visregFit ~ atxLf, xlim=xlimL, ylim=ylimcLA, pch='—', cex=cxbar, xaxt='n', yaxt='n') mtext(side=2, line=lylab, text='Standardized fitness', cex=cxylab) text(x=xmin-.2, y=0, labels='source x treatment P = 0.026', pos=4, cex=cxsig) #P= 0.026 mtext(at=xmin, side=3, line=llet, text='E', cex=1.2) abline(v=c( (xct+xbp)/2, (xbothp+xbn)/2), col='black', lwd=abw/2) #if doing by effect then type #axis(side=1, at=atxLA, labels=NA); xtl=-.3; #no ick for 'both' axis(side=1, at=c(xct,xbp,xap,xbothp,xbn,xan), labels=NA); xtl=-.3; #tick for 'both' atxL=c(xct,xbp,xap,xbothp,xbn,xan) mtext(side=1, at=atxL, line=xtl+1, text=c('control','biotic','abiotic','both','abiotic','biotic'), cex=cxxlab) mtext(side=1, at=c(xct,3.5,6.5), line=xtl+2.3, text=c('','increase fitness','decrease fitness'), cex=cxxlab) #text(x=xbothp+.3, y=.99, labels='*', cex=2.4) #signif across treatments if treat x source NS text(x=c(xbp,xbn), y=.998, labels='*', cex=2.4, col=c(colpB,colgreystar)) #signif across treatments if treat x source NS #__F (bottom right) std_fit LA --------- #bottom row: fitness / right col: ct all studies winFit.all <- fitall100km.ct summary(tpCTfit$std_fit_site) Anova(winFit.all) #interaction signif (if use all data points) winFit.alllsm <- data.frame(summary(lsmeans(winFit.all, ~ local100 | unCTmanip.env))); winFit.alllsm$btmean <- exp(winFit.alllsm$lsmean)/(1+exp(winFit.alllsm$lsmean)) winFit.alllsm$btLCL <- exp(winFit.alllsm$asymp.LCL)/(1+exp(winFit.alllsm$asymp.LCL)) winFit.alllsm$btUCL <- exp(winFit.alllsm$asymp.UCL)/(1+exp(winFit.alllsm$asymp.UCL)) winFit.alllsm #visreg(winFit.all, xvar="local",by="unCTmanip") vFitall <- visreg(winFit.all, xvar='local100', by="unCTmanip.env", partial=T, scale='response', plot=F) head(vFitall$res) #residuals vFitall$res$treatn <- ifelse(vFitall$res$unCTmanip.env=='nat' & vFitall$res$local100=='yes', xctL, ifelse(vFitall$res$unCTmanip.env=='nat' & vFitall$res$local100=='no', xctF, ifelse(vFitall$res$unCTmanip.env=='biotic' & vFitall$res$local100=='yes', xbpL, ifelse(vFitall$res$unCTmanip.env=='biotic' & vFitall$res$local100=='no', xbpF, ifelse(vFitall$res$unCTmanip.env=='abiotic' & vFitall$res$local100=='yes', xap-adj, ifelse(vFitall$res$unCTmanip.env=='abiotic' & vFitall$res$local100=='no', xap+adj, ifelse(vFitall$res$unCTmanip.env=='both' & vFitall$res$local100=='yes', xbothp-adj, xbothp+adj))))))) winFit.alllsm vFitall$fit par(mar=marR); ymin=0 xfitR <- c(xctL,xctF,xbpL,xbpF,xap-adj,xap+adj,xbothp-adj,xbothp+adj) plot(vFitall$fit$visregFit ~ xfitR, xlim=xlimR, ylim=ylimcLA, las=1, pch='—', cex=cxbar, yaxt='n', xaxt='n') axis(side=2, at=seq(0,1,by=.2), labels=NA) # nat local points(jitter(x=vFitall$res$treatn[vFitall$res$unCTmanip.env=='nat' & vFitall$res$local100=='yes'], factor=8.5), y=vFitall$res$visregRes[vFitall$res$unCTmanip.env=='nat'& vFitall$res$local100=='yes'], cex=cxpf, pch=16, col=colCTCI) rect(xleft=xctL-r, xright=xctL+r, ytop=winFit.alllsm$btUCL[1], ybottom =winFit.alllsm$btLCL[1], col=colCTCI, border=NA) # nat foreign points(jitter(x=vFitall$res$treatn[vFitall$res$unCTmanip.env=='nat' & vFitall$res$local100=='no'], factor=6.1), y=vFitall$res$visregRes[vFitall$res$unCTmanip.env=='nat'& vFitall$res$local100=='no'], cex=cxpf, pch=16, col=colCTCI) rect(xleft=xctF-r, xright=xctF+r, ytop=winFit.alllsm$btUCL[2], ybottom =winFit.alllsm$btLCL[2], col=colCTCI, border=NA) # biotic local points(jitter(x=vFitall$res$treatn[vFitall$res$unCTmanip.env=='biotic' & vFitall$res$local100=='yes'], factor=3.4), y=vFitall$res$visregRes[vFitall$res$unCTmanip.env=='biotic'& vFitall$res$local100=='yes'], cex=cxpf, pch=16, col=colpBCI) rect(xleft=xbpL-r, xright=xbpL+r, ytop=winFit.alllsm$btUCL[3], ybottom =winFit.alllsm$btLCL[3], col=colpBCI, border=NA) # biotic foreign points(jitter(x=vFitall$res$treatn[vFitall$res$unCTmanip.env=='biotic' & vFitall$res$local100=='no'], factor=2.5), y=vFitall$res$visregRes[vFitall$res$unCTmanip.env=='biotic'& vFitall$res$local100=='no'], cex=cxpf, pch=16, col=colpBCI) rect(xleft=xbpF-r, xright=xbpF+r, ytop=winFit.alllsm$btUCL[4], ybottom =winFit.alllsm$btLCL[4], col=colpBCI, border=NA) # rest local points(jitter(x=vFitall$res$treatn[(vFitall$res$unCTmanip=='both'|vFitall$res$unCTmanip=='abiotic')], factor=1.8), y=vFitall$res$visregRes[(vFitall$res$unCTmanip=='both'|vFitall$res$unCTmanip=='abiotic')], cex=cxpf, pch=16, col=colgreyp) xrest <- c(xapL,xapF,xbothpL,xbothpF) xfitR; winFit.alllsm rect(xleft=xfitR[5:8]-r, xright=xfitR[5:8]+r, ytop=winFit.alllsm$btUCL[5:8], ybottom =winFit.alllsm$btLCL[5:8], col=colgreyCI, border=NA) par(new=T); plot(vFitall$fit$visregFit ~ xfitR, xlim=xlimR, ylim=ylimcLA, pch='—', cex=cxbar, yaxt='n', xaxt='n') text(x=xmin-.2, y=ymin, labels='source x alteration P = 0.004', pos=4, cex=cxsig) #interaction p=.00363 text(x=c(xct,xbp,xbothp)+.02, y=rep(1,3), labels=rep('*',3), cex=2.4, col=c(colCT,colpB,colgreyp)) mtext(at=xmin+.1, side=3, line=llet, text='F', cex=1.2) abline(v=(1+2+gap)/2, col='black', lwd=abw/2) axis(side=1, at=atxR, labels=NA) #if dont want line for local/foreign mtext(side=1, at=atxR, line=xtl+1, text=c('natural','biotic','abiotic','both'), cex=cxxlab) mtext(side=1, at=c(3+gap), line=xtl+2.4, text=c('increase fitness'), cex=cxxlab)