#### 1. Getting ready #### rm(list=ls()) getwd() # Use output from this line to modify the next one ... #setwd("PATHWAY TO YOUR WORKING DIRECTORY") Hc_data <- read.table("hypericum_data_94_07.txt", header=T) Hc_pop <- read.table("popmeanHc.txt", header=T) Height_data <-Hc_data\$ht_init[Hc_data\$stage !="sg" & Hc_data\$ht_init<90] #### 2. Plots and calculation of averages using a frequentist approach #### hist(Height_data,100,main="Histogram of Hypericum cumulicola height (cm)",xlab= "Height(cm)") abline(v=Hc_pop\$pop_mean, col="blue") abline(v=mean(Height_data), col="red",lwd=3) mean(Height_data) var(Height_data) round(sqrt(var(Height_data))/sqrt(length(Height_data)),4) summary(lm(Height_data~1)) mean.pop.mean <- mean(Hc_pop\$pop_mean) mean.pop.var <- mean(Hc_pop\$pop_sd^2) par(mfrow=c (1,2)) hist(Hc_pop\$pop_mean,8,main="population mean", xlab="height mean") abline(v = mean(Hc_pop\$pop_mean), col="blue") abline(v=mean(Height_data), col="red",lwd=3) hist(Hc_pop\$pop_sd^2,8,main="population variance",xlab="height variance") abline(v = mean(Hc_pop\$pop_sd^2), col="blue") abline(v=var(Height_data), col="red",lwd=3) #### 3. Sampling from the dataset #### size <-10 x <- sample(Height_data,size) n <- length(x) pop.mean.mean <- mean(Hc_pop\$pop_mean) pop.mean.var <- var(Hc_pop\$pop_mean) log.pop.var.mean <- mean(log(Hc_pop\$pop_sd^2)) log.pop.var.var <- var(log(Hc_pop\$pop_sd^2)) #### 4. Calculation of averages using a Bayesian approach #### library(rjags) library(ggplot2) ## a) Bayesian analysis with uninformed priors ## model = "meanmodel_jags.R" # Inits function inits=list( list(var_height=100,mean_height=100), list(var_height=80,mean_height=60), list(var_height=120,mean_height=90) ) # Bundle data mean.data=list( Y=x, nobs=n ) # Parameters to estimate params <- c("mean_height","var_height") # MCMC settings, start Gibbs sampler and plot results and diagnostics jm=jags.model(model,data=mean.data,inits=inits,n.chains=length(inits),n.adapt=5000) update(jm, n.iter=10000) zc=coda.samples(jm,variable.names=params,n.iter=10000) gelman.diag(zc) plot(zc) summary(zc) R11<- summary(zc) ## b) Bayesian analysis with informed priors ## model = "meanmodel_informative_jags.R" # Bundle data mean.data=list( Y=x, nobs=n, m=pop.mean.mean, v=1/pop.mean.var, mv=log.pop.var.mean, vv=1/log.pop.var.var ) # MCMC settings, start Gibbs sampler and plot results and diagnostics jm=jags.model(model, data=mean.data, inits=inits,n.chains = length(inits), n.adapt=5000) update(jm, n.iter=10000) zc=coda.samples(jm, variable.names=params, n.iter=10000) gelman.diag(zc) plot(zc) summary(zc) R22<- summary(zc) par(mfrow=c(1,2)) hist(rgamma(1000,0.001,0.001),10,xlim=c(0,30),xlab="var height",main ="uninformed") hist(rlnorm(10000,mean(log(Hc_pop\$pop_sd^2)),sqrt(var(log(Hc_pop\$pop_sd^2)))), ylim=c(0,200),400,xlab="var height",main ="informed") #### 5. Comparing all approaches #### Results <- array(0,c(5,6)) colnames(Results) <- c("mean","variance", "M 2.5%", "M 97.5","V 2.5%", "V 97.5") rownames(Results) <- c("Frequentist", "B without priors","B with priors", "Population","Overall") #Frequentist# Results[1,1] <- sample_mean <- mean(x) Results[1,2] <- var(x) sample_std_error <- sd(x)/sqrt(size) Results[1,3] <- sample_mean - qt(0.975, size-1)*sqrt(var(x)/size) Results[1,4] <- sample_mean + qt(0.975, size-1)*sqrt(var(x)/size) Results[1,5] <- NA Results[1,6] <- NA #Bayesian - uninformed priors# Results[2,1] <- R11[1]\$statistics[1] Results[2,2] <- R11[1]\$statistics[2] Results[2,3] <- R11[2]\$quantiles[1] Results[2,4] <- R11[2]\$quantiles[9] Results[2,5] <- R11[2]\$quantiles[2] Results[2,6] <- R11[2]\$quantiles[10] #Bayesian - informed priors# Results[3,1] <- R22[1]\$statistics[1] Results[3,2] <- R22[1]\$statistics[2] Results[3,3] <- R22[2]\$quantiles[1] Results[3,4] <- R22[2]\$quantiles[9] Results[3,5] <- R22[2]\$quantiles[2] Results[3,6] <- R22[2]\$quantiles[10] #Population data# Results[4,1] <- mean(Hc_pop\$pop_mean) Results[4,2] <- mean(Hc_pop\$pop_sd^2) Results[4,3] <- NA Results[4,4] <- NA Results[4,5] <- NA Results[4,6] <- NA #Individual data# Results[5,1] <- mean(Height_data) Results[5,2] <- var(Height_data) Results[5,3] <- mean(Height_data) - qt(0.975, size-1)*sample_std_error Results[5,4] <- mean(Height_data) + qt(0.975, size-1)*sample_std_error Results[5,5] <- NA Results[5,6] <- NA Results par(mfrow=c(1,4)) b=seq(-10000,10000,5) b1=seq(-50,100,1) hist(rnorm(1000,0,sqrt(1e+06)),breaks=b, ylim=c(0,700),xlab="height",main ="uninformed") hist(rnorm(10000,Hc_pop\$pop_mean,sqrt(var(Hc_pop\$pop_mean))),breaks=b1, ylim=c(0,700),xlab="height",main ="informed") hist(Height_data,breaks=b1, ylim=c(0,700),xlab="height",main ="all data") hist(sample(Height_data,10),breaks=b1, ylim=c(0,10),xlab="height",main ="sample of data") ############################################################## s_p <- 100000 prior_m_no <-rnorm(s_p,0, sqrt(1/1.0E-6)) prior_m_in <-rnorm(s_p,mean(Hc_pop\$pop_mean),sqrt(pop.mean.var)) post_m_no <- rnorm(s_p, R11[1]\$statistics[1], R11[1]\$statistics[3]) post_m_in <- rnorm(s_p, R22[1]\$statistics[1], R22[1]\$statistics[3]) 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(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)) 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("informed","posterior_no","posterior_in") qplot(p.graph\$priors,geom="density",colour=p.graph\$type,size=p.graph\$type, xlab="Height")+ scale_color_manual(values = c("red","darkgreen","blue")) + scale_size_manual(values=c(0.7,0.7,0.7))+ geom_vline(xintercept =sample_mean,lwd=2.0) #####################################################################