### ### ### Part I: Preparing and plotting the data ### ### ### rm(list=ls()) getwd() #setwd("pathway") reg_data <- read.table("wg_data.txt", header=T) par(mfrow=c(1,2)) hist(reg_data$lg,xlab="height",main="height") hist(reg_data$lf,xlab="number of tillers",main="number of tillers") par(mfrow=c(1,1)) plot(reg_data$lf,reg_data$lg,pch=16,cex=0.55,xlab="number of tillers",ylab="height") max <- 17 x<- reg_data$lf[reg_data$lf one") abline(h=0, col="red") hist(residuals[reg_data$tp==1],xlab="residuals",main="predicted one") hist(residuals[reg_data$tp==2],xlab="residuals",main="predicted > one") ### ### ### Part III: Bayesian analysis with uninformed priors ### ### ### library(rjags) n <- length(y) # Write model model = "grasses_dif_int.R" # Bundle data win.data=list( x = x, x2 = x2, y = y, tp = tp, n=n ) # Inits function inits <- function(){list(b=c(runif(1),runif(1),runif(1),runif(1),runif(1),runif(1)),prec=100)} # Inits inits=list( list(b=c(runif(1),runif(1),runif(1),runif(1),runif(1),runif(1)),prec=100), list(b=c(runif(1),runif(1),runif(1),runif(1),runif(1),runif(1)),prec=100), list(b=c(runif(1),runif(1),runif(1),runif(1),runif(1),runif(1)),prec=100) ) # Parameters to estimate params <- c("b","prec") # MCMC settings nc=3 ni=10000 nb=100 nt=70000 # MCMC settings, start Gibbs sampler and plot results and diagnostics jm=jags.model(model,data=win.data,inits=inits,n.chains=nc,n.adapt=nb) update(jm, n.iter=ni) zc=coda.samples(jm,variable.names=params,n.iter=nt) # Display results gelman.diag(zc) plot(zc) results<-summary(zc) restat<-results$statistics resquant<-results$quantiles s<-round(runif(100)*500,0) estim1 <- estim2 <- array(0,c(max,100)) for(i in 1:100){ for(j in 1:max){ estim1[j,i] <- as.numeric(zc[s[i],1][1]) + as.numeric(zc[s[i],2][1])*j + as.numeric(zc[s[i],3][1])*j^2+ as.numeric(zc[s[i],4][1]) + as.numeric(zc[s[i],5][1])*j + as.numeric(zc[s[i],6][1])*j^2 estim2[j,i] <- as.numeric(zc[s[i],1][1]) + as.numeric(zc[s[i],2][1])*j + as.numeric(zc[s[i],3][1])*j^2+ as.numeric(zc[s[i],4][1])*2 + as.numeric(zc[s[i],5][1])*j*2 + as.numeric(zc[s[i],6][1])*j^2*2 } } par(mfrow=c(1,1)) a <- seq(1,max,1) plot(reg_data$lf, reg_data$lg,type="n", main="plants per pot",xlab="number of tillers", ylab="height") points(reg_data$lf[reg_data$tp==1], reg_data$lg[reg_data$tp==1],col="darkblue",pch=1,cex=0.55) points(reg_data$lf[reg_data$tp==2], reg_data$lg[reg_data$tp==2],col="lightblue",pch=1,cex=0.55) for (k in 1:100) { lines(a,estim1[,k],col="blue") lines(a,estim2[,k],col="cyan") } lines(xp,preone,col="red") lines(xp,pretwo,col="black") results ### NOTE: It may be necessary to end your session here, and restart a new one to run the next analysis ### ### ### Part IV: Bayesian analysis with informed priors ### ### ### model = "grasses_inf_int.R" # MCMC settings, start Gibbs sampler and plot results and diagnostics jm=jags.model(model,data=win.data,inits=inits,n.chains=nc,n.adapt=nb) update(jm, n.iter=ni) zc=coda.samples(jm,variable.names=params,n.iter=nt) # Display results gelman.diag(zc) plot(zc) results<-summary(zc) restat<-results$statistics resquant<-results$quantiles s<-round(runif(100)*500,0) estim1 <- estim2 <- array(0,c(max,100)) for(i in 1:100){ for(j in 1:max){ estim1[j,i] <- as.numeric(zc[s[i],1][1]) + as.numeric(zc[s[i],2][1])*j + as.numeric(zc[s[i],3][1])*j^2+ as.numeric(zc[s[i],4][1]) + as.numeric(zc[s[i],5][1])*j + as.numeric(zc[s[i],6][1])*j^2 estim2[j,i] <- as.numeric(zc[s[i],1][1]) + as.numeric(zc[s[i],2][1])*j + as.numeric(zc[s[i],3][1])*j^2+ as.numeric(zc[s[i],6][1])*j^2*2 } } par(mfrow=c(1,1)) a <- seq(1,max,1) plot(reg_data$lf, reg_data$lg,type="n", main="plants per pot",xlab="number of tillers", ylab="height") points(reg_data$lf[reg_data$tp==1], reg_data$lg[reg_data$tp==1],col="darkblue",pch=1,cex=0.55) points(reg_data$lf[reg_data$tp==2], reg_data$lg[reg_data$tp==2],col="lightblue",pch=1,cex=0.55) for (k in 1:100) { lines(a,estim1[,k],col="blue") lines(a,estim2[,k],col="cyan") } lines(xp,preone,col="red") lines(xp,pretwo,col="black") results ### ### ### Part V: Compare approaches ### ### ### ## NOTE: the code is currently set up to make the graph for beta2, but can be modified to make those of b0 and b1 library(ggplot2) s_p <- 100000 prior_m_no <-rnorm(s_p,0,sqrt(1/1.0E-6)) prior_m_in <-rnorm(s_p,-0.163,0.019) post_m_no <- rnorm(s_p, 0.0004, 0.13) post_m_in <- rnorm(s_p,-0.266 ,0.078 ) priors <-c(prior_m_no,prior_m_in,post_m_no,post_m_in) length(priors) type <- c(rep(1,s_p),rep(2,s_p),rep(3,s_p),rep(4,s_p)) levels(type) <- c("not_informed","informed","posterior_no","posterior_in") length(type) p.graph <- cbind(type,priors) p.graph <-as.data.frame(p.graph) p.graph$type <- factor(p.graph$type) qplot(p.graph$priors,geom="density",colour=p.graph$type,size=p.graph$type)+ scale_color_manual(values = c("black","red","darkgreen","blue")) + scale_size_manual(values=c(0.7,0.7,0.7,0.7)) priors <-c(post_m_no,post_m_in) length(priors) type <- c(rep(1,s_p),rep(2,s_p)) length(type) p.graph <- cbind(type,priors) p.graph <-as.data.frame(p.graph) p.graph$type <- factor(p.graph$type) levels(p.graph$type) <- c("posterior_no","posterior_in") qplot(p.graph$priors,geom="density",colour=p.graph$type,size=p.graph$type, xlab="beta 2")+ scale_color_manual(values = c("red","darkgreen","blue")) + scale_size_manual(values=c(0.7,0.7,0.7))+ geom_vline(xintercept = -0.017,lwd=2.0)