#CROPS545_Lecture29_Bayes_alphabet (A/B and Cpi) #Zhiwu Zhang #April 25, 2016 rm(list=ls()) #Import GAPIT #source("http://www.bioconductor.org/biocLite.R") #biocLite("multtest") #install.packages("EMMREML") #install.packages("gplots") #install.packages("scatterplot3d") library('MASS') # required for ginv library(multtest) library(gplots) library(compiler) #required for cmpfun library("scatterplot3d") library("EMMREML") source("http://zzlab.net/GAPIT/emma.txt") source("http://zzlab.net/GAPIT/emma.txt") source("http://zzlab.net/GAPIT/gapit_functions.txt") #Prepare BAGS source('http://zzlab.net/StaGen/R/BAGS.R') #source('~/Dropbox/Current/ZZLab/WSUCourse/CROPS545/Demo/BAGS.R') #Import demo data myGD=read.table(file="http://zzlab.net/GAPIT/data/mdp_numeric.txt",head=T) myGM=read.table(file="http://zzlab.net/GAPIT/data/mdp_SNP_information.txt",head=T) myCV=read.table(file="http://zzlab.net/GAPIT/data/mdp_env.txt",head=T) #Preparing data X=myGD[,-1] taxa=myGD[,1] index1to5=myGM[,2]<6 X1to5 = X[,index1to5] GD.candidate=cbind(as.data.frame(taxa),X1to5) set.seed(99164) mySim=GAPIT.Phenotype.Simulation(GD=GD.candidate,GM=myGM[index1to5,],h2=.5,NQTN=100, effectunit =.95,QTNDist="normal",CV=myCV,cveff=c(.0002,.0002),a2=.5,adim=3,category=1,r=.4) n=nrow(X) m=ncol(X) setwd("~/Desktop/temp") #Change the directory to yours set.seed(99164) ref=sample(n,round(n/2),replace=F) GR=myGD[ref,-1];YR=as.matrix(mySim$Y[ref,2]) GI=myGD[-ref,-1];YI=as.matrix(mySim$Y[-ref,2]) #running BAGS source('~/Dropbox/Current/ZZLab/WSUCourse/CROPS545/Demo/BAGS.R') burn.in=100;burn.out=100;niter=burn.in+burn.out set.seed(99164) myBayes=BAGS(X=GR,y=YR,pi=0,burn.in=100,burn.out=100,recording=T) str(myBayes) #predicting YRP = as.matrix(GR) %*% myBayes$effect #reference YIP = as.matrix(GI) %*% myBayes$effect #inference ### step 7: calculating accuracy fit = cor(YR,YRP);accuracy = cor(YI,YIP) fit;accuracy par(mfrow=c(2,2), mar = c(3,4,1,1)) plot(myBayes$mcmc.p[,1],type="b") plot(myBayes$mcmc.p[,2],type="b") plot(myBayes$mcmc.p[,3],type="b") plot(myBayes$mcmc.p[,4],type="b") ylim=c(min(myBayes$mcmc.b),max(myBayes$mcmc.b)) par(mfrow=c(1,1), mar = c(3,4,1,1)) plot(myBayes$mcmc.b[,1],type="p",ylim=c(-.03,.03)) for(i in 2:m){ points(myBayes$mcmc.b[,i],type="p",col=i) } myVar=myBayes$mcmc.v av=myVar for (j in 1:m){ for(i in 1:niter){ av[i,j]=mean(myVar[1:i,j]) }} ylim=c(min(av),max(av)) plot(av[,1],type="l",ylim=c(0,.005)) for(i in 2:m){ points(av[,i],type="l",col=i) }