library(phytools) library(geiger) ############################################## ############################################## ######## ######## Arguments: ######## nsim - class numeric; number of simulations to run ######## ######## test char - class character, either "total", "premax" or "max". Determines ######## whether the model fitting is carried out on a character ######## representing number of teeth in the whole tooth row, the ######## premaxillary tooth row or maxillary tooth row ######## ######## All other arguments as for the tooth.sim() function ######## ############################################## ############################################## multi.tooth.sim<-function(nsim=100,test.char="total",max.teeth=40,param.H.1=0.002,param.H.2=0.002,param.H.3=0.002,param.T=0.5,param.R=0.04) { res.matrix<-matrix(ncol=4,nrow=nsim) colnames(res.matrix)<-c("ER AICc","Meristic AICc","ER weights","Meristic weights") for(i in 1:nsim) { ######simulate a phylogeny with the required number of taxa suit.sim<-F while(suit.sim==F) { tree<-pbtree(b=0.1,d=0.05,n=100,scale=100) if(length(tree$tip.label)>=50) { suit.sim<-T } } #####simulate a tooth row. tooth.reg<-tooth.sim(tree,max.teeth=max.teeth,param.H.1=param.H.1,param.H.2=param.H.2,param.H.3=param.H.3,param.T=param.T,param.R=param.R) ######find number of teeth in region of interest if(test.char=="total") { if(is.vector(tooth.reg)) { tot.tooth<-tooth.reg } else(tot.tooth<-rowSums(tooth.reg))#### calculate total number of teeth } else if(test.char=="premax") { tot.tooth<-tooth.reg[,1]#### calculate number of teeth in premaxilla } if(test.char=="max") { tot.tooth<-rowSums(tooth.reg[,2:4])#### calculate number of teeth in maxilla } #####convert to a character with five states, or if less than five teeth #####in the region of interest, as many states as there are teeth if(length(unique(tot.tooth))>4) { diff<-max(tot.tooth)-min(tot.tooth) gap<-round(diff/4) tooth.char<-round(tot.tooth/gap) } else(tooth.char<-tot.tooth) tooth.char<-as.character(tooth.char) names(tooth.char)<-names(tot.tooth) ##### model fitting ER<-fitDiscrete(tree,tooth.char,model="ER")###equal rates ########ordered q matrix q.mat<-matrix(nrow=length(unique(tooth.char)),ncol=length(unique(tooth.char)),data=0) for(j in 1:ncol(q.mat)) { if(j!=1) { q.mat[j-1,j]<-1 } if(j!=ncol(q.mat)) { q.mat[j+1,j]<-1 } q.mat[j,i]<-0 } rownames(q.mat)<-colnames(q.mat)<-1:length(unique(tooth.char)) merist<-fitDiscrete(tree,tooth.char,model=q.mat) ####ordered ####extract AIC scores and weights res.matrix[i,1]<-ER$opt$aicc res.matrix[i,2]<-merist$opt$aicc res.matrix[i,3:4]<-aicw(res.matrix[i,1:2])[,3] print(i) } return(res.matrix) }