lizzieinvancouver
commented
1 year ago And here's the code, just so I have a record:
> options(mc.cores = parallel::detectCores())
>
> #'###############################
> # Flags for how to run the code #
> #'###############################
> runmodels <- FALSE
> runbbstanleadin <- T # leave as false to speed up Supp and ms. compilation
> runwithchillports <- T
>
> #'######################################
> #### get data through bbstanleadin ####
> #'######################################
>
> if(runbbstanleadin){
+ # Flags to choose for bbstanleadin.R #
+ setwd("..//bb_analysis")
+ getwd()
+ # Master flags! Here you pick if you want the flags for the main model (figure in main text) versus the all spp model (supp)
+ use.flags.for.mainmodel <- FALSE
+ use.flags.for.allsppmodel <- TRUE
+ use.yourown.flagdesign <- FALSE
+ nocrops <- TRUE
+ agiosponly <- TRUE
+
+ if(use.flags.for.mainmodel==TRUE & use.flags.for.allsppmodel | use.flags.for.mainmodel==TRUE & use.yourown.flagdesign |
+ use.yourown.flagdesign & use.flags.for.allsppmodel | use.flags.for.mainmodel==TRUE & use.flags.for.allsppmodel
+ & use.yourown.flagdesign) print("ALERT! You have set too many master flags to true, you must pick only one!")
+
+ if(use.flags.for.mainmodel){
+ use.chillports = FALSE
+ use.zscore = TRUE
+ use.allspp =FALSE # for the main model this is false
+ use.multcuespp = FALSE
+ use.cropspp = FALSE
+ # Default is species complex use alltypes of designs
+ use.expramptypes.fp = TRUE
+ use.exptypes.fp = FALSE
+ use.expchillonly = FALSE
+ }
+
+ if(use.flags.for.allsppmodel){
+ use.chillports = FALSE
+ use.zscore = TRUE
+ use.allspp = TRUE
+ use.multcuespp = FALSE
+ use.cropspp = TRUE
+ use.expramptypes.fp = FALSE
+ use.exptypes.fp = FALSE
+ use.expchillonly = FALSE
+ }
+
+ if(use.yourown.flagdesign){
+ use.chillports = F # change to false for using utah instead of chill portions (most models use chill portions z)
+ use.zscore = TRUE # change to false to use raw predictors
+
+ # Default is species complex and no crops
+ use.allspp = F
+ use.multcuespp = FALSE
+ use.cropspp = FALSE
+
+ # Default is species complex use alltypes of designs
+ use.expramptypes.fp = TRUE
+ use.exptypes.fp = FALSE
+
+ #Default is all chilling data
+ use.expchillonly = FALSE # change to true for only experimental chilling
+ #note: with only exp chilling, there is only exp photo and force too.
+ #also: subsetting to exp chill only reduces dataset to 3 species, <9 studies
+ }
+
+ source("..//bb_analysis/source/bbstanleadin.R")
+
+ namesdat <- unique(paste(bb.stan$genus,bb.stan$species,sep="_"))
+ bb.stan$spps <- paste(bb.stan$genus,bb.stan$species,sep="_")
+ bb.stan$phylo <- paste(bb.stan$genus,bb.stan$species,sep="_")
+
+
+ #'###################################
+ #### get phylogeny ####
+ #'###################################
+
+ setwd("..//phylogeny")
+ source("source/get_phylo_models.R")
+
+ ## read and pre-process phylogeny
+ #phylo <- read.tree("../../data/phylogeny/SBphylo_62complex.tre")
+ #phylo <- read.tree("../../data/phylogeny/SBphylo_101sps.tre")
+ phylo <- phy.plants.ospree
+
+
+ namesphy <- phylo$tip.label
+ phylo <- force.ultrametric(phylo, method="extend")
+ phylo$node.label <- seq(1,length(phylo$node.label),1)
+ is.ultrametric(phylo)
+ #plot(phylo, cex=0.7)
+ VCVPHY <- vcv.phylo(phylo,corr=TRUE)
+
+
+
+ ## deal with subgrouping
+
+ if(nocrops & agiosponly){
+ gymno <- c("Metasequoia_glyptostroboides", "Pseudotsuga_menziesii","Larix_laricina",
+ "Larix_gmelinii", "Larix_decidua" ,"Larix_kaempferi",
+ "Pinus_nigra","Pinus_sylvestris","Pinus_banksiana",
+ "Pinus_contorta","Pinus_wallichiana","Pinus_strobus",
+ "Picea_abies" ,"Picea_mariana" ,"Picea_glauca" ,
+ "Cedrus_libani" ,"Abies_alba" ,"Abies_homolepis","Ginkgo_biloba")
+ croplist <- read.csv("../../data/croplist/agricultural_species.csv")
+ cropgymno <- c(croplist$Species_name,gymno)
+ bb.stan$crops <- ifelse(bb.stan$spps %in% cropgymno, "cropgymno","nocrop")
+ cropspps <- unique(bb.stan$spps[which(bb.stan$crops=="cropgymno")])
+ bb.stan <- subset(bb.stan, crops == "nocrop")
+ phylo <- drop.tip(phylo, cropspps)
+ VCVPHY<-vcv.phylo(phylo,corr=T)
+ }
+
+ if(nocrops & !agiosponly){
+ croplist <- read.csv("../../data/croplist/agricultural_species.csv")
+ bb.stan$crops <- ifelse(bb.stan$spps %in% croplist$Species_name, "crop","nocrop")
+ cropspps <- unique(bb.stan$spps[which(bb.stan$crops=="crop")])
+ bb.stan <- subset(bb.stan, crops == "nocrop")
+ phylo <- drop.tip(phylo, cropspps)
+ VCVPHY<-vcv.phylo(phylo,corr=T)
+ }
+
+
+ if(!nocrops & agiosponly){
+ gymno <- c("Metasequoia_glyptostroboides", "Pseudotsuga_menziesii","Larix_laricina",
+ "Larix_gmelinii", "Larix_decidua" ,"Larix_kaempferi",
+ "Pinus_nigra","Pinus_sylvestris","Pinus_banksiana",
+ "Pinus_contorta","Pinus_wallichiana","Pinus_strobus",
+ "Picea_abies" ,"Picea_mariana" ,"Picea_glauca" ,
+ "Cedrus_libani" ,"Abies_alba" ,"Abies_homolepis","Ginkgo_biloba")
+ cropgymno <- c(gymno)
+ bb.stan$crops <- ifelse(bb.stan$spps %in% cropgymno, "cropgymno","nocrop")
+ cropspps <- unique(bb.stan$spps[which(bb.stan$crops=="cropgymno")])
+ bb.stan <- subset(bb.stan, crops == "nocrop")
+ phylo <- drop.tip(phylo, cropspps)
+ VCVPHY<-vcv.phylo(phylo,corr=T)
+ }
+
+
+ # Get spps and VCVPHY in same order
+ # bb.stan$spps[phylo$tip.label]
+ phylo$tip.label
+ d <- bb.stan[match(phylo$tip.label, bb.stan$spps),] # hmmm, only gives ONE match
+
+ phymatch <- data.frame(tip=phylo$tip.label, sppnum=c(1:length(phylo$tip.label)))
+ d <- merge(bb.stan, phymatch, by.x="spps", by.y="tip")
+ d <- d[order(d$sppnum),]
+ # Tilia_cordata versus Tilia_Cordata in phylo
+ nspecies <- max(d$sppnum)
+
+ ## remove outliers
+ # d$resp
+ head(d)
+ ff = subset(d,latbi %in% c("Populus_balsamifera","Populus_tremuloides"))
+ d = subset(d,!latbi %in% c("Populus_balsamifera","Populus_tremuloides"))
+ nspecies = 192
+ phylo <- drop.tip(phylo, c("Populus_balsamifera","Populus_tremuloides"))
+ d$sppnum <- as.numeric(as.factor(d$sppnum))
+
+
+ ## remove names of species that are wrong (e.g. Acer pseudolatanus) Malyshev2018
+ idswrong = which(d$spps == "Acer_pseudolatauns")
+ d$spps[idswrong] = "Acer_pseudoplatanus"
+ d$species[idswrong] = "pseudoplatanus"
+ d$latbi[idswrong] = "Acer_pseudoplatanus"
+ d$phylo[idswrong] = "Acer_pseudoplatanus"
+
+ #d$sppnum[which(d$latbi=="Acer_pseudoplatanus")]
+ d$sppnum[idswrong] = 127
+ d$sppnum[which(d$sppnum>137)] = d$sppnum[which(d$sppnum>137)]-1
+
+ nspecies = 191
+ phylo <- drop.tip(phylo, "Acer_pseudolatauns")
+
+
+
+ ## remove names of species that are wrong (e.g. Juglans spp)
+ idswrong = which(d$spps == "Juglans_spp")
+ d <- d[-idswrong,]
+ phylo <- drop.tip(phylo, "Juglans_spp")
+
+ nspecies <- length(phylo$tip.label)
+ phymatch2 <- data.frame(tip=phylo$tip.label, sppnum=c(1:length(phylo$tip.label)))
+ d2 <- merge(d, phymatch2, by.x="spps", by.y="tip")
+ d2 <- d2[order(d2$sppnum.y),]
+ d2$sppnum <- d2$sppnum.y
+ d <- d2
+
+ # save for faster loading-compiling
+ #write.csv(d,file = "input/datasetforphyloms.csv")
+ #write.tree(phylo, file = "input/phyloforphyloms.tre")
+
+ } else {
+
+ d = read.csv("input/datasetforphyloms.csv")
+ phylo = read.tree("input/phyloforphyloms.tre")
+
+ }
Joining with `by = join_by(genus, species)`
Joining with `by = join_by(genus, species)`
Joining with `by = join_by(genus, species)`
Joining with `by = join_by(genus, species)`
List of 7
$ y : num [1:4442] 56.4 52.2 59.6 49.2 61.5 ...
$ chill: num [1:4442] 0.05553 0.00235 0.05553 0.00235 0.05553 ...
$ force: num [1:4442] -1.8 -1.8 -1.8 -1.8 -1.8 ...
$ photo: num [1:4442] -0.725 -0.438 -0.725 -0.438 -0.725 ...
$ sp : num [1:4442] 1 1 1 1 1 1 1 1 1 1 ...
$ N : int 4442
$ n_sp : int 233
[1] "Unique forcing types are ..."
[1] "ramped" "exp" "amb"
[1] "Unique photo types are ..."
[1] "ramped" "exp" "amb" "none"
[1] "Unique chilling types are ..."
[1] "fldest" "bothest" "bothrep" "exp"
[1] "Range of chilling is ..."
[1] -3.274617 4.558866
Number of species in tree before: 0
Number of species in tree now: 24
***************************************************************
* Note: *
* force.ultrametric does not include a formal method to *
* ultrametricize a tree & should only be used to coerce *
* a phylogeny that fails is.ultramtric due to rounding -- *
* not as a substitute for formal rate-smoothing methods. *
***************************************************************
Warning messages:
1: In eval(ei, envir) : NAs introduced by coercion
2: In eval(ei, envir) : NAs introduced by coercion
3: In eval(ei, envir) : NAs introduced by coercion
4: In eval(ei, envir) : NAs introduced by coercion
>
> if(runwithchillports){
+
+ d$chill.z = as.numeric(scale(d$chill.ports))
+
+ }
> runwithchillports
[1] TRUE
> fit <- stan("stan/uber_threeslopeintercept_modified_cholesky_updatedpriors.stan",
+ data=list(N=nrow(d),
+ n_sp=nspecies,
+ sp=d$sppnum,
+ x1=d$force.z,
+ x2 = d$chill.z,
+ x3=d$photo.z,
+ y=d$resp,
+ Vphy=vcv(phylo, corr = TRUE)),
+ iter = 4000,
+ warmup = 2000, # half the iter as warmp is default, but leaving in case we want to change
+ chains = 4,
+ seed = 117
+ )
Warning message:
Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#bulk-ess
> 
MoralesCastilla
@MoralesCastilla I am fairly sure you lose a little data in the switch to chill portions, but not so much.
I ran the model using chill portions, I'll paste my code below so you can see if I ran it correctly, but I did not get output similar to yours. Here's what I got: