Coral Analysis

The following code generate the parameter table for coral analysis

#install.packages("kable")
library(knitr)
load(url("https://tonyjhwueng.info/ououcir/coralparams.RData"))
kable(round(coralmodelBFreorder,3))

	ououcir	ououbm	oubmcir	oubmbm
ououcir	1.000	3.589	6.119	7.610
ououbm	0.279	1.000	1.705	2.121
oubmcir	0.163	0.587	1.000	1.244
oubmbm	0.131	0.472	0.804	1.000

kable(round(table.out,3))

	alpha.y	alpha.x	alpha.tau	theta.x	theta.tau	sigma.x	tau	sigma.tau	b0	b1	b2
OUBMBM	3.022	NA	NA	NA	NA	0.904	3.490	NA	0.890	1.147	0.179
OUOUBM	2.677	4.430	NA	-0.150	NA	1.241	4.609	NA	0.836	1.216	0.230
OUBMCIR	3.585	NA	3.097	NA	1.577	0.951	NA	1.823	0.899	1.153	0.150
OUOUCIR	2.801	4.073	2.950	-0.121	1.639	1.242	NA	1.884	0.825	1.165	0.253
GLS	NA	NA	NA	NA	NA	NA	NA	NA	0.900	1.206	0.201

The raw code can be accessed in the following (it run for a while, the size of RData is big–500 MB)

rm(list=ls())
library(xtable)
library(abc)
load(url("https://tonyjhwueng.info/ououcir/EmpiricalMaincodeV2/treetraitV2/sanchez.Lasker_coral/sanchez.Lasker_coral_log_thickness.RData"))

summary(modsel.mnlog)

coralogthick<-summary(modsel.mnlog)


modelorder<-c("ououcir","ououbm","oubmcir","oubmbm")
coralmodelBFreorder<-coralogthick$mnlogistic$BayesF[,modelorder]
coralmodelBFreorder<-coralmodelBFreorder[modelorder,]
coralmodelBFreorder



xtable(coralmodelBFreorder,digits=3)



post.oubmbm.mean <- round(apply(post.oubmbm,2,mean),digits = 4)
post.ououbm.mean <- round(apply(post.ououbm,2,mean),digits = 4)
post.oubmcir.mean <- round(apply(post.oubmcir,2,mean),digits = 4)
post.ououcir.mean <- round(apply(post.ououcir,2,mean),digits = 4)


table.params <- data.frame(matrix(NA,4,8))
rownames(table.params) <- c("OUBMBM","OUOUBM","OUBMCIR","OUOUCIR")
colnames(table.params)<- c("alpha.y","sigmasq.x","tau","alpha.x","theta.x","alpha.tau","theta.tau","sigmasq.tau")

table.params["OUBMBM",c("alpha.y","sigmasq.x","tau")]<-post.oubmbm.mean[c("alpha.y","sigmasq.x","tau")]


table.params["OUOUBM",c("alpha.y","alpha.x","theta.x", "sigmasq.x", "tau")]<-post.ououbm.mean[c("alpha.y","alpha.x","theta.x", "sigmasq.x", "tau")] 

table.params["OUBMCIR",c("alpha.y", "sigmasq.x","alpha.tau","theta.tau", "sigmasq.tau")]<-post.oubmcir.mean[c("alpha.y", "sigmasq.x","alpha.tau","theta.tau", "sigmasq.tau")]

table.params["OUOUCIR",c("alpha.y","alpha.x","theta.x","sigmasq.x","alpha.tau","theta.tau", "sigmasq.tau")]<-post.ououcir.mean[c("alpha.y","alpha.x","theta.x","sigmasq.x","alpha.tau","theta.tau", "sigmasq.tau")]
table.params[,"sigmasq.x"]<-sqrt(table.params[,"sigmasq.x"])
table.params[,"sigmasq.tau"]<-sqrt(table.params[,"sigmasq.tau"])

colnames(table.params)<- c("alpha.y","sigma.x","tau","alpha.x","theta.x","alpha.tau","theta.tau","sigma.tau")
table.params<-table.params[, c("alpha.y","alpha.x", "alpha.tau", "theta.x","theta.tau", "sigma.x","tau","sigma.tau")]

xtable(table.params,digits=4)


table.b<- matrix(NA,4,3)
rownames(table.b) <- c("OUBMBM","OUOUBM","OUBMCIR","OUOUCIR")
colnames(table.b)<-c("b0","b1","b2")
table.b[1,]<-post.oubmbm.mean[c(4:6)]
table.b[2,]<-post.ououbm.mean[c(6:8)]
table.b[3,]<-post.oubmcir.mean[c(6:8)]
table.b[4,]<-post.ououcir.mean[c(8:10)]
table.b


ols<-lm(resptrait~predtrait1+predtrait2)
ols$coefficients
OLS<-rep(NA,11)
OLS[9:11]<-ols$coefficients
table.out<-rbind(cbind(table.params,table.b),OLS)
xtable(table.out,digits=3)