##EGF Modeling, hopefully final, April 2017 #Packages library(ape) library(caper) library(MCMCglmm) library(mulTree) library(MuMIn) library(fmsb) #Loading and prepping data setwd("D:/copyofgdrive/Konstanz/Inventory") #setwd("E:/Inventory") #setwd("C:/Users/Ha/Google Drive/Konstanz/Inventory") #load data df00<-read.csv("modelingEGFupdatedfulldata_2017.csv", header=T, na.strings=c("NA", "", "NO MATCH")) #should just be called species df00$Species<-df00$EGF.Name #remove extra columns #df<-df00[c(4:6, 10:13, 15, 17:19, 22:23, 26:27, 29:46, 50:52)] #rm(df00) #remove remaining european natives df<-df00[which(df00$Eur_incl=="other"),] rm(df00) #Naturalized in europe yes or no df$Nat_in_E <- df$Glonaf_Europe_region_num df$Nat_in_E[df$Glonaf_Europe_region_num > 0] <- 1 ###########Standardize continuous vars #stdize glonaf regions #z<-df$Non.Europe_region_num #x<-mean(na.omit(z)) #y<-sd(na.omit(z)) #df$NE_gl_reg_sd<-(z-x)/y df$NE_gl_reg<-df$Non.Europe_region_num #stdiz height df$Height<-as.numeric(df$Height_max_max) #z<-df$Height #x<-mean(na.omit(z)) #y<-sd(na.omit(z)) #df$Height_sd<-(z-x)/y #stdiz nursery sales #z<-df$VDV_all #x<-mean(na.omit(z)) #y<-sd(na.omit(z)) #df$VDV_sd<-(z-x)/y #stdiz CSI curmedian #z<-df$curmedian #x<-mean(na.omit(z)) #y<-sd(na.omit(z)) #df$curmedian_sd<-(z-x)/y #stdiz CSI curmean #z<-df$curmean #x<-mean(na.omit(z)) #y<-sd(na.omit(z)) #df$curmean_sd<-(z-x)/y #stdiz CSI curmax #z<-df$curmax #x<-mean(na.omit(z)) #y<-sd(na.omit(z)) #df$curmax_sd<-(z-x)/y #stdiz CSI curquant75 #z<-df$curquant75 #x<-mean(na.omit(z)) #y<-sd(na.omit(z)) #df$curquant75_sd<-(z-x)/y #stdiz CSI curquant95 #z<-df$curquant95 #x<-mean(na.omit(z)) #y<-sd(na.omit(z)) #df$curquant95_sd<-(z-x)/y #native range size #z<-df$Combo_NRS #x<-mean(na.omit(z)) #y<-sd(na.omit(z)) #df$NRS_sd<-(z-x)/y #hardiness zone df$HZ<-df$Hard1 levels(df$Hard1)<-c(2,1,7,6,5,4,3) df$HZ<-as.numeric(as.factor(df$Hard1)) #z<-df$HZ #x<-mean(na.omit(z)) #y<-sd(na.omit(z)) ###should I standardize HZ (divide by sd)? Or merely center? #df$HZ_sd<-(z-x)/y ########################### #centering categorical variables ###for some reason all the zeros are registering as 2's, so just swap the 0's and 1's? #propagation method-- converted to "seed y/n" and "veg y/n" df$PropSeed<-df$Propag df$PropVeg<-df$Propag levels(df$PropSeed)<-c(1, 1, 0) levels(df$PropVeg)<-c(1, 0, 1) df$PropSeed<-as.numeric(levels(df$PropSeed))[df$PropSeed] df$PropVeg<-as.numeric(levels(df$PropVeg))[df$PropVeg] #z<-df$PropSeed #x<-mean(na.omit(z)) #df$PropSeed_sd<-z-x #z<-df$PropVeg #x<-mean(na.omit(z)) #df$PropVeg_sd<-z-x ###Life form--herbaceous v wody df$Herb<-df$LF_APW levels(df$Herb)<-c(1, 0, 0) df$Herb<-as.numeric(levels(df$Herb))[df$Herb] #z<-df$Herb #x<-mean(na.omit(z)) #df$Herb_sd<-z-x ##Sex--mono v dio ###Life form--herbaceous v wody df$Mono<-df$Sex levels(df$Mono)<-c(0, 1, 1, 1) df$Mono<-as.numeric(levels(df$Mono))[df$Mono] #z<-df$Mono #x<-mean(na.omit(z)) #df$Mono_sd<-z-x ##Storage organ presence-- y/n df$SO<-df$Organ_simp levels(df$SO)<-c(0, 1) df$SO<-as.numeric(levels(df$SO))[df$SO] #z<-df$SO #x<-mean(na.omit(z)) #df$SO_sd<-z-x ##Hybridized-- y/n (might not use this, not as much data, not many hybrids) df$hybr<-df$hybrid levels(df$hybr)<-c(1, 0) df$hybr<-as.numeric(levels(df$hybr))[df$hybr] #z<-df$hybr #x<-mean(na.omit(z)) #df$hybr_sd<-z-x #hist(df$hybr_sd) ##RarelyPlanted-- y/n (might not use this, not as much data, not many rps) df$RP<-df$RarelyPlanted levels(df$RP)<-c(0, 1) df$RP<-as.numeric(levels(df$RP))[df$RP] #z<-df$RP #x<-mean(na.omit(z)) #df$RP_sd<-z-x #Should also try it without: Mono_sd, (and native range unless I can find more) #all others only have a few NAs ##removing genera not available in the phylogeny df1<-df df1<-df1[which(df1$Genus!="Anoda"),] df1<-df1[which(df1$Genus!="Balsamorhiza"),] df1<-df1[which(df1$Genus!="Barnardia"),] df1<-df1[which(df1$Genus!="Bartlettina"),] df1<-df1[which(df1$Genus!="Beschorneria"),] df1<-df1[which(df1$Genus!="Chondrosum"),] df1<-df1[which(df1$Genus!="Fittonia"),] df1<-df1[which(df1$Genus!="Haplopappus"),] df1<-df1[which(df1$Genus!="Hoodia"),] df1<-df1[which(df1$Genus!="Huernia"),] df1<-df1[which(df1$Genus!="Orbea"),] df1<-df1[which(df1$Genus!="Pachystachys"),] df1<-df1[which(df1$Genus!="Parahebe"),] df1<-df1[which(df1$Genus!="Pediocactus"),] df1<-df1[which(df1$Genus!="Piaranthus"),] df1<-df1[which(df1$Genus!="Sanvitalia"),] df1<-df1[which(df1$Genus!="Sidalcea"),] df1<-df1[which(df1$Genus!="Synthyris"),] ###editing to fit in tree df1$Species<-gsub(" ", "_", df1$Species) ###########Remove excess columns df2<-df1[c(5, 19,21, 24, 26, 29, 31, 34:46)] #complete species list for creating full tree #df4<-df1[c(35:36)] ##Further edits for accomodating the phylogeny #setwd("D:/copyofgdrive/Konstanz/Inventory/Trees") #setwd("E:/Inventory/Trees") #tree1<-read.tree(file = "C:/Users/Ha/Google Drive/Konstanz/Inventory/Trees/whoisnext_chronogram_sp_polytomy1.tre", tree1<-read.tree(file = "D:/copyofgdrive/Konstanz/Inventory/Trees/whoisnext_chronogram_sp_polytomy1.tre", text = NULL, tree.names = NULL, skip = 0, comment.char = "#", keep.multi = FALSE) tree1<-as.phylo(tree1) tree1$node.label<-1:1242 #any(duplicated(tree1$node.label)) ###################### #creating the data for MCMCglmm #comp_data1 <- comparative.data(phy = tree1, # data =yyy, # names.col = Species, # vcv=TRUE) yyy<-df2[which(df2$Species!="Persicaria_virginiana"),] #for some reason there's a problem and this one needs to be removed? mcmc_data1 <- cbind(animal = yyy$Species, yyy) #mcmc_tree1 <- comp_data1$phy mcmc_tree1 <- tree1 #Create prior for the binomial mcmcglmm's prior1= list(R = list(V = 1, fix = 1), G = list(G1 = list(V = 1, nu = 0.002))) #update MCMCglmm function so that I can dredge it later uMCMCglmm <- updateable(MCMCglmm) #MCMCglmm <- updateable(MCMCglmm) ############# ####creating NA-free datasets with which to create these models dat_a<-na.omit(mcmc_data1[c(1,7,2)]) dat_b<-na.omit(mcmc_data1[c(1,7,3)]) dat_c<-na.omit(mcmc_data1[c(1,7,4)]) dat_d<-na.omit(mcmc_data1[c(1,7,5)]) dat_e<-na.omit(mcmc_data1[c(1,7,8)]) dat_f<-na.omit(mcmc_data1[c(1,7,9)]) dat_g<-na.omit(mcmc_data1[c(1,7,10)]) dat_h<-na.omit(mcmc_data1[c(1,7,11)]) dat_i<-na.omit(mcmc_data1[c(1,7,12)]) dat_j<-na.omit(mcmc_data1[c(1,7,13)]) dat_k<-na.omit(mcmc_data1[c(1,7,14)]) dat_l<-na.omit(mcmc_data1[c(1,7,15)]) dat_m<-na.omit(mcmc_data1[c(1,7,16)]) dat_n<-na.omit(mcmc_data1[c(1,7,17)]) dat_a0<-na.omit(mcmc_data1[c(1,7,2,3)]) dat_b0<-na.omit(mcmc_data1[c(1,7,2,4)]) dat_c0<-na.omit(mcmc_data1[c(1,7,2,5)]) dat_d0<-na.omit(mcmc_data1[c(1,7,2,8)]) dat_e0<-na.omit(mcmc_data1[c(1,7,2,9)]) dat_f0<-na.omit(mcmc_data1[c(1,7,2,10)]) dat_g0<-na.omit(mcmc_data1[c(1,7,2,11)]) dat_h0<-na.omit(mcmc_data1[c(1,7,2,12)]) dat_i0<-na.omit(mcmc_data1[c(1,7,2,13)]) dat_j0<-na.omit(mcmc_data1[c(1,7,2,14)]) dat_k0<-na.omit(mcmc_data1[c(1,7,2,15)]) dat_l0<-na.omit(mcmc_data1[c(1,7,2,16)]) dat_n0<-na.omit(mcmc_data1[c(1,7,2,17)]) dat_a1<-na.omit(mcmc_data1[c(1,7,3,4)]) dat_b1<-na.omit(mcmc_data1[c(1,7,3,5)]) dat_c1<-na.omit(mcmc_data1[c(1,7,3,8)]) dat_d1<-na.omit(mcmc_data1[c(1,7,3,9)]) dat_e1<-na.omit(mcmc_data1[c(1,7,3,10)]) dat_f1<-na.omit(mcmc_data1[c(1,7,3,11)]) dat_g1<-na.omit(mcmc_data1[c(1,7,3,12)]) dat_h1<-na.omit(mcmc_data1[c(1,7,3,13)]) dat_i1<-na.omit(mcmc_data1[c(1,7,3,14)]) dat_j1<-na.omit(mcmc_data1[c(1,7,3,15)]) dat_k1<-na.omit(mcmc_data1[c(1,7,3,16)]) dat_l1<-na.omit(mcmc_data1[c(1,7,3,17)]) dat_a2<-na.omit(mcmc_data1[c(1,7,4,5)]) dat_b2<-na.omit(mcmc_data1[c(1,7,4,8)]) dat_c2<-na.omit(mcmc_data1[c(1,7,4,9)]) dat_d2<-na.omit(mcmc_data1[c(1,7,4,10)]) dat_e2<-na.omit(mcmc_data1[c(1,7,4,11)]) dat_f2<-na.omit(mcmc_data1[c(1,7,4,12)]) dat_g2<-na.omit(mcmc_data1[c(1,7,4,13)]) dat_h2<-na.omit(mcmc_data1[c(1,7,4,14)]) dat_i2<-na.omit(mcmc_data1[c(1,7,4,15)]) dat_j2<-na.omit(mcmc_data1[c(1,7,4,16)]) dat_k2<-na.omit(mcmc_data1[c(1,7,4,17)]) dat_a3<-na.omit(mcmc_data1[c(1,7,5, 8)]) dat_b3<-na.omit(mcmc_data1[c(1,7,5, 9)]) dat_c3<-na.omit(mcmc_data1[c(1,7,5,10)]) dat_d3<-na.omit(mcmc_data1[c(1,7,5,11)]) dat_e3<-na.omit(mcmc_data1[c(1,7,5,12)]) dat_f3<-na.omit(mcmc_data1[c(1,7,5,13)]) dat_g3<-na.omit(mcmc_data1[c(1,7,5,14)]) dat_h3<-na.omit(mcmc_data1[c(1,7,5,15)]) dat_i3<-na.omit(mcmc_data1[c(1,7,5,16)]) dat_j3<-na.omit(mcmc_data1[c(1,7,5,17)]) dat_a4<-na.omit(mcmc_data1[c(1,7,8,9)]) dat_b4<-na.omit(mcmc_data1[c(1,7,8,10)]) dat_c4<-na.omit(mcmc_data1[c(1,7,8,11)]) dat_d4<-na.omit(mcmc_data1[c(1,7,8,12)]) dat_e4<-na.omit(mcmc_data1[c(1,7,8,13)]) dat_f4<-na.omit(mcmc_data1[c(1,7,8,14)]) dat_g4<-na.omit(mcmc_data1[c(1,7,8,15)]) dat_h4<-na.omit(mcmc_data1[c(1,7,8,16)]) dat_i4<-na.omit(mcmc_data1[c(1,7,8,17)]) dat_a5<-na.omit(mcmc_data1[c(1,7,9,10)]) dat_b5<-na.omit(mcmc_data1[c(1,7,9,11)]) dat_c5<-na.omit(mcmc_data1[c(1,7,9,12)]) dat_d5<-na.omit(mcmc_data1[c(1,7,9,13)]) dat_e5<-na.omit(mcmc_data1[c(1,7,9,14)]) dat_f5<-na.omit(mcmc_data1[c(1,7,9,15)]) dat_g5<-na.omit(mcmc_data1[c(1,7,9,16)]) dat_h5<-na.omit(mcmc_data1[c(1,7,9,17)]) dat_a6<-na.omit(mcmc_data1[c(1,7,10,11)]) dat_b6<-na.omit(mcmc_data1[c(1,7,10,12)]) dat_c6<-na.omit(mcmc_data1[c(1,7,10,13)]) dat_d6<-na.omit(mcmc_data1[c(1,7,10,14)]) dat_e6<-na.omit(mcmc_data1[c(1,7,10,15)]) dat_f6<-na.omit(mcmc_data1[c(1,7,10,16)]) dat_g6<-na.omit(mcmc_data1[c(1,7,10,17)]) dat_a7<-na.omit(mcmc_data1[c(1,7,11,13)]) dat_b7<-na.omit(mcmc_data1[c(1,7,11,14)]) dat_c7<-na.omit(mcmc_data1[c(1,7,11,15)]) dat_d7<-na.omit(mcmc_data1[c(1,7,11,16)]) dat_e7<-na.omit(mcmc_data1[c(1,7,11,17)]) dat_a8<-na.omit(mcmc_data1[c(1,7,12,13)]) dat_b8<-na.omit(mcmc_data1[c(1,7,12,14)]) dat_c8<-na.omit(mcmc_data1[c(1,7,12,15)]) dat_d8<-na.omit(mcmc_data1[c(1,7,12,16)]) dat_e8<-na.omit(mcmc_data1[c(1,7,12,17)]) dat_a9<-na.omit(mcmc_data1[c(1,7,13,14)]) dat_b9<-na.omit(mcmc_data1[c(1,7,13,15)]) dat_c9<-na.omit(mcmc_data1[c(1,7,13,16)]) dat_d9<-na.omit(mcmc_data1[c(1,7,13,17)]) dat_a10<-na.omit(mcmc_data1[c(1,7,14,15)]) dat_b10<-na.omit(mcmc_data1[c(1,7,14,16)]) dat_c10<-na.omit(mcmc_data1[c(1,7,14,17)]) dat_a11<-na.omit(mcmc_data1[c(1,7,15,16)]) dat_b11<-na.omit(mcmc_data1[c(1,7,15,17)]) dat_a12<-na.omit(mcmc_data1[c(1,7,16,17)]) ################## #Formulas for singleterm/interaction models for dredging/AIC comparisons formula_NULL<-Nat_in_E~1 #Null model for comparing AIC's against formula_a<-Nat_in_E~scale(VDV_all) formula_b<-Nat_in_E~scale(curmax) formula_c<-Nat_in_E~scale(curmedian) formula_d<-Nat_in_E~scale(Combo_NRS) formula_e<-Nat_in_E~scale(NE_gl_reg) formula_f<-Nat_in_E~scale(Height) formula_g<-Nat_in_E~scale(HZ) formula_h<-Nat_in_E~scale(PropSeed) formula_i<-Nat_in_E~scale(PropVeg) formula_j<-Nat_in_E~scale(Herb) formula_k<-Nat_in_E~scale(Mono) formula_l<-Nat_in_E~scale(SO) formula_m<-Nat_in_E~scale(hybr) formula_n<-Nat_in_E~scale(RP) formula_a0<-Nat_in_E~scale(VDV_all)*scale(curmax) formula_b0<-Nat_in_E~scale(VDV_all)*scale(curmedian) formula_c0<-Nat_in_E~scale(VDV_all)*scale(Combo_NRS) formula_d0<-Nat_in_E~scale(VDV_all)*scale(NE_gl_reg) formula_e0<-Nat_in_E~scale(VDV_all)*scale(Height) formula_f0<-Nat_in_E~scale(VDV_all)*scale(HZ) formula_g0<-Nat_in_E~scale(VDV_all)*scale(PropSeed) formula_h0<-Nat_in_E~scale(VDV_all)*scale(PropVeg) formula_i0<-Nat_in_E~scale(VDV_all)*scale(Herb) formula_j0<-Nat_in_E~scale(VDV_all)*scale(Mono) formula_k0<-Nat_in_E~scale(VDV_all)*scale(SO) formula_l0<-Nat_in_E~scale(VDV_all)*scale(hybr) formula_n0<-Nat_in_E~scale(VDV_all)*scale(RP) formula_a1<-Nat_in_E~scale(curmax)*scale(curmedian) formula_b1<-Nat_in_E~scale(curmax)*scale(Combo_NRS) formula_c1<-Nat_in_E~scale(curmax)*scale(NE_gl_reg) formula_d1<-Nat_in_E~scale(curmax)*scale(Height) formula_e1<-Nat_in_E~scale(curmax)*scale(HZ) formula_f1<-Nat_in_E~scale(curmax)*scale(PropSeed) formula_g1<-Nat_in_E~scale(curmax)*scale(PropVeg) formula_h1<-Nat_in_E~scale(curmax)*scale(Herb) formula_i1<-Nat_in_E~scale(curmax)*scale(Mono) formula_j1<-Nat_in_E~scale(curmax)*scale(SO) formula_k1<-Nat_in_E~scale(curmax)*scale(hybr) formula_l1<-Nat_in_E~scale(curmax)*scale(RP) formula_a2<-Nat_in_E~scale(curmedian)*scale(Combo_NRS) formula_b2<-Nat_in_E~scale(curmedian)*scale(NE_gl_reg) formula_c2<-Nat_in_E~scale(curmedian)*scale(Height) formula_d2<-Nat_in_E~scale(curmedian)*scale(HZ) formula_e2<-Nat_in_E~scale(curmedian)*scale(PropSeed) formula_f2<-Nat_in_E~scale(curmedian)*scale(PropVeg) formula_g2<-Nat_in_E~scale(curmedian)*scale(Herb) formula_h2<-Nat_in_E~scale(curmedian)*scale(Mono) formula_i2<-Nat_in_E~scale(curmedian)*scale(SO) formula_j2<-Nat_in_E~scale(curmedian)*scale(hybr) formula_k2<-Nat_in_E~scale(curmedian)*scale(RP) formula_a3<-Nat_in_E~scale(Combo_NRS)*scale(NE_gl_reg) formula_b3<-Nat_in_E~scale(Combo_NRS)*scale(Height) formula_c3<-Nat_in_E~scale(Combo_NRS)*scale(HZ) formula_d3<-Nat_in_E~scale(Combo_NRS)*scale(PropSeed) formula_e3<-Nat_in_E~scale(Combo_NRS)*scale(PropVeg) formula_f3<-Nat_in_E~scale(Combo_NRS)*scale(Herb) formula_g3<-Nat_in_E~scale(Combo_NRS)*scale(Mono) formula_h3<-Nat_in_E~scale(Combo_NRS)*scale(SO) formula_i3<-Nat_in_E~scale(Combo_NRS)*scale(hybr) formula_j3<-Nat_in_E~scale(Combo_NRS)*scale(RP) formula_a4<-Nat_in_E~scale(NE_gl_reg)*scale(Height) formula_b4<-Nat_in_E~scale(NE_gl_reg)*scale(HZ) formula_c4<-Nat_in_E~scale(NE_gl_reg)*scale(PropSeed) formula_d4<-Nat_in_E~scale(NE_gl_reg)*scale(PropVeg) formula_e4<-Nat_in_E~scale(NE_gl_reg)*scale(Herb) formula_f4<-Nat_in_E~scale(NE_gl_reg)*scale(Mono) formula_g4<-Nat_in_E~scale(NE_gl_reg)*scale(SO) formula_h4<-Nat_in_E~scale(NE_gl_reg)*scale(hybr) formula_i4<-Nat_in_E~scale(NE_gl_reg)*scale(RP) formula_a5<-Nat_in_E~scale(Height)*scale(HZ) formula_b5<-Nat_in_E~scale(Height)*scale(PropSeed) formula_c5<-Nat_in_E~scale(Height)*scale(PropVeg) formula_d5<-Nat_in_E~scale(Height)*scale(Herb) formula_e5<-Nat_in_E~scale(Height)*scale(Mono) formula_f5<-Nat_in_E~scale(Height)*scale(SO) formula_g5<-Nat_in_E~scale(Height)*scale(hybr) formula_h5<-Nat_in_E~scale(Height)*scale(RP) formula_a6<-Nat_in_E~scale(HZ)*scale(PropSeed) formula_b6<-Nat_in_E~scale(HZ)*scale(PropVeg) formula_c6<-Nat_in_E~scale(HZ)*scale(Herb) formula_d6<-Nat_in_E~scale(HZ)*scale(Mono) formula_e6<-Nat_in_E~scale(HZ)*scale(SO) formula_f6<-Nat_in_E~scale(HZ)*scale(hybr) formula_g6<-Nat_in_E~scale(HZ)*scale(RP) formula_a7<-Nat_in_E~scale(PropSeed)*scale(Herb) formula_b7<-Nat_in_E~scale(PropSeed)*scale(Mono) formula_c7<-Nat_in_E~scale(PropSeed)*scale(SO) formula_d7<-Nat_in_E~scale(PropSeed)*scale(hybr) formula_e7<-Nat_in_E~scale(PropSeed)*scale(RP) formula_a8<-Nat_in_E~scale(PropVeg)*scale(Herb) formula_b8<-Nat_in_E~scale(PropVeg)*scale(Mono) formula_c8<-Nat_in_E~scale(PropVeg)*scale(SO) formula_d8<-Nat_in_E~scale(PropVeg)*scale(hybr) formula_e8<-Nat_in_E~scale(PropVeg)*scale(RP) formula_a9<-Nat_in_E~scale(Herb)*scale(Mono) formula_b9<-Nat_in_E~scale(Herb)*scale(SO) formula_c9<-Nat_in_E~scale(Herb)*scale(hybr) formula_d9<-Nat_in_E~scale(Herb)*scale(RP) formula_a10<-Nat_in_E~scale(Mono)*scale(SO) formula_b10<-Nat_in_E~scale(Mono)*scale(hybr) formula_c10<-Nat_in_E~scale(Mono)*scale(RP) formula_a11<-Nat_in_E~scale(SO)*scale(hybr) formula_b11<-Nat_in_E~scale(SO)*scale(RP) formula_a12<-Nat_in_E~scale(hybr)*scale(RP) formula_a13<-Nat_in_E~scale(VDV_all) + I(scale(VDV_all)^2) formula_b13<-Nat_in_E~scale(curmax) + I(scale(curmax)^2) formula_c13<-Nat_in_E~scale(curmedian) + I(scale(curmedian)^2) formula_d13<-Nat_in_E~scale(Combo_NRS) + I(scale(Combo_NRS)^2) formula_e13<-Nat_in_E~scale(NE_gl_reg) + I(scale(NE_gl_reg)^2) formula_f13<-Nat_in_E~scale(Height) + I(scale(Height)^2) formula_g13<-Nat_in_E~scale(HZ) + I(scale(HZ)^2) ################### #RUNNING THE MODELS mod_NULL<- uMCMCglmm(fixed = formula_NULL, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = mcmc_data1, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_a<- uMCMCglmm(fixed = formula_a, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_a, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_b<- uMCMCglmm(fixed = formula_b, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_b, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_c<- uMCMCglmm(fixed = formula_c, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_c, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_d<- uMCMCglmm(fixed = formula_d, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_d, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_e<- uMCMCglmm(fixed = formula_e, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_e, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_f<- uMCMCglmm(fixed = formula_f, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_f, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_g<- uMCMCglmm(fixed = formula_g, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_g, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_h<- uMCMCglmm(fixed = formula_h, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_h, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_i<- uMCMCglmm(fixed = formula_i, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_i, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_j<- uMCMCglmm(fixed = formula_j, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_j, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_k<- uMCMCglmm(fixed = formula_k, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_k, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_l<- uMCMCglmm(fixed = formula_l, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_l, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_m<- uMCMCglmm(fixed = formula_m, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_m, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_n<- uMCMCglmm(fixed = formula_n, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_n, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_a0<- uMCMCglmm(fixed = formula_a0, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_a0, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_b0<- uMCMCglmm(fixed = formula_b0, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_b0, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_c0<- uMCMCglmm(fixed = formula_c0, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_c0, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_d0<- uMCMCglmm(fixed = formula_d0, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_d0, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_e0<- uMCMCglmm(fixed = formula_e0, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_e0, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_f0<- uMCMCglmm(fixed = formula_f0, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_f0, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_g0<- uMCMCglmm(fixed = formula_g0, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_g0, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_h0<- uMCMCglmm(fixed = formula_h0, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_h0, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_i0<- uMCMCglmm(fixed = formula_i0, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_i0, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_j0<- uMCMCglmm(fixed = formula_j0, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_j0, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_k0<- uMCMCglmm(fixed = formula_k0, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_k0, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_l0<- uMCMCglmm(fixed = formula_l0, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_l0, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_n0<- uMCMCglmm(fixed = formula_n0, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_n0, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_a1<- uMCMCglmm(fixed = formula_a1, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_a1, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_b1<- uMCMCglmm(fixed = formula_b1, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_b1, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_c1<- uMCMCglmm(fixed = formula_c1, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_c1, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_d1<- uMCMCglmm(fixed = formula_d1, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_d1, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_e1<- uMCMCglmm(fixed = formula_e1, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_e1, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_f1<- uMCMCglmm(fixed = formula_f1, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_f1, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_g1<- uMCMCglmm(fixed = formula_g1, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_g1, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_h1<- uMCMCglmm(fixed = formula_h1, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_h1, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_i1<- uMCMCglmm(fixed = formula_i1, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_i1, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_j1<- uMCMCglmm(fixed = formula_j1, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_j1, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_k1<- uMCMCglmm(fixed = formula_k1, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_k1, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_l1<- uMCMCglmm(fixed = formula_l1, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_l1, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_a2<- uMCMCglmm(fixed = formula_a2, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_a2, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_b2<- uMCMCglmm(fixed = formula_b2, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_b2, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_c2<- uMCMCglmm(fixed = formula_c2, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_c2, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_d2<- uMCMCglmm(fixed = formula_d2, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_d2, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_e2<- uMCMCglmm(fixed = formula_e2, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_e2, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_f2<- uMCMCglmm(fixed = formula_f2, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_f2, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_g2<- uMCMCglmm(fixed = formula_g2, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_g2, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_h2<- uMCMCglmm(fixed = formula_h2, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_h2, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_i2<- uMCMCglmm(fixed = formula_i2, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_i2, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_j2<- uMCMCglmm(fixed = formula_j2, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_j2, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_k2<- uMCMCglmm(fixed = formula_k2, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_k2, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_a3<- uMCMCglmm(fixed = formula_a3, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_a3, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_b3<- uMCMCglmm(fixed = formula_b3, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_b3, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_c3<- uMCMCglmm(fixed = formula_c3, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_c3, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_d3<- uMCMCglmm(fixed = formula_d3, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_d3, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_e3<- uMCMCglmm(fixed = formula_e3, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_e3, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_f3<- uMCMCglmm(fixed = formula_f3, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_f3, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_g3<- uMCMCglmm(fixed = formula_g3, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_g3, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_h3<- uMCMCglmm(fixed = formula_h3, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_h3, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_i3<- uMCMCglmm(fixed = formula_i3, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_i3, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_j3<- uMCMCglmm(fixed = formula_j3, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_j3, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_a4<- uMCMCglmm(fixed = formula_a4, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_a4, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_b4<- uMCMCglmm(fixed = formula_b4, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_b4, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_c4<- uMCMCglmm(fixed = formula_c4, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_c4, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_d4<- uMCMCglmm(fixed = formula_d4, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_d4, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_e4<- uMCMCglmm(fixed = formula_e4, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_e4, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_f4<- uMCMCglmm(fixed = formula_f4, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_f4, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_g4<- uMCMCglmm(fixed = formula_g4, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_g4, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_h4<- uMCMCglmm(fixed = formula_h4, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_h4, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_i4<- uMCMCglmm(fixed = formula_i4, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_i4, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_a5<- uMCMCglmm(fixed = formula_a5, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_a5, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_b5<- uMCMCglmm(fixed = formula_b5, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_b5, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_c5<- uMCMCglmm(fixed = formula_c5, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_c5, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_d5<- uMCMCglmm(fixed = formula_d5, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_d5, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_e5<- uMCMCglmm(fixed = formula_e5, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_e5, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_f5<- uMCMCglmm(fixed = formula_f5, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_f5, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_g5<- uMCMCglmm(fixed = formula_g5, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_g5, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_h5<- uMCMCglmm(fixed = formula_h5, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_h5, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_a6<- uMCMCglmm(fixed = formula_a6, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_a6, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_b6<- uMCMCglmm(fixed = formula_b6, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_b6, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_c6<- uMCMCglmm(fixed = formula_c6, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_c6, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_d6<- uMCMCglmm(fixed = formula_d6, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_d6, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_e6<- uMCMCglmm(fixed = formula_e6, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_e6, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_f6<- uMCMCglmm(fixed = formula_f6, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_f6, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_g6<- uMCMCglmm(fixed = formula_g6, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_g6, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_a7<- uMCMCglmm(fixed = formula_a7, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_a7, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_b7<- uMCMCglmm(fixed = formula_b7, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_b7, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_c7<- uMCMCglmm(fixed = formula_c7, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_c7, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_d7<- uMCMCglmm(fixed = formula_d7, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_d7, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_e7<- uMCMCglmm(fixed = formula_e7, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_e7, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_a8<- uMCMCglmm(fixed = formula_a8, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_a8, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_b8<- uMCMCglmm(fixed = formula_b8, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_b8, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_c8<- uMCMCglmm(fixed = formula_c8, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_c8, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_d8<- uMCMCglmm(fixed = formula_d8, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_d8, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_e8<- uMCMCglmm(fixed = formula_e8, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_e8, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_a9<- uMCMCglmm(fixed = formula_a9, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_a9, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_b9<- uMCMCglmm(fixed = formula_b9, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_b9, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_c9<- uMCMCglmm(fixed = formula_c9, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_c9, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_d9<- uMCMCglmm(fixed = formula_d9, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_d9, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_a10<- uMCMCglmm(fixed = formula_a10, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_a10, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_b10<- uMCMCglmm(fixed = formula_b10, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_b10, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_c10<- uMCMCglmm(fixed = formula_c10, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_c10, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_a11<- uMCMCglmm(fixed = formula_a11, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_a11, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_b11<- uMCMCglmm(fixed = formula_b11, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_b11, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_a12<- uMCMCglmm(fixed = formula_a12, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_a12, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_a13<- uMCMCglmm(fixed = formula_a13, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_a, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_b13<- uMCMCglmm(fixed = formula_b13, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_b, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_c13<- uMCMCglmm(fixed = formula_c13, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_c, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_d13<- uMCMCglmm(fixed = formula_d13, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_d, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_e13<- uMCMCglmm(fixed = formula_e13, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_e, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_f13<- uMCMCglmm(fixed = formula_f13, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_f, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_g13<- uMCMCglmm(fixed = formula_g13, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_g, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) #save workspace image save.image("D:/copyofgdrive/Konstanz/Inventory/EGF_initial_models_06062017.RData") ############################## #DREDGING ############################## ###Parallelizing the dredge function, to speed it up? library(parallel) no_cores <- detectCores() - 1 cl<-makeCluster(no_cores) clusterEvalQ(cl, library(MuMIn)) clusterEvalQ(cl, library(MCMCglmm)) clusterExport(cl, c('formula_a0', 'mcmc_data1', 'prior1', 'mcmc_tree1', 'uMCMCglmm', ###Data "dat_a", "dat_b", "dat_c", "dat_d", "dat_e", "dat_f", "dat_g", "dat_h", "dat_i", "dat_j", "dat_k", "dat_l", "dat_m", "dat_n", "dat_a0", "dat_b0", "dat_c0", "dat_d0", "dat_e0", "dat_f0", "dat_g0", "dat_h0", "dat_i0", "dat_j0", "dat_k0", "dat_l0", "dat_n0", "dat_a1", "dat_b1", "dat_c1", "dat_d1", "dat_e1", "dat_f1", "dat_g1", "dat_h1", "dat_i1", "dat_j1", "dat_k1", "dat_l1", "dat_a2", "dat_b2", "dat_c2", "dat_d2", "dat_e2", "dat_f2", "dat_g2", "dat_h2", "dat_i2", "dat_j2", "dat_k2", "dat_a3", "dat_b3", "dat_c3", "dat_d3", "dat_e3", "dat_f3", "dat_g3", "dat_h3", "dat_i3", "dat_j3", "dat_a4", "dat_b4", "dat_c4", "dat_d4", "dat_e4", "dat_f4", "dat_g4", "dat_h4", "dat_i4", "dat_a5", "dat_b5", "dat_c5", "dat_d5", "dat_e5", "dat_f5", "dat_g5", "dat_h5", "dat_a6", "dat_b6", "dat_c6", "dat_d6", "dat_e6", "dat_f6", "dat_g6", "dat_a7", "dat_b7", "dat_c7", "dat_d7", "dat_e7", "dat_a8", "dat_b8", "dat_c8", "dat_d8", "dat_e8", "dat_a9", "dat_b9", "dat_c9", "dat_d9", "dat_a10", "dat_b10", "dat_c10", "dat_a11", "dat_b11", "dat_a12", ######Formulas "formula_a", "formula_b", "formula_c", "formula_d", "formula_e", "formula_f", "formula_g", "formula_h", "formula_i", "formula_j", "formula_k", "formula_l", "formula_m", "formula_n", "formula_a0", "formula_b0", "formula_c0", "formula_d0", "formula_e0", "formula_f0", "formula_g0", "formula_h0", "formula_i0", "formula_j0", "formula_k0", "formula_l0", "formula_n0", "formula_a1", "formula_b1", "formula_c1", "formula_d1", "formula_e1", "formula_f1", "formula_g1", "formula_h1", "formula_i1", "formula_j1", "formula_k1", "formula_l1", "formula_a2", "formula_b2", "formula_c2", "formula_d2", "formula_e2", "formula_f2", "formula_g2", "formula_h2", "formula_i2", "formula_j2", "formula_k2", "formula_a3", "formula_b3", "formula_c3", "formula_d3", "formula_e3", "formula_f3", "formula_g3", "formula_h3", "formula_i3", "formula_j3", "formula_a4", "formula_b4", "formula_c4", "formula_d4", "formula_e4", "formula_f4", "formula_g4", "formula_h4", "formula_i4", "formula_a5", "formula_b5", "formula_c5", "formula_d5", "formula_e5", "formula_f5", "formula_g5", "formula_h5", "formula_a6", "formula_b6", "formula_c6", "formula_d6", "formula_e6", "formula_f6", "formula_g6", "formula_a7", "formula_b7", "formula_c7", "formula_d7", "formula_e7", "formula_a8", "formula_b8", "formula_c8", "formula_d8", "formula_e8", "formula_a9", "formula_b9", "formula_c9", "formula_d9", "formula_a10", "formula_b10", "formula_c10", "formula_a11", "formula_b11", "formula_a12", "formula_a13", "formula_b13", "formula_c13", "formula_d13", "formula_e13", "formula_f13", "formula_g13" )) ###################################################### pdr_a0<-pdredge(mod_a0, cluster=cl) pdr_b0<-pdredge(mod_b0, cluster=cl) pdr_c0<-pdredge(mod_c0, cluster=cl) pdr_d0<-pdredge(mod_d0, cluster=cl) pdr_e0<-pdredge(mod_e0, cluster=cl) pdr_f0<-pdredge(mod_f0, cluster=cl) pdr_g0<-pdredge(mod_g0, cluster=cl) pdr_h0<-pdredge(mod_h0, cluster=cl) pdr_i0<-pdredge(mod_i0, cluster=cl) pdr_j0<-pdredge(mod_j0, cluster=cl) pdr_k0<-pdredge(mod_k0, cluster=cl) pdr_l0<-pdredge(mod_l0, cluster=cl) pdr_n0<-pdredge(mod_l0, cluster=cl) pdr_a1<-pdredge(mod_a1, cluster=cl) pdr_b1<-pdredge(mod_b1, cluster=cl) pdr_c1<-pdredge(mod_c1, cluster=cl) pdr_d1<-pdredge(mod_d1, cluster=cl) pdr_e1<-pdredge(mod_e1, cluster=cl) pdr_f1<-pdredge(mod_f1, cluster=cl) pdr_g1<-pdredge(mod_g1, cluster=cl) pdr_h1<-pdredge(mod_h1, cluster=cl) pdr_i1<-pdredge(mod_i1, cluster=cl) pdr_j1<-pdredge(mod_j1, cluster=cl) pdr_k1<-pdredge(mod_k1, cluster=cl) pdr_l1<-pdredge(mod_l1, cluster=cl) save.image("D:/copyofgdrive/Konstanz/Inventory/EGF_dredge_initial_models_06062017.RData") pdr_a2<-pdredge(mod_a2, cluster=cl) pdr_b2<-pdredge(mod_b2, cluster=cl) pdr_c2<-pdredge(mod_c2, cluster=cl) pdr_d2<-pdredge(mod_d2, cluster=cl) pdr_e2<-pdredge(mod_e2, cluster=cl) pdr_f2<-pdredge(mod_f2, cluster=cl) pdr_g2<-pdredge(mod_g2, cluster=cl) pdr_h2<-pdredge(mod_h2, cluster=cl) pdr_i2<-pdredge(mod_i2, cluster=cl) pdr_j2<-pdredge(mod_j2, cluster=cl) pdr_k2<-pdredge(mod_k2, cluster=cl) pdr_a3<-pdredge(mod_a3, cluster=cl) pdr_b3<-pdredge(mod_b3, cluster=cl) pdr_c3<-pdredge(mod_c3, cluster=cl) pdr_d3<-pdredge(mod_d3, cluster=cl) pdr_e3<-pdredge(mod_e3, cluster=cl) pdr_f3<-pdredge(mod_f3, cluster=cl) pdr_g3<-pdredge(mod_g3, cluster=cl) pdr_h3<-pdredge(mod_h3, cluster=cl) pdr_i3<-pdredge(mod_i3, cluster=cl) pdr_j3<-pdredge(mod_j3, cluster=cl) save.image("D:/copyofgdrive/Konstanz/Inventory/EGF_dredge_initial_models_06062017.RData") pdr_a4<-pdredge(mod_a4, cluster=cl) pdr_b4<-pdredge(mod_b4, cluster=cl) pdr_c4<-pdredge(mod_c4, cluster=cl) pdr_d4<-pdredge(mod_d4, cluster=cl) pdr_e4<-pdredge(mod_e4, cluster=cl) pdr_f4<-pdredge(mod_f4, cluster=cl) pdr_g4<-pdredge(mod_g4, cluster=cl) pdr_h4<-pdredge(mod_h4, cluster=cl) pdr_i4<-pdredge(mod_i4, cluster=cl) pdr_a5<-pdredge(mod_a5, cluster=cl) pdr_b5<-pdredge(mod_b5, cluster=cl) pdr_c5<-pdredge(mod_c5, cluster=cl) pdr_d5<-pdredge(mod_d5, cluster=cl) pdr_e5<-pdredge(mod_e5, cluster=cl) pdr_f5<-pdredge(mod_f5, cluster=cl) pdr_g5<-pdredge(mod_g5, cluster=cl) pdr_h5<-pdredge(mod_h5, cluster=cl) save.image("D:/copyofgdrive/Konstanz/Inventory/EGF_dredge_initial_models_06062017.RData") pdr_a6<-pdredge(mod_a6, cluster=cl) pdr_b6<-pdredge(mod_b6, cluster=cl) pdr_c6<-pdredge(mod_c6, cluster=cl) pdr_d6<-pdredge(mod_d6, cluster=cl) pdr_e6<-pdredge(mod_e6, cluster=cl) pdr_f6<-pdredge(mod_f6, cluster=cl) pdr_g6<-pdredge(mod_g6, cluster=cl) pdr_a7<-pdredge(mod_a7, cluster=cl) pdr_b7<-pdredge(mod_b7, cluster=cl) pdr_c7<-pdredge(mod_c7, cluster=cl) pdr_d7<-pdredge(mod_d7, cluster=cl) pdr_e7<-pdredge(mod_e7, cluster=cl) pdr_a8<-pdredge(mod_a8, cluster=cl) pdr_b8<-pdredge(mod_b8, cluster=cl) pdr_c8<-pdredge(mod_c8, cluster=cl) pdr_d8<-pdredge(mod_d8, cluster=cl) pdr_e8<-pdredge(mod_e8, cluster=cl) save.image("D:/copyofgdrive/Konstanz/Inventory/EGF_dredge_initial_models_06062017.RData") pdr_a9<-pdredge(mod_a9, cluster=cl) pdr_b9<-pdredge(mod_b9, cluster=cl) pdr_c9<-pdredge(mod_c9, cluster=cl) pdr_d9<-pdredge(mod_d9, cluster=cl) pdr_a10<-pdredge(mod_a10, cluster=cl) pdr_b10<-pdredge(mod_b10, cluster=cl) pdr_c10<-pdredge(mod_c10, cluster=cl) pdr_a11<-pdredge(mod_a11, cluster=cl) pdr_b11<-pdredge(mod_b11, cluster=cl) pdr_a12<-pdredge(mod_a12, cluster=cl) pdr_a13<-pdredge(mod_a13, cluster=cl) pdr_b13<-pdredge(mod_b13, cluster=cl) pdr_c13<-pdredge(mod_c13, cluster=cl) pdr_d13<-pdredge(mod_d13, cluster=cl) pdr_e13<-pdredge(mod_e13, cluster=cl) pdr_f13<-pdredge(mod_f13, cluster=cl) pdr_g13<-pdredge(mod_g13, cluster=cl) pdr_a<-pdredge(mod_a, cluster=cl) pdr_b<-pdredge(mod_b, cluster=cl) pdr_c<-pdredge(mod_c, cluster=cl) pdr_d<-pdredge(mod_d, cluster=cl) pdr_e<-pdredge(mod_e, cluster=cl) pdr_f<-pdredge(mod_f, cluster=cl) pdr_g<-pdredge(mod_g, cluster=cl) pdr_h<-pdredge(mod_h, cluster=cl) pdr_i<-pdredge(mod_i, cluster=cl) pdr_j<-pdredge(mod_j, cluster=cl) pdr_k<-pdredge(mod_k, cluster=cl) pdr_l<-pdredge(mod_l, cluster=cl) pdr_m<-pdredge(mod_m, cluster=cl) pdr_n<-pdredge(mod_n, cluster=cl) ############################################### stopCluster(cl) remove(cl) save.image("D:/copyofgdrive/Konstanz/Inventory/EGF_dredge_initial_models_06062017.RData") ############################################# #########Building the saturated model: ##Saturated Model: formula_sat<-Nat_in_E~ # scale(VDV_all)+ scale(curmedian)+ scale(NE_gl_reg)+ scale(Height)+ # scale(Herb)+ # scale(SO)+ scale(Mono)+ # scale(Combo_NRS)+ # scale(PropVeg)+ # scale(HZ)+ # scale(VDV_all)*scale(Height)+ scale(VDV_all)*scale(Mono) # scale(Combo_NRS)*scale(NE_gl_reg)+ # scale(NE_gl_reg)*scale(Height)+ # scale(NE_gl_reg)*scale(PropVeg) # scale(Height)*scale(HZ) #mcmc_data2<-mcmc_data1[c(1:2, 4:15)] #mcmc_data2$Mono[is.na(mcmc_data2$Mono)]<-1 #dat_sat<-na.omit(mcmc_data2) #mod_sat<- uMCMCglmm(fixed = formula_sat, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_sat, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) #mod_sat1<- uMCMCglmm(fixed = formula_sat, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_sat, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) #mod_sat2<- uMCMCglmm(fixed = formula_sat, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_sat, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) #mod_sat3<- uMCMCglmm(fixed = formula_sat, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_sat, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) #mod_sat4<- uMCMCglmm(fixed = formula_sat, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_sat, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) #mod_sat5<- uMCMCglmm(fixed = formula_sat, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_sat, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) #mod_sat6<- uMCMCglmm(fixed = formula_sat, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_sat, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) #mod_sat7<- uMCMCglmm(fixed = formula_sat, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_sat, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) #mod_sat8<- uMCMCglmm(fixed = formula_sat, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_sat, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) #mod_sat9<- uMCMCglmm(fixed = formula_sat, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_sat, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) #mod_sat10<- uMCMCglmm(fixed = formula_sat, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_sat, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) #mod_sat11<- uMCMCglmm(fixed = formula_sat, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_sat, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) #mod_sat12<- uMCMCglmm(fixed = formula_sat, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_sat, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) #mod_sat13<- uMCMCglmm(fixed = formula_sat, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_sat, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) #mod_sat14<- uMCMCglmm(fixed = formula_sat, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_sat, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) #mod_sat15<- uMCMCglmm(fixed = formula_sat, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_sat, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) #mod_sat16<- uMCMCglmm(fixed = formula_sat, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_sat, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) #mod_sat17<- uMCMCglmm(fixed = formula_sat, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_sat, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) mod_sat18<- uMCMCglmm(fixed = formula_sat, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_sat, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) AIC(mod_sat12, mod_sat18) summary(mod_sat12) save.image("D:/copyofgdrive/Konstanz/Inventory/EGF_sat_reduction_nomax_12062017.RData") ############################################# ####Final model: dat_minad<-mcmc_data2[c(1:3, 6:8, 13)] dat_minad<-na.omit(dat_minad) formula_minad<-Nat_in_E~ scale(VDV_all)+ scale(curmedian)+ scale(NE_gl_reg)+ scale(Height)+ scale(Mono)+ scale(VDV_all)*scale(Mono) mod_minad<- uMCMCglmm(fixed = formula_minad, prior=prior1, random= ~ animal,family="categorical",pedigree = mcmc_tree1,scale=F,data = dat_minad, nitt = 110000,burnin = 10000,thin = 10, singular.ok = TRUE) ############################################### #Setting up data for predictions: df3<-df1[c(5, 19,21, 24, 26, 29, 31, 34:46)] yyyz<-df3[which(df3$Species!="Persicaria_virginiana"),] #for some reason there's a problem and this one needs to be removed? mcmc_data3 <- cbind(animal = yyyz$Species, yyyz) dat_pred<-mcmc_data3[c(1:2,6,8,11:13,18 )] dat_pred$Mono[is.na(dat_pred$Mono)]<-1 dat_pred<-na.omit(dat_pred) dat_pred45<-dat_pred colnames(dat_pred45)[3]<-"curmedian" dat_pred85<-dat_pred colnames(dat_pred85)[4]<-"curmedian" dat_predcur<-dat_minad pred45<-predict(mod_minad, newdata=dat_pred45, ###future fits pedigree = mcmc_tree1, marginal=mod_minad$Random$formula, type="response", interval="confidence", level=0.95, it=NULL, posterior="all", verbose=FALSE) save.image("D:/copyofgdrive/Konstanz/Inventory/EGF_predictions_nomax_12062017.RData") pred85<-predict(mod_minad, newdata=dat_pred85, ###future fits pedigree = mcmc_tree1, marginal=mod_minad$Random$formula, type="response", interval="confidence", level=0.95, it=NULL, posterior="all", verbose=FALSE) save.image("D:/copyofgdrive/Konstanz/Inventory/EGF_predictions_nomax_12062017.RData") predcur<-predict(mod_minad, newdata=dat_predcur, ###future fits pedigree = mcmc_tree1, marginal=mod_minad$Random$formula, type="response", interval="confidence", level=0.95, it=NULL, posterior="all", verbose=FALSE) save.image("D:/copyofgdrive/Konstanz/Inventory/EGF_predictions_nomax_12062017.RData") predall<-cbind(dat_minad, pred45) colnames(predall)[8]<-"pred45" colnames(predall)[9]<-"pred45lwr" colnames(predall)[10]<-"pred45upr" predall<-cbind(predall, pred85) colnames(predall)[11]<-"pred85" colnames(predall)[12]<-"pred85lwr" colnames(predall)[13]<-"pred85upr" predall<-cbind(predall, predcur) colnames(predall)[14]<-"predcur" colnames(predall)[15]<-"predcurlwr" colnames(predall)[16]<-"predcurupr" write.csv(predall, "D:/copyofgdrive/Konstanz/Inventory/EGF_predictions_nomax_12062017.csv")