#' summary_single
#'
#' Summarize output of single trajectory model
#'
#' @param model: List, output from dualtraj or dualtrajMS
#' @param X: Matrix, design matrix for series 1. 1st column should be the id.
#' @param y: Vector, outcomes for series 1
#' @param z: Matrix, K x dim(X)[2] indicator matrix indicating which variables to inlcude in each group.
#' @param burn: float, fraction of draws to keep for post burn-in period.
#'
#' @export
summary_single = function(model,X,y,z,burn) {
#number of observations and covariates
n = dim(model$c)[1]
K = length(model$beta)
#initialize dataframe
df = data.frame(A= numeric(0), B= numeric(0), C= numeric(0), D= numeric(0),E= numeric(0))
#Get posterior quanttities for beta
beta.e = matrix(0,nrow=K,ncol=max(do.call(rbind,lapply(model$beta,dim))[,2]))
for (k in 1:K) {
A = cbind(colMeans(tail(model$beta[[k]],n*burn)),
apply(tail(model$beta[[k]],n*burn), 2, sd),
t(apply(tail(model$beta[[k]],n*burn),2,quantile,probs=c(0.025,0.5,0.975))))
beta.e[k,z[k,]==1] = A[,1]
rownames(A) = sprintf(paste("beta_",k,"[%d]",sep=''),seq(1:dim(model$beta[[k]])[2]))
df = rbind(df,A)
}
#Get posterior quantities for standard deviation (NOTE THIS IS NOT THE VARIANCE!!)
A = as.data.frame(t(c(mean(tail(sqrt(model$sigma),n*burn)),
sd(tail(sqrt(model$sigma),n*burn)),
quantile(tail(sqrt(model$sigma),n*burn),probs=c(0.025,0.5,0.975)))))
sigma.e = as.numeric(A[1])^2
rownames(A) = "sigma"
df = rbind(df,A)
#Get posterior quantities for group memberships
A = cbind(colMeans(tail(100*model$pi,n*burn)),
apply(tail(100*model$pi,n*burn),2,sd),
t(apply(tail(100*model$pi,n*burn),2,quantile,probs=c(0.025,0.5,0.975))))
pi.e = A[,1]/100
rownames(A) = sprintf("pi[%d]",seq(1:K))
df = rbind(df,A)
#relabel dataframe
colnames(df)[1:2] = c("Estimate","Standard Deviation")
#Get likelihood and BIC at posterior median
id = X[,1]
X[,1]=1
ll = log_lik(X,y,pi,beta.e,sigma.e,id)
BIC = BIC_calc(X,y,pi,beta.e,sigma.e,id,z)
return(list(estimates = df,
log.likelihood = ll,
BIC = BIC))
}
#' summary_single_MS
#'
#' Summarize output of single trajectory model with Bayesian model averaging.
#'
#' @param model: List, output from dualtraj or dualtrajMS
#' @param X: Matrix, design matrix for series 1. 1st column should be the id.
#' @param y: Vector, outcomes for series 1
#' @param burn: float, fraction of draws to keep for post burn-in period.
#'
#' @export
summary_single_MS = function(model,X,y,burn) {
#number of observations and covariates
n = dim(model$c)[1]
K = length(model$beta)
#initialize dataframe
df = data.frame(A= numeric(0), B= numeric(0), C= numeric(0), D= numeric(0),E= numeric(0),F= numeric(0))
beta.e = matrix(0,nrow=K,ncol=dim(X)[2])
sigma.e = rep(NA,K)
for (k in 1:K) {
#Get posterior quanttities for beta
A = cbind(colMeans(tail(model$beta[[k]],n*burn)),
apply(tail(model$beta[[k]],n*burn), 2, sd),
t(apply(tail(model$beta[[k]],n*burn),2,quantile,probs=c(0.025,0.5,0.975))),
colMeans(model$z[[k]]))
beta.e[k,] = A[,4]
rownames(A) = sprintf(paste("beta_",k,"[%d]",sep=''),seq(1:dim(model$beta[[k]])[2]))
df = rbind(df,A)
#Get posterior quantities for standard deviation (NOTE THIS IS NOT THE VARIANCE!!)
A = as.data.frame(t(c(mean(tail(sqrt(model$sigma[,k]),n*burn)),
sd(tail(sqrt(model$sigma[,k]),n*burn)),
quantile(tail(sqrt(model$sigma[,k]),n*burn),probs=c(0.025,0.5,0.975)),
NA)))
sigma.e[k] = as.numeric(A[1])^2
rownames(A) = paste("sigma_",k,sep='')
df = rbind(df,A)
}
#Get posterior quantities for group memberships
A = cbind(colMeans(tail(100*model$pi,n*burn)),
apply(tail(100*model$pi,n*burn),2,sd),
t(apply(tail(100*model$pi,n*burn),2,quantile,probs=c(0.025,0.5,0.975))),
NA)
pi1.e = A[,1]/100
rownames(A) = sprintf("pi[%d]",seq(1:K))
df = rbind(df,A)
#relabel dataframe
colnames(df)[1:2] = c("Estimate","Standard Deviation")
colnames(df)[6] = "Inclusion Prob."
#Get likelihood and BIC at posterior median
id = X[,1]
X[,1]=1
z = (beta.e != 0)*1
ll = log_lik(X,y,pi,beta.e,sigma.e,id)
BIC = BIC_calc(X,y,pi,beta.e,sigma.e,id,z)
return(list(estimates = df,
log.likelihood = ll,
BIC = BIC))
}
#' summary_dual
#'
#' Summarize output of dual trajectory model
#'
#' @param model: List, output from dualtraj or dualtrajMS
#' @param X1: Matrix, design matrix for series 1. 1st column should be the id.
#' @param X2: Matrix, design matrix for series 2. 1st column should be the id.
#' @param y1: Vector, outcomes for series 1
#' @param y2: Vector, outcomes for series 2
#' @param z1: Matrix, K1 x dim(X1)[2] indicator matrix indicating which variables to inlcude in each group.
#' @param z2: Matrix, K2 x dim(X2)[2] indicator matrix indicating which variables to inlcude in each group.
#' @param burn: float, fraction of draws to keep for post burn-in period.
#'
#' @export
summary_dual = function(model,X1,X2,y1,y2,z1,z2,burn) {
#number of observations and covariates
n1 = dim(model$c1)[1]
K1 = length(model$beta1)
K2 = length(model$beta2)
#initialize dataframe
df = data.frame(A= numeric(0), B= numeric(0), C= numeric(0), D= numeric(0),E= numeric(0))
#group 1
#Get posterior quanttities for beta
beta1.e = matrix(0,nrow=K1,ncol=max(do.call(rbind,lapply(model$beta1,dim))[,2]))
for (k in 1:K1) {
A = cbind(colMeans(tail(model$beta1[[k]],n1*burn)),
apply(tail(model$beta1[[k]],n1*burn), 2, sd),
t(apply(tail(model$beta1[[k]],n1*burn),2,quantile,probs=c(0.025,0.5,0.975))))
beta1.e[k,z1[k,]==1] = A[,1]
rownames(A) = sprintf(paste("beta1_",k,"[%d]",sep=''),seq(1:dim(model$beta1[[k]])[2]))
df = rbind(df,A)
}
#Get posterior quantities for standard deviation (NOTE THIS IS NOT THE VARIANCE!!)
A = as.data.frame(t(c(mean(tail(sqrt(model$sigma1),n1*burn)),
sd(tail(sqrt(model$sigma1),n1*burn)),
quantile(tail(sqrt(model$sigma1),n1*burn),probs=c(0.025,0.5,0.975)))))
sigma1.e = as.numeric(A[1])^2
rownames(A) = "sigma1"
df = rbind(df,A)
#Get posterior quantities for group memberships
A = cbind(colMeans(tail(100*model$pi1,n1*burn)),
apply(tail(100*model$pi1,n1*burn),2,sd),
t(apply(tail(100*model$pi1,n1*burn),2,quantile,probs=c(0.025,0.5,0.975))))
pi1.e = A[,1]/100
rownames(A) = sprintf("pi1[%d]",seq(1:K1))
df = rbind(df,A)
#group 2
#Get posterior quanttities for beta
beta2.e = matrix(0,nrow=K2,ncol=max(do.call(rbind,lapply(model$beta2,dim))[,2]))
for (k in 1:K2) {
A = cbind(colMeans(tail(model$beta2[[k]],n1*burn)),
apply(tail(model$beta2[[k]],n1*burn), 2, sd),
t(apply(tail(model$beta2[[k]],n1*burn),2,quantile,probs=c(0.025,0.5,0.975))))
beta2.e[k,z2[k,]==1] = A[,1]
rownames(A) = sprintf(paste("beta2_",k,"[%d]",sep=''),seq(1:dim(model$beta2[[k]])[2]))
df = rbind(df,A)
}
#Get posterior quantities for standard deviation (NOTE THIS IS NOT THE VARIANCE!!)
A = as.data.frame(t(c(mean(tail(sqrt(model$sigma2),n1*burn)),
sd(tail(sqrt(model$sigma2),n1*burn)),
quantile(tail(sqrt(model$sigma2),n1*burn),probs=c(0.025,0.5,0.975)))))
sigma2.e = as.numeric(A[1])^2
rownames(A) = "sigma2"
df = rbind(df,A)
#Get posterior quantities for group memberships
A = cbind(colMeans(tail(100*model$pi2,n1*burn)),
apply(tail(100*model$pi2,n1*burn),2,sd),
t(apply(tail(100*model$pi2,n1*burn),2,quantile,probs=c(0.025,0.5,0.975))))
rownames(A) = sprintf("pi2[%d]",seq(1:K2))
df = rbind(df,A)
#transitions
A = cbind(as.vector(t(apply(100*model$pi1_2[(n1*(1-burn)+1):n1,,], c(2,3), mean))),
as.vector(t(apply(100*model$pi1_2[(n1*(1-burn)+1):n1,,], c(2,3), sd))),
as.vector(t(apply(100*model$pi1_2[(n1*(1-burn)+1):n1,,], c(2,3), quantile,probs=0.025))),
as.vector(t(apply(100*model$pi1_2[(n1*(1-burn)+1):n1,,], c(2,3), quantile,probs=0.5))),
as.vector(t(apply(100*model$pi1_2[(n1*(1-burn)+1):n1,,], c(2,3), quantile,probs=0.975))))
pi1_2.e = matrix(A[,1],nrow=K1,ncol=K2,byrow=TRUE)/100
colnames(A) = colnames(df)
rn = c()
for (k1 in 1:K1) {
for (k2 in 1:K2) {
rn = c(rn,paste("1->2,",k1,"->",k2,sep=''))
}
}
rownames(A) = rn
df = rbind(df,A)
#transitions
A = cbind(as.vector(apply(100*model$pi2_1[(n1*(1-burn)+1):n1,,], c(2,3), mean)),
as.vector(apply(100*model$pi2_1[(n1*(1-burn)+1):n1,,], c(2,3), sd)),
as.vector(apply(100*model$pi2_1[(n1*(1-burn)+1):n1,,], c(2,3), quantile,probs=0.025)),
as.vector(apply(100*model$pi2_1[(n1*(1-burn)+1):n1,,], c(2,3), quantile,probs=0.5)),
as.vector(apply(100*model$pi2_1[(n1*(1-burn)+1):n1,,], c(2,3), quantile,probs=0.975)))
colnames(A) = colnames(df)
rn = c()
for (k2 in 1:K2) {
for (k1 in 1:K1) {
rn = c(rn,paste("2->1,",k2,"->",k1,sep=''))
}
}
rownames(A) = rn
df = rbind(df,A)
#joint
A = cbind(as.vector(t(apply(100*model$pi12[(n1*(1-burn)+1):n1,,], c(2,3), mean))),
as.vector(t(apply(100*model$pi12[(n1*(1-burn)+1):n1,,], c(2,3), sd))),
as.vector(t(apply(100*model$pi12[(n1*(1-burn)+1):n1,,], c(2,3), quantile,probs=0.025))),
as.vector(t(apply(100*model$pi12[(n1*(1-burn)+1):n1,,], c(2,3), quantile,probs=0.5))),
as.vector(t(apply(100*model$pi12[(n1*(1-burn)+1):n1,,], c(2,3), quantile,probs=0.975))))
colnames(A) = colnames(df)
rn = c()
for (k1 in 1:K1) {
for (k2 in 1:K2) {
rn = c(rn,paste(k1,k2,sep=','))
}
}
rownames(A) = rn
df = rbind(df,A)
#relabel dataframe
colnames(df)[1:2] = c("Estimate","Standard Deviation")
#Get likelihood and BIC at posterior median
id1 = X1[,1]
id2 = X2[,1]
X1[,1]=1
X2[,1]=1
ll = log_lik_dual(X1,X2,y1,y2,pi1.e,pi1_2.e,beta1.e,beta2.e,sigma1.e,sigma2.e,id1,id2)
BIC = BIC_calc_dual(X1,X2,y1,y2,pi1.e,pi1_2.e,beta1.e,beta2.e,sigma1.e,sigma2.e,id1,id2,z1,z2)
return(list(estimates = df,
log.likelihood = ll,
BIC = BIC))
}
#' summary_dual_constrained
#'
#' Summarize output of dual trajectory model
#'
#' @param model: List, output from dualtraj or dualtrajMS
#' @param X1: Matrix, design matrix for series 1. 1st column should be the id.
#' @param X2: Matrix, design matrix for series 2. 1st column should be the id.
#' @param y1: Vector, outcomes for series 1
#' @param y2: Vector, outcomes for series 2
#' @param z: Matrix, K x dim(X1)[2] indicator matrix indicating which variables to inlcude
#' @param burn: float, fraction of draws to keep for post burn-in period.
#'
#' @export
summary_dual_constrained = function(model,X1,X2,y1,y2,z,burn) {
n = dim(model$c1)[1]
K = length(model$beta)
df = data.frame(A= numeric(0), B= numeric(0), C= numeric(0), D= numeric(0),E= numeric(0))
#group 1
beta.e = matrix(0,nrow=K,ncol=max(do.call(rbind,lapply(model$beta,dim))[,2]))
for (k in 1:K) {
A = cbind(colMeans(tail(model$beta[[k]],n*burn)),
apply(tail(model$beta[[k]],n*burn), 2, sd),
t(apply(tail(model$beta[[k]],n*burn),2,quantile,probs=c(0.025,0.5,0.975))))
beta.e[k,z[k,]==1] = A[,1]
rownames(A) = sprintf(paste("beta_",k,"[%d]",sep=''),seq(1:dim(model$beta[[k]])[2]))
df = rbind(df,A)
}
A = as.data.frame(t(c(mean(tail(sqrt(model$sigma),n*burn)),
sd(tail(sqrt(model$sigma),n*burn)),
quantile(tail(sqrt(model$sigma),n*burn),probs=c(0.025,0.5,0.975)))))
sigma.e = as.numeric(A[1])^2
rownames(A) = "sigma"
df = rbind(df,A)
A = cbind(colMeans(tail(100*model$pi1,n*burn)),
apply(tail(100*model$pi1,n*burn),2,sd),
t(apply(tail(100*model$pi1,n*burn),2,quantile,probs=c(0.025,0.5,0.975))))
pi1.e = A[,1]/100
rownames(A) = sprintf("pi1[%d]",seq(1:K))
df = rbind(df,A)
#transitions
A = cbind(as.vector(t(apply(100*model$pi1_2[(n*(1-burn)+1):n,,], c(2,3), mean))),
as.vector(t(apply(100*model$pi1_2[(n*(1-burn)+1):n,,], c(2,3), sd))),
as.vector(t(apply(100*model$pi1_2[(n*(1-burn)+1):n,,], c(2,3), quantile,probs=0.025))),
as.vector(t(apply(100*model$pi1_2[(n*(1-burn)+1):n,,], c(2,3), quantile,probs=0.5))),
as.vector(t(apply(100*model$pi1_2[(n*(1-burn)+1):n,,], c(2,3), quantile,probs=0.975))))
pi1_2.e = matrix(A[,1],nrow=K,ncol=K,byrow=TRUE)/100
colnames(A) = colnames(df)
rn = c()
for (k1 in 1:K) {
for (k2 in 1:K) {
rn = c(rn,paste("1->2,",k1,"->",k2,sep=''))
}
}
rownames(A) = rn
df = rbind(df,A)
#transitions
A = cbind(as.vector(apply(100*model$pi2_1[(n*(1-burn)+1):n,,], c(2,3), mean)),
as.vector(apply(100*model$pi2_1[(n*(1-burn)+1):n,,], c(2,3), sd)),
as.vector(apply(100*model$pi2_1[(n*(1-burn)+1):n,,], c(2,3), quantile,probs=0.025)),
as.vector(apply(100*model$pi2_1[(n*(1-burn)+1):n,,], c(2,3), quantile,probs=0.5)),
as.vector(apply(100*model$pi2_1[(n*(1-burn)+1):n,,], c(2,3), quantile,probs=0.975)))
colnames(A) = colnames(df)
rn = c()
for (k2 in 1:K) {
for (k1 in 1:K) {
rn = c(rn,paste("2->1,",k2,"->",k1,sep=''))
}
}
rownames(A) = rn
df = rbind(df,A)
#joint
A = cbind(as.vector(t(apply(100*model$pi12[(n*(1-burn)+1):n,,], c(2,3), mean))),
as.vector(t(apply(100*model$pi12[(n*(1-burn)+1):n,,], c(2,3), sd))),
as.vector(t(apply(100*model$pi12[(n*(1-burn)+1):n,,], c(2,3), quantile,probs=0.025))),
as.vector(t(apply(100*model$pi12[(n*(1-burn)+1):n,,], c(2,3), quantile,probs=0.5))),
as.vector(t(apply(100*model$pi12[(n*(1-burn)+1):n,,], c(2,3), quantile,probs=0.975))))
colnames(A) = colnames(df)
rn = c()
for (k1 in 1:K) {
for (k2 in 1:K) {
rn = c(rn,paste(k1,k2,sep=','))
}
}
rownames(A) = rn
df = rbind(df,A)
colnames(df)[1:2] = c("Estimate","Standard Deviation")
id1 = X1[,1]
id2 = X2[,1]
X1[,1]=1
X2[,1]=1
ll = log_lik_dual(X1,X2,y1,y2,pi1.e,pi1_2.e,beta.e,beta.e,sigma.e,sigma.e,id1,id2)
BIC = BIC_calc_dual(X1,X2,y1,y2,pi1.e,pi1_2.e,beta.e,beta.e,sigma.e,sigma.e,id1,id2,z,z,constrain=TRUE)
return(list(estimates = df,
log.likelihood = ll,
BIC = BIC))
}
#' summary_dual_MS
#'
#' Summarize output of dual trajectory model with Bayesian model averaging.
#'
#' @param model: List, output from dualtraj or dualtrajMS
#' @param X1: Matrix, design matrix for series 1. 1st column should be the id.
#' @param X2: Matrix, design matrix for series 2. 1st column should be the id.
#' @param y1: Vector, outcomes for series 1
#' @param y2: Vector, outcomes for series 2
#' @param burn: float, fraction of draws to keep for post burn-in period.
#'
#' @export
summary_dual_MS = function(model,X1,X2,y1,y2,burn) {
#number of observations and covariates
n1 = dim(model$c1)[1]
K1 = length(model$beta1)
K2 = length(model$beta2)
#initialize dataframe
df = data.frame(A= numeric(0), B= numeric(0), C= numeric(0), D= numeric(0),E= numeric(0),F= numeric(0))
#group 1
beta1.e = matrix(0,nrow=K1,ncol=dim(X1)[2])
sigma1.e = rep(NA,K1)
for (k in 1:K1) {
#Get posterior quanttities for beta
A = cbind(colMeans(tail(model$beta1[[k]],n1*burn)),
apply(tail(model$beta1[[k]],n1*burn), 2, sd),
t(apply(tail(model$beta1[[k]],n1*burn),2,quantile,probs=c(0.025,0.5,0.975))),
colMeans(model$z1[[k]]))
beta1.e[k,] = A[,4]
rownames(A) = sprintf(paste("beta1_",k,"[%d]",sep=''),seq(1:dim(model$beta1[[k]])[2]))
df = rbind(df,A)
#Get posterior quantities for standard deviation (NOTE THIS IS NOT THE VARIANCE!!)
A = as.data.frame(t(c(mean(tail(sqrt(model$sigma1[,k]),n1*burn)),
sd(tail(sqrt(model$sigma1[,k]),n1*burn)),
quantile(tail(sqrt(model$sigma1[,k]),n1*burn),probs=c(0.025,0.5,0.975)),
NA)))
sigma1.e[k] = as.numeric(A[1])^2
rownames(A) = paste("sigma1_",k,sep='')
df = rbind(df,A)
}
#Get posterior quantities for group memberships
A = cbind(colMeans(tail(100*model$pi1,n1*burn)),
apply(tail(100*model$pi1,n1*burn),2,sd),
t(apply(tail(100*model$pi1,n1*burn),2,quantile,probs=c(0.025,0.5,0.975))),
NA)
pi1.e = A[,1]/100
rownames(A) = sprintf("pi1[%d]",seq(1:K1))
df = rbind(df,A)
#group 2
beta2.e = matrix(0,nrow=K2,ncol=dim(X2)[2])
sigma2.e = rep(NA,K2)
for (k in 1:K2) {
#Get posterior quanttities for beta
A = cbind(colMeans(tail(model$beta2[[k]],n1*burn)),
apply(tail(model$beta2[[k]],n1*burn), 2, sd),
t(apply(tail(model$beta2[[k]],n1*burn),2,quantile,probs=c(0.025,0.5,0.975))),
colMeans(model$z2[[k]]))
beta2.e[k,] = A[,4]
rownames(A) = sprintf(paste("beta2_",k,"[%d]",sep=''),seq(1:dim(model$beta2[[k]])[2]))
df = rbind(df,A)
#Get posterior quantities for standard deviation (NOTE THIS IS NOT THE VARIANCE!!)
A = as.data.frame(t(c(mean(tail(sqrt(model$sigma2[,k]),n1*burn)),
sd(tail(sqrt(model$sigma2[,k]),n1*burn)),
quantile(tail(sqrt(model$sigma2[,k]),n1*burn),probs=c(0.025,0.5,0.975)),
NA)))
sigma2.e[k] = as.numeric(A[1])^2
rownames(A) = paste("sigma2_",k,sep='')
df = rbind(df,A)
}
#Get posterior quantities for group memberships
A = cbind(colMeans(tail(100*model$pi2,n1*burn)),
apply(tail(100*model$pi2,n1*burn),2,sd),
t(apply(tail(100*model$pi2,n1*burn),2,quantile,probs=c(0.025,0.5,0.975))),
NA)
rownames(A) = sprintf("pi2[%d]",seq(1:K2))
df = rbind(df,A)
#transitions
A = cbind(as.vector(t(apply(100*model$pi1_2[(n1*(1-burn)+1):n1,,], c(2,3), mean))),
as.vector(t(apply(100*model$pi1_2[(n1*(1-burn)+1):n1,,], c(2,3), sd))),
as.vector(t(apply(100*model$pi1_2[(n1*(1-burn)+1):n1,,], c(2,3), quantile,probs=0.025))),
as.vector(t(apply(100*model$pi1_2[(n1*(1-burn)+1):n1,,], c(2,3), quantile,probs=0.5))),
as.vector(t(apply(100*model$pi1_2[(n1*(1-burn)+1):n1,,], c(2,3), quantile,probs=0.975))),
NA)
pi1_2.e = matrix(A[,1],nrow=K1,ncol=K2,byrow=TRUE)/100
colnames(A) = colnames(df)
rn = c()
for (k1 in 1:K1) {
for (k2 in 1:K2) {
rn = c(rn,paste("1->2,",k1,"->",k2,sep=''))
}
}
rownames(A) = rn
df = rbind(df,A)
#transitions
A = cbind(as.vector(apply(100*model$pi2_1[(n1*(1-burn)+1):n1,,], c(2,3), mean)),
as.vector(apply(100*model$pi2_1[(n1*(1-burn)+1):n1,,], c(2,3), sd)),
as.vector(apply(100*model$pi2_1[(n1*(1-burn)+1):n1,,], c(2,3), quantile,probs=0.025)),
as.vector(apply(100*model$pi2_1[(n1*(1-burn)+1):n1,,], c(2,3), quantile,probs=0.5)),
as.vector(apply(100*model$pi2_1[(n1*(1-burn)+1):n1,,], c(2,3), quantile,probs=0.975)),
NA)
colnames(A) = colnames(df)
rn = c()
for (k2 in 1:K2) {
for (k1 in 1:K1) {
rn = c(rn,paste("2->1,",k2,"->",k1,sep=''))
}
}
rownames(A) = rn
df = rbind(df,A)
#joint
A = cbind(as.vector(t(apply(100*model$pi12[(n1*(1-burn)+1):n1,,], c(2,3), mean))),
as.vector(t(apply(100*model$pi12[(n1*(1-burn)+1):n1,,], c(2,3), sd))),
as.vector(t(apply(100*model$pi12[(n1*(1-burn)+1):n1,,], c(2,3), quantile,probs=0.025))),
as.vector(t(apply(100*model$pi12[(n1*(1-burn)+1):n1,,], c(2,3), quantile,probs=0.5))),
as.vector(t(apply(100*model$pi12[(n1*(1-burn)+1):n1,,], c(2,3), quantile,probs=0.975))),
NA)
colnames(A) = colnames(df)
rn = c()
for (k1 in 1:K1) {
for (k2 in 1:K2) {
rn = c(rn,paste(k1,k2,sep=','))
}
}
rownames(A) = rn
df = rbind(df,A)
#relabel dataframe
colnames(df)[1:2] = c("Estimate","Standard Deviation")
colnames(df)[6] = "Inclusion Prob."
#Get likelihood and BIC at posterior median
id1 = X1[,1]
id2 = X2[,1]
X1[,1]=1
X2[,1]=1
z1 = (beta1.e != 0)*1
z2 = (beta2.e != 0)*1
ll = log_lik_dual(X1,X2,y1,y2,pi1.e,pi1_2.e,beta1.e,beta2.e,sigma1.e,sigma2.e,id1,id2)
BIC = BIC_calc_dual(X1,X2,y1,y2,pi1.e,pi1_2.e,beta1.e,beta2.e,sigma1.e,sigma2.e,id1,id2,z1,z2)
return(list(estimates = df,
log.likelihood = ll,
BIC = BIC))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.