library(rstan) library(foreign) library(ape) rstan_options(auto_write = TRUE) #options(width=Sys.getenv("COLUMNS")) # fix the number of columns # a value used for missing data -- when this occurs, # stan knows that the data point must be imputed sentinel = -9999; t <- read.nexus("WSpeciesTree.nex") t$edge.length = t$edge.length / max(t$edge.length) phylocov=vcv(t) d = read.csv('NumberData.csv', header = TRUE) d = subset(d, species <= 33) #for this model reads predictor file only for species number coding s <- read.csv('SpeciesPredictors.csv') m <- stan_model('model1.stan') #added to prevent error 'recompiling to avoid crashing R' SW <- NULL # species W data frame D <- NULL # mycol="rbv.no" # Construct the data to send data <- list( sentinel=sentinel, N_SPECIES=length(unique(d$species)), N_STUDIES=length(unique(d$study)), N_SUBJECTS=length(unique(d$subject)), N_ROWS=nrow(d), N_TASKS=3, subject=d$subject, species=d$species, study=d$study, task=d$task, n1=d$n1, n2=d$n2, ncorrect=as.integer(d$ncorrect), ntrials=d$ntrials, phylocov=phylocov ) fit <- sampling(m, data=data, iter=10000, thin=5, chains=5, cores=8, pars=c("subject_effect"), include=FALSE, control=list(adapt_delta=0.999, max_treedepth=15)) print("----------------------------------------------------------") print(fit) saveRDS(fit, file = "fitmodel1.rds") # main predictors of interest here b0v <- extract(fit, "B0")$B0 lv <- extract(fit, "lambda_w")$lambda_w # extract all of the species samples for(s in 1:length(unique(d$species))) { SW <- rbind(SW, data.frame(species=s, value=extract(fit, paste0("species_W[",s,"]"))$species_W )) } D <- rbind(D, data.frame(b0.5=quantile(b0v,0.05), b0.25=quantile(b0v,0.25), b0.50=quantile(b0v,.5), b0.75=quantile(b0v,0.75), b0.95=quantile(b0v,0.95), l.5=quantile(lv,0.05), l.25=quantile(lv,0.25), l.50=quantile(lv,.5), l.75=quantile(lv,0.75), l.95=quantile(lv,0.95))) write.csv(D, "D.csv", row.names=TRUE) D <- rstan::summary(fit) write.csv(D,"summaryfitmodel1.csv") #Plot the species Ws (no predictors) #plotting in order of phylogenetic tree, node numbers determined by ggtree package SW$species <- as.factor(SW$species) m <- ggplot(SW, aes(x=species, y=value)) + geom_boxplot(outlier.shape=NA) + scale_x_discrete(limits = c("5","1","33","4","2","3","30","31", "12","13","11","6","9","10","8", "7","32","14","15","16","17","18","19","20", "21","22","23","24","25","29","26","27","28")) + scale_y_continuous(limits = c(0,4)) + coord_flip() + theme_classic() png("speciesW.png") print(m) dev.off()