R/BMSC-main.R
In michelescandola/BMSC: BMSC: Bayesian Multilevel Single Case models using 'Stan'
Defines functions
BMSC
.building.model
.building.data.list
summary.BMSC
print.summary.BMSC
pp_check.BMSC
plot.BMSC
Documented in
BMSC
pp_check.BMSC
#' BMSC: Bayesian Multilevel Single Case models using 'Stan'
#'
#' The \strong{BMSC} package provides an interface to fit Bayesian Multilevel Single Case models.
#' These models compare the performance of a single patient against a control group, combining
#' the flexibility of multilevel models and the potentiality of Bayesian Statistics.
#'
#' The package is now limited to gaussian data only, but we will further expand it to cover
#' binomial and ordinal (Likert scales) data.
#'
#' By means of BMSC the effects of the control group and the effects of the deviance between the
#' patient and the group will be esimated.
#'
#' The model to estimate the controls parameters is:
#'
#' \deqn{y~\mathcal{N}(\beta X + b Z, \sigma^2)}
#'
#' where $y$ is the controls' dependent variable, $X$ the contrast matrix for Population-level (or Fixed)
#' Effects, and $\beta$ are the unknown coefficients to be estimate. $Z$ is the contrast matrix for the
#' Varying (or Random, or Group-level) effects, and $b$ are the unknown estimates for the varying effects.
#' $\sigma^2$ is the variance.
#'
#' In order to estimate the coefficients of the patient, the formula is the following:
#'
#' \deqn{y_{pt}~\mathcal{N}(\phi X_{pt}, \sigma_{pt}^2)}
#'
#' where $\phi = \beta + \delta$.
#'
#' @section Details:
#' The main function of \strong{BMSC} is \code{\link{BMSC}}, which uses formula syntax to
#' specify your model.
#'
#' @docType package
#' @name BMSC
NULL
#' Fit Bayesian Multilevel Single Case models
#'
#' \code{BMSC} fits the Bayesian Multilevel Single Case models.
#'
#' @usage BMSC(formula, data.ctrl, data.pt, cores = 1, chains = 4, warmup = 2000, iter = 4000, seed = NA, control = list(adapt_delta=0.8))
#'
#' @param formula An object of class \code{formula}: a symbolic description of the model to be fitted.
#' @param data.ctrl An object of class \code{data.frame} (or one that can be coerced to that class)
#' containing data of all variables used in the model for the control group.
#' @param data.pt An object of class \code{data.frame} (or one that can be coerced to that class)
#' containing data of all variables used in the model for the patient.
#' @param chains Number of Markov chains (defaults to 4).
#' @param iter Number of total iterations per chain (including warmup; defaults to 4000).
#' @param warmup A positive integer specifying number of warmup (aka burnin) iterations.
#' This also specifies the number of iterations used for stepsize adaptation,
#' so warmup samples should not be used for inference. The number of warmup should
#' not be larger than iter and the default is 2000.
#' @param seed The seed for random number generation to make results reproducible.
#' If NA (the default), Stan will set the seed randomly.
#' @param typeprior Set the desired prior distribution for the fixed effects.
#' \describe{
#' \item{normal} a normal distribution with \mu = 0 and \sigma = 10
#' \item{cauchy} a cauchy distribution with \mu = 0 and scale \sqrt(2)/2
#' \item{student} a Student's T distribution, with \mu = 0, \nu = 3 and \sigma = 10
#' }
#' The normal distribution is the default.
#' @param control A named list of parameters to control the sampler's behavior.
#' It defaults to NULL so all the default values are
#' used. For a comprehensive overview see \strong{stan}.
#'
#' @export
BMSC = function(formula,data.ctrl,data.pt,
cores = 1, chains = 4, warmup = 2000,
iter = 4000, seed = NA, typeprior = "normal",
control = list(adapt_delta=0.8)){
mypaste = function(x){
out = NULL
if(length(x)>1){
for(xx in 1:(length(x)-1)) out = paste(out,x[xx],"+")
}
out = paste(out,x[length(x)])
return(out)
}
if(typeprior!="normal"&&typeprior!="cauchy"&&typeprior!="student")
stop("Not a valid typeprior")
# extract formula's terms
form.terms = attributes(terms(formula))$term.labels
# build contrasts matrices of fixed effects
fix.terms = form.terms[!(grepl("\\|",form.terms))]
fix.formula = paste0(" ~",mypaste(fix.terms))
matrix.fix.ctrl = model.matrix(as.formula(fix.formula),data.ctrl)
matrix.fix.pt = model.matrix(as.formula(fix.formula),data.pt)
# build contrasts matrices for random effects
ran.terms = form.terms[(grepl("\\|",form.terms))]
ran.matrices = list()
grouping = list()
for(ran in ran.terms){
tmp = unlist(strsplit(ran,"\\|"))
grouping[[ran]] = trimws(tmp[2])
ran.matrices[[ran]] = model.matrix(as.formula(paste0(" ~",mypaste(tmp[1]))),data.ctrl)
}
stancode = .building.model(ran.matrices,typeprior)
datalist = .building.data.list(ran.matrices,grouping,matrix.fix.ctrl,
matrix.fix.pt,data.ctrl,data.pt,formula)
mdl = stan(model_code = stancode, data = datalist, iter = iter,
chains = chains,cores = cores, warmup = warmup,
seed = seed, control = control)
out = list(formula,mdl,data.pt,data.ctrl,datalist,stancode,typeprior)
class(out) = append(class(out),"BMSC")
return(out)
}
.building.model = function(ran.matrices=NULL,typeprior){
code.data = "
data {
int<lower=1> Obs_Controls; // the total number of observations in the control group
int<lower=1> Obs_Patients; // the total numebr of observation in the patient
int<lower=1> Nparameters; // the total number of parameters for the independent variables
real y_Ctrl[Obs_Controls]; // the dependent variable for the control group
real y_Pts[Obs_Patients]; // the patient d.v.
matrix[Obs_Controls,Nparameters] XF_Ctrl; // the control matrix
matrix[Obs_Patients,Nparameters] XF_Pts; // the patient matrix"
code.parameter = "parameters {
vector[Nparameters] b_Ctrl; //the regression parameters for controls
vector[Nparameters] b_Delta; //the regression parameters for the controls - patients difference
real<lower=0> sigmaC; //the standard deviation for controls
real<lower=0> sigmaP; //the standard deviation for patient"
code.transformed.parameter = "transformed parameters{
real mu_Pts[Obs_Patients];
real mu_Ctrl[Obs_Controls];"
last.string.code.transformed.parameter = "
for(i in 1:Obs_Patients){
mu_Pts[i] = dot_product(b_Ctrl+b_Delta,XF_Pts[i,]);
}
for(i in 1:Obs_Controls){
mu_Ctrl[i] = dot_product(b_Ctrl,XF_Ctrl[i,])"
if(typeprior=="normal"){
code.model = "
target += cauchy_lpdf(sigmaC|0,1000);
target += cauchy_lpdf(sigmaP|0,1000);
target += normal_lpdf(b_Ctrl | 0, 10);
target += normal_lpdf(b_Delta | 0, 10);
target += cauchy_lpdf(y_Pts|mu_Pts,sigmaP);
target += cauchy_lpdf(y_Ctrl|mu_Ctrl,sigmaC);
}"
}else if(typeprior=="cauchy"){
code.model = "
target += cauchy_lpdf(sigmaC|0,1000);
target += cauchy_lpdf(sigmaP|0,1000);
target += cauchy_lpdf(b_Ctrl | 0, sqrt(2)/2);
target += cauchy_lpdf(b_Delta | 0, sqrt(2)/2);
target += cauchy_lpdf(y_Pts|mu_Pts,sigmaP);
target += cauchy_lpdf(y_Ctrl|mu_Ctrl,sigmaC);
}"
}else if(typeprior=="student"){
code.model = "
target += cauchy_lpdf(sigmaC|0,1000);
target += cauchy_lpdf(sigmaP|0,1000);
target += student_t_lpdf(b_Ctrl | 3, 0, 10);
target += student_t_lpdf(b_Delta | 3, 0, 10);
target += cauchy_lpdf(y_Pts|mu_Pts,sigmaP);
target += cauchy_lpdf(y_Ctrl|mu_Ctrl,sigmaC);
}"
}
code.generated.quantities = "generated quantities {
real y_pt_rep[Obs_Patients];
real y_ct_rep[Obs_Controls];
for(i in 1:Obs_Patients){
y_pt_rep[i] = cauchy_rng(mu_Pts[i], sigmaP);
}
for(i in 1:Obs_Controls){
y_ct_rep[i] = cauchy_rng(mu_Ctrl[i], sigmaC);
}
}"
if(!is.null(ran.matrices)){
ir = 1
for(ran in ran.matrices){
if(ncol(ran)>1){
code.data = paste(code.data,
paste0(" int<lower=1> Nrands",ir,"; //number of random coefficients for grouping factor ",ir),
paste0(" matrix[Obs_Controls,Nrands",ir,"] XR_Ctrl",ir,"; //the control random matrix for grouping factor ",ir),
paste0(" int grouping",ir,"[Obs_Controls]; //the index vector for the grouping factor ",ir),
paste0(" int<lower=1> Ngrouping",ir,"; // the total number of levels for grouping",ir),
sep ="\n"
)
code.parameter = paste(code.parameter,
paste0(" vector<lower=0>[Nrands",ir,"] sigma_u",ir,"; // random effects sd for grouping factor ",ir),
paste0(" cholesky_factor_corr[Nrands",ir,"] L_Omega",ir,";"),
paste0(" matrix[Nrands",ir,",Ngrouping",ir,"] z_u",ir,";"),
sep="\n"
)
code.transformed.parameter = paste(code.transformed.parameter,
paste0(" matrix[Nrands",ir,",Ngrouping",ir,"] u",ir,";"),
paste0(" u",ir," = (diag_pre_multiply(sigma_u",ir,", L_Omega",ir,") * z_u",ir,"); //random effects for grouping factor",ir),
sep="\n")
last.string.code.transformed.parameter = paste(last.string.code.transformed.parameter,
paste0("+ dot_product(u",ir,"[,grouping",ir,"[i]],XR_Ctrl",ir,"[i,])"))
code.model = paste(paste0("target += lkj_corr_cholesky_lpdf(L_Omega",ir," | 1);"),
paste0("target += normal_lpdf(to_vector(z_u",ir,") | 0, 1);"),
code.model,sep="\n")
}else if(dim(ran)[2]==1){
code.data = paste(code.data,
paste0(" int grouping",ir,"[Obs_Controls]; //the index vector for the grouping factor ",ir),
paste0(" int<lower=1> Ngrouping",ir,"; // the total number of levels for grouping",ir),
sep ="\n"
)
code.parameter = paste(code.parameter,
paste0(" real<lower=0> sigma_u",ir,"; // random effects sd for grouping factor ",ir),
paste0(" real[Ngrouping",ir,"] u",ir,";"),
sep="\n"
)
last.string.code.transformed.parameter = paste(last.string.code.transformed.parameter,
paste0("+ u",ir,"[grouping",ir,"[i]]"))
code.model = paste(paste0("target += normal_lpdf(u",ir," | 0, 10);"),
code.model,sep="\n")
}
ir = ir +1
}
}
code.transformed.parameter = paste(code.transformed.parameter,
paste0(last.string.code.transformed.parameter,";"),
" }",sep="\n")
code.model = paste(" model{
//priors",code.model,sep="\n")
out = paste(code.data," }",
code.parameter," }",
code.transformed.parameter," }",
code.model,
code.generated.quantities,sep="\n")
return(out)
}
.building.data.list = function(ran.matrices=NULL,grouping,
matrix.fix.ctrl,matrix.fix.pt,
data.ctrl,data.pt,formula){
data.list = list(
Nparameters = ncol(matrix.fix.ctrl),
y_Ctrl=data.ctrl[,as.character(formula[2])],
y_Pts =data.pt[,as.character(formula[2])],
XF_Ctrl=matrix.fix.ctrl,
XF_Pts =matrix.fix.pt,
Obs_Controls = nrow(matrix.fix.ctrl),
Obs_Patients = nrow(matrix.fix.pt)
)
if(!is.null(ran.matrices)){
ir = 1
for(ran in ran.matrices){
nn = length(data.list)
if(ncol(ran)>1){
data.list[[paste0("Nrands",ir)]] = ncol(ran)
data.list[[paste0("XR_Ctrl",ir)]] = ran
data.list[[paste0("grouping",ir)]] = as.numeric(data.ctrl[,grouping[[ir]]])
data.list[[paste0("Ngrouping",ir)]] = length(unique(data.ctrl[,grouping[[ir]]]))
}else if(dim(ran)[2]==1){
data.list[[paste0("grouping",ir)]] = as.numeric(data.ctrl[,grouping[[ir]]])
data.list[[paste0("Ngrouping",ir)]] = length(unique(data.ctrl[,grouping[[ir]]]))
}
ir = ir + 1
}
}
return(data.list)
}
#' @export
summary.BMSC = function(x){
if(class(x)[2]!="BMSC") stop("Not a valid BMSC object.")
se <- function(x){sd(x)/sqrt(length(x))}
if(x[[7]]=="normal"){
d0 <- dnorm(0,0,10)
}else if(x[[7]]=="cauchy"){
d0 <- dcauchy(0,0,sqrt(2)/2)
}else if(x[[7]]=="student"){
d0 <- LaplacesDemon::dst(0,10,3)
}
delta = extract(x[[2]],pars="b_Delta")
delta_logspl = apply(delta$b_Delta,2,logspline)
BF10_delta = lapply(delta_logspl,FUN=function(x){d0/dlogspline(0,x)})
beta = extract(x[[2]],pars="b_Ctrl")
beta_logspl = apply(beta$b_Ctrl,2,logspline)
BF10_beta = lapply(beta_logspl,FUN=function(x){d0/dlogspline(0,x)})
pts = beta$b_Ctrl+delta$b_Delta
pts_logspl = apply(pts,2,logspline)
BF10_pts = lapply(pts_logspl,FUN=function(x){d0/dlogspline(0,x)})
sum01 = as.data.frame(summary(x[[2]],pars="b_Ctrl")[[1]])
sum02 = as.data.frame(summary(x[[2]],pars="sigmaC")[[1]])
sum03 = as.data.frame(summary(x[[2]],pars="b_Delta")[[1]])
sum04 = as.data.frame(summary(x[[2]],pars="sigmaP")[[1]])
sum05 = as.data.frame(cbind(apply(pts,2,mean),apply(pts,2,se),apply(pts,2,sd),
apply(pts,2,quantile,probs=2.5/100),
apply(pts,2,quantile,probs=25/100),
apply(pts,2,quantile,probs=50/100),
apply(pts,2,quantile,probs=75/100),
apply(pts,2,quantile,probs=97.5/100)))
colnames(sum05) = c("mean","se_mean","sd","2.5%","25%","50%","75%","97.5%")
rownames(sum01) <- colnames(x[[5]]$XF_Ctrl)
rownames(sum03) <- colnames(x[[5]]$XF_Pts)
rownames(sum05) <- colnames(x[[5]]$XF_Pts)
sum01$BF10 <- BF10_beta
sum03$BF10 <- BF10_delta
sum05$BF10 <- BF10_pts
out = list(sum01,sum02,sum03,sum04,x,sum05,x[[7]])
class(out) = append(class(out),"summary.BMSC")
return(out)
}
#' @export
print.summary.BMSC = function(x, ...){
cat("\nBayesian Multilevel Single Case model\n\n")
print(x[[5]][[1]], ...)
cat("\n")
print(paste("Prior:",x[[7]]))
cat("\n\n Fixed Effects for the Control Group\n\n")
print(x[[1]], ...)
cat("\n")
print(x[[2]], ...)
cat("\n\n Fixed Effects for the Patient\n\n")
print(x[[6]], ...)
cat("\n\n Fixed Effects for the difference between the Patient and the Control Group\n\n")
print(x[[3]], ...)
cat("\n")
print(x[[4]], ...)
}
#' Posterior predictive checks for BMSC objects
#'
#' \code{pp_check()} plots the posterior predictive check for BMSC objects.
#'
#' @param object a BMSC object
#' @param type
#' \describe{
#' \item{dens} density overlay plot
#' \item{hist} histogram plot
#' \item{mean} the distribution of the mean statistic, over the simulated datasets, compared to the mean of the real data
#' }
#' @return a ggplot2 object
#' @export
pp_check.BMSC = function(object,type="dens"){
if(class(x)[2]!="BMSC") stop("Not a valid BMSC object.")
y_ct_rep = extract(object[[2]], pars = "y_ct_rep")
y_pt_rep = extract(object[[2]], pars = "y_pt_rep")
y_ct = object[[4]][,as.character(object[[1]][2])]
y_pt = object[[3]][,as.character(object[[1]][2])]
ans=list()
if ((requireNamespace("bayesplot", quietly = TRUE))&&(requireNamespace("gridExtra", quietly = TRUE))) {
if(type=="hist"){
ct = ppc_hist(y_ct, y_ct_rep[[1]][1:8, ])+ggtitle("Control Group")+xlim(c(-sd(y_ct)*4,sd(y_ct)*4))
pt = ppc_hist(y_pt, y_pt_rep[[1]][1:8, ])+ggtitle("Patient")+xlim(c(-sd(y_pt)*4,sd(y_pt)*4))
ans = gridExtra::grid.arrange(ct,pt)
}else if(type=="mean"){
ct = ppc_stat(y_ct, y_ct_rep[[1]])+ggtitle("Control Group")+xlim(c(-sd(y_ct)*4,sd(y_ct)*4))
pt = ppc_stat(y_pt, y_pt_rep[[1]])+ggtitle("Patient")+xlim(c(-sd(y_pt)*4,sd(y_pt)*4))
ans = gridExtra::grid.arrange(ct,pt)
}else{
ct = ppc_dens_overlay(y_ct, y_ct_rep[[1]][1:200, ])+ggtitle("Control Group")+xlim(c(-sd(y_ct)*4,sd(y_ct)*4))
pt = ppc_dens_overlay(y_pt, y_pt_rep[[1]][1:200, ])+ggtitle("Patient")+xlim(c(-sd(y_pt)*4,sd(y_pt)*4))
ans = gridExtra::grid.arrange(ct,pt)
}
}
return(ans)
}
#' Plot estimates from a \code{BMSC} object.
#'
#' @usage plot(object, who = "both", type = "interval")
#'
#' @param object An object of class \code{BMSC}.
#' @param who parameter to choose the estimates to plot
#' \describe{
#' \item{both} plot in the same graph both controls and the patient
#' \item{control} only the controls
#' \item{patient} only the patient (\beta + \delta)
#' \item{delta} only the difference between the patient and controls
#' }
#' @param type a parameter to select the typology of graph
#' \describe{
#' \item{interval} the estimates will be represented by means of pointrange, with median and the boundaries of the credible interval
#' \item{area} a density plot
#' \item{hist} a density histogram
#' }
#' @param CI the dimension of the Credible Interval
#'
#' @return a plot, a ggplot2 object, or a bayesplot object
#'
#' @export
plot.BMSC = function(mdl, who = "both", type = "interval", CI = 0.95){
if(class(mdl)[2]!="BMSC") stop("Not a valid BMSC object.")
limits = c((1-CI)/2,1-((1-CI)/2))
if(type=="interval"){
if(who == "both"){
beta <- extract(mdl[[2]],pars="b_Ctrl")
controls <- apply(beta[[1]],2,quantile,probs=c(limits[1],0.5,limits[2]))
delta <- extract(mdl[[2]],pars="b_Delta")
patient <- apply(delta[[1]]+beta[[1]],2,quantile,probs=c(limits[1],0.5,limits[2]))
namC <- colnames(mdl[[5]]$XF_Ctrl)
namP <- colnames(mdl[[5]]$XF_Pts)
dat <- data.frame(
low = c(patient[1,],controls[1,]),
med = c(patient[2,],controls[2,]),
high = c(patient[3,],controls[3,]),
group = c(rep("Patient",ncol(patient)),rep("Controls",ncol(controls))),
i = ordered(c(namP,namC),
level=namC[length(namC):1])
)
ans <- ggplot(dat,aes(x=i,ymin=low,ymax=high,y=med,colour=group))+
geom_pointrange()+
coord_flip()+xlab("coefficients")+ylab("")
}else if(who=="control"){
beta <- extract(mdl[[2]],pars="b_Ctrl")
controls <- apply(beta[[1]],2,quantile,probs=c(limits[1],0.5,limits[2]))
namC <- colnames(mdl[[5]]$XF_Ctrl)
dat <- data.frame(
low = controls[1,],
med = controls[2,],
high = controls[3,],
i = ordered(namC,level=namC[length(namC):1])
)
ans <- ggplot(dat,aes(x=i,ymin=low,ymax=high,y=med))+
geom_pointrange()+
coord_flip()+xlab("coefficients")+ylab("")
}else if(who=="patient"){
beta <- extract(mdl[[2]],pars="b_Ctrl")
delta <- extract(mdl[[2]],pars="b_Delta")
patient <- apply(delta[[1]]+beta[[1]],2,quantile,probs=c(limits[1],0.5,limits[2]))
namP <- colnames(mdl[[5]]$XF_Pts)
dat <- data.frame(
low = patient[1,],
med = patient[2,],
high = patient[3,],
i = ordered(namP,level=namP[length(namP):1])
)
ans <- ggplot(dat,aes(x=i,ymin=low,ymax=high,y=med))+
geom_pointrange()+
coord_flip()+xlab("coefficients")+ylab("")
}else if(who=="delta"){
delta <- extract(mdl[[2]],pars="b_Delta")
patient <- apply(delta[[1]],2,quantile,probs=c(limits[1],0.5,limits[2]))
namP <- colnames(mdl[[5]]$XF_Pts)
dat <- data.frame(
low = patient[1,],
med = patient[2,],
high = patient[3,],
i = ordered(namP,level=namP[length(namP):1])
)
ans <- ggplot(dat,aes(x=i,ymin=low,ymax=high,y=med))+
geom_pointrange()+
coord_flip()+xlab("coefficients")+ylab("")
}
} else if (type == "area") {
if(requireNamespace("reshape2", quietly = TRUE)){
if(who == "both"){
beta <- extract(mdl[[2]],pars="b_Ctrl")
controls <- reshape2::melt(beta[[1]])
namC <- colnames(mdl[[5]]$XF_Ctrl)
controls$Var2 <- factor(controls$Var2)
levels(controls$Var2) <- namC
controls$Var2 <- ordered(controls$Var2,level=namC[length(namC):1])
controls$group = "Controls"
delta <- extract(mdl[[2]],pars="b_Delta")
tmp <- delta[[1]]+beta[[1]]
patient <- reshape2::melt(tmp)
namP <- colnames(mdl[[5]]$XF_Pts)
patient$Var2 <- factor(patient$Var2)
levels(patient$Var2) <- namP
patient$Var2 <- ordered(patient$Var2,level=namP[length(namP):1])
patient$group = "Patient"
dat <- rbind(
patient,controls
)
ans <- ggplot(dat,aes(x=value,colour=group))+
geom_density()+
facet_grid(Var2~.)+
xlab("coefficients")+ylab("")
}else if(who=="control"){
beta <- extract(mdl[[2]],pars="b_Ctrl")
controls <- reshape2::melt(beta[[1]])
namC <- colnames(mdl[[5]]$XF_Ctrl)
controls$Var2 <- factor(controls$Var2)
levels(controls$Var2) <- namC
controls$Var2 <- ordered(controls$Var2,level=namC[length(namC):1])
dat <- controls
ans <- ggplot(dat,aes(x=value))+
geom_density()+
facet_grid(Var2~.)+
xlab("coefficients")+ylab("")
}else if(who=="patient"){
beta <- extract(mdl[[2]],pars="b_Ctrl")
delta <- extract(mdl[[2]],pars="b_Delta")
tmp <- delta[[1]]+beta[[1]]
patient <- reshape2::melt(tmp)
namP <- colnames(mdl[[5]]$XF_Pts)
patient$Var2 <- factor(patient$Var2)
levels(patient$Var2) <- namP
patient$Var2 <- ordered(patient$Var2,level=namP[length(namP):1])
dat <- patient
ans <- ggplot(dat,aes(x=value))+
geom_density()+
facet_grid(Var2~.)+
xlab("coefficients")+ylab("")
}else if(who=="delta"){
delta <- extract(mdl[[2]],pars="b_Delta")
patient <- reshape2::melt(delta[[1]])
namP <- colnames(mdl[[5]]$XF_Pts)
patient$Var2 <- factor(patient$Var2)
levels(patient$Var2) <- namP
patient$Var2 <- ordered(patient$Var2,level=namP[length(namP):1])
dat <- patient
ans <- ggplot(dat,aes(x=value))+
geom_density()+
facet_grid(Var2~.)+
xlab("coefficients")+ylab("")
}
}
} else if (type == "hist") {
if(requireNamespace("reshape2", quietly = TRUE)){
if(who == "both"){
beta <- extract(mdl[[2]],pars="b_Ctrl")
controls <- reshape2::melt(beta[[1]])
namC <- colnames(mdl[[5]]$XF_Ctrl)
controls$Var2 <- factor(controls$Var2)
levels(controls$Var2) <- namC
controls$Var2 <- ordered(controls$Var2,level=namC[length(namC):1])
controls$group = "Controls"
delta <- extract(mdl[[2]],pars="b_Delta")
tmp <- delta[[1]]+beta[[1]]
patient <- reshape2::melt(tmp)
namP <- colnames(mdl[[5]]$XF_Pts)
patient$Var2 <- factor(patient$Var2)
levels(patient$Var2) <- namP
patient$Var2 <- ordered(patient$Var2,level=namP[length(namP):1])
patient$group = "Patient"
dat <- rbind(
patient,controls
)
ans <- ggplot(dat,aes(x=value,fill=group))+
geom_histogram()+
facet_grid(Var2~.)+
xlab("coefficients")+ylab("")
}else if(who=="control"){
beta <- extract(mdl[[2]],pars="b_Ctrl")
controls <- reshape2::melt(beta[[1]])
namC <- colnames(mdl[[5]]$XF_Ctrl)
controls$Var2 <- factor(controls$Var2)
levels(controls$Var2) <- namC
controls$Var2 <- ordered(controls$Var2,level=namC[length(namC):1])
dat <- controls
ans <- ggplot(dat,aes(x=value))+
geom_histogram()+
facet_grid(Var2~.)+
xlab("coefficients")+ylab("")
}else if(who=="patient"){
beta <- extract(mdl[[2]],pars="b_Ctrl")
delta <- extract(mdl[[2]],pars="b_Delta")
tmp <- delta[[1]]+beta[[1]]
patient <- reshape2::melt(tmp)
namP <- colnames(mdl[[5]]$XF_Pts)
patient$Var2 <- factor(patient$Var2)
levels(patient$Var2) <- namP
patient$Var2 <- ordered(patient$Var2,level=namP[length(namP):1])
dat <- patient
ans <- ggplot(dat,aes(x=value))+
geom_histogram()+
facet_grid(Var2~.)+
xlab("coefficients")+ylab("")
}else if(who=="delta"){
delta <- extract(mdl[[2]],pars="b_Delta")
patient <- reshape2::melt(delta[[1]])
namP <- colnames(mdl[[5]]$XF_Pts)
patient$Var2 <- factor(patient$Var2)
levels(patient$Var2) <- namP
patient$Var2 <- ordered(patient$Var2,level=namP[length(namP):1])
dat <- patient
ans <- ggplot(dat,aes(x=value))+
geom_histogram()+
facet_grid(Var2~.)+
xlab("coefficients")+ylab("")
}
}
}
return(ans)
}
#' Computes log marginal likelihood via bridge sampling.
#'
#' @param object a BMSC object
#' @return an "psis_loo" "loo" object
#' @export
bridge_sampler.BMSC = function(object, ...){
if (requireNamespace("bridgesampling", quietly = TRUE)) {
ans = bridgesampling::bridge_sampler(object[[2]],...)
}
return(ans)
}
michelescandola/BMSC documentation built on Dec. 7, 2019, 7:46 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions BMSC .building.model .building.data.list summary.BMSC print.summary.BMSC pp_check.BMSC plot.BMSC
Documented in BMSC pp_check.BMSC
#' BMSC: Bayesian Multilevel Single Case models using 'Stan'
#'
#' The \strong{BMSC} package provides an interface to fit Bayesian Multilevel Single Case models.
#' These models compare the performance of a single patient against a control group, combining
#' the flexibility of multilevel models and the potentiality of Bayesian Statistics.
#'
#' The package is now limited to gaussian data only, but we will further expand it to cover
#' binomial and ordinal (Likert scales) data.
#'
#' By means of BMSC the effects of the control group and the effects of the deviance between the
#' patient and the group will be esimated.
#'
#' The model to estimate the controls parameters is:
#'
#' \deqn{y~\mathcal{N}(\beta X + b Z, \sigma^2)}
#'
#' where $y$ is the controls' dependent variable, $X$ the contrast matrix for Population-level (or Fixed)
#' Effects, and $\beta$ are the unknown coefficients to be estimate. $Z$ is the contrast matrix for the
#' Varying (or Random, or Group-level) effects, and $b$ are the unknown estimates for the varying effects.
#' $\sigma^2$ is the variance.
#'
#' In order to estimate the coefficients of the patient, the formula is the following:
#'
#' \deqn{y_{pt}~\mathcal{N}(\phi X_{pt}, \sigma_{pt}^2)}
#'
#' where $\phi = \beta + \delta$.
#'
#' @section Details:
#' The main function of \strong{BMSC} is \code{\link{BMSC}}, which uses formula syntax to
#' specify your model.
#'
#' @docType package
#' @name BMSC
NULL
#' Fit Bayesian Multilevel Single Case models
#'
#' \code{BMSC} fits the Bayesian Multilevel Single Case models.
#'
#' @usage BMSC(formula, data.ctrl, data.pt, cores = 1, chains = 4, warmup = 2000, iter = 4000, seed = NA, control = list(adapt_delta=0.8))
#'
#' @param formula An object of class \code{formula}: a symbolic description of the model to be fitted.
#' @param data.ctrl An object of class \code{data.frame} (or one that can be coerced to that class)
#' containing data of all variables used in the model for the control group.
#' @param data.pt An object of class \code{data.frame} (or one that can be coerced to that class)
#' containing data of all variables used in the model for the patient.
#' @param chains Number of Markov chains (defaults to 4).
#' @param iter Number of total iterations per chain (including warmup; defaults to 4000).
#' @param warmup A positive integer specifying number of warmup (aka burnin) iterations.
#' This also specifies the number of iterations used for stepsize adaptation,
#' so warmup samples should not be used for inference. The number of warmup should
#' not be larger than iter and the default is 2000.
#' @param seed The seed for random number generation to make results reproducible.
#' If NA (the default), Stan will set the seed randomly.
#' @param typeprior Set the desired prior distribution for the fixed effects.
#' \describe{
#' \item{normal} a normal distribution with \mu = 0 and \sigma = 10
#' \item{cauchy} a cauchy distribution with \mu = 0 and scale \sqrt(2)/2
#' \item{student} a Student's T distribution, with \mu = 0, \nu = 3 and \sigma = 10
#' }
#' The normal distribution is the default.
#' @param control A named list of parameters to control the sampler's behavior.
#' It defaults to NULL so all the default values are
#' used. For a comprehensive overview see \strong{stan}.
#'
#' @export
BMSC = function(formula,data.ctrl,data.pt,
cores = 1, chains = 4, warmup = 2000,
iter = 4000, seed = NA, typeprior = "normal",
control = list(adapt_delta=0.8)){
mypaste = function(x){
out = NULL
if(length(x)>1){
for(xx in 1:(length(x)-1)) out = paste(out,x[xx],"+")
}
out = paste(out,x[length(x)])
return(out)
}
if(typeprior!="normal"&&typeprior!="cauchy"&&typeprior!="student")
stop("Not a valid typeprior")
# extract formula's terms
form.terms = attributes(terms(formula))$term.labels
# build contrasts matrices of fixed effects
fix.terms = form.terms[!(grepl("\\|",form.terms))]
fix.formula = paste0(" ~",mypaste(fix.terms))
matrix.fix.ctrl = model.matrix(as.formula(fix.formula),data.ctrl)
matrix.fix.pt = model.matrix(as.formula(fix.formula),data.pt)
# build contrasts matrices for random effects
ran.terms = form.terms[(grepl("\\|",form.terms))]
ran.matrices = list()
grouping = list()
for(ran in ran.terms){
tmp = unlist(strsplit(ran,"\\|"))
grouping[[ran]] = trimws(tmp[2])
ran.matrices[[ran]] = model.matrix(as.formula(paste0(" ~",mypaste(tmp[1]))),data.ctrl)
}
stancode = .building.model(ran.matrices,typeprior)
datalist = .building.data.list(ran.matrices,grouping,matrix.fix.ctrl,
matrix.fix.pt,data.ctrl,data.pt,formula)
mdl = stan(model_code = stancode, data = datalist, iter = iter,
chains = chains,cores = cores, warmup = warmup,
seed = seed, control = control)
out = list(formula,mdl,data.pt,data.ctrl,datalist,stancode,typeprior)
class(out) = append(class(out),"BMSC")
return(out)
}
.building.model = function(ran.matrices=NULL,typeprior){
code.data = "
data {
int<lower=1> Obs_Controls; // the total number of observations in the control group
int<lower=1> Obs_Patients; // the total numebr of observation in the patient
int<lower=1> Nparameters; // the total number of parameters for the independent variables
real y_Ctrl[Obs_Controls]; // the dependent variable for the control group
real y_Pts[Obs_Patients]; // the patient d.v.
matrix[Obs_Controls,Nparameters] XF_Ctrl; // the control matrix
matrix[Obs_Patients,Nparameters] XF_Pts; // the patient matrix"
code.parameter = "parameters {
vector[Nparameters] b_Ctrl; //the regression parameters for controls
vector[Nparameters] b_Delta; //the regression parameters for the controls - patients difference
real<lower=0> sigmaC; //the standard deviation for controls
real<lower=0> sigmaP; //the standard deviation for patient"
code.transformed.parameter = "transformed parameters{
real mu_Pts[Obs_Patients];
real mu_Ctrl[Obs_Controls];"
last.string.code.transformed.parameter = "
for(i in 1:Obs_Patients){
mu_Pts[i] = dot_product(b_Ctrl+b_Delta,XF_Pts[i,]);
}
for(i in 1:Obs_Controls){
mu_Ctrl[i] = dot_product(b_Ctrl,XF_Ctrl[i,])"
if(typeprior=="normal"){
code.model = "
target += cauchy_lpdf(sigmaC|0,1000);
target += cauchy_lpdf(sigmaP|0,1000);
target += normal_lpdf(b_Ctrl | 0, 10);
target += normal_lpdf(b_Delta | 0, 10);
target += cauchy_lpdf(y_Pts|mu_Pts,sigmaP);
target += cauchy_lpdf(y_Ctrl|mu_Ctrl,sigmaC);
}"
}else if(typeprior=="cauchy"){
code.model = "
target += cauchy_lpdf(sigmaC|0,1000);
target += cauchy_lpdf(sigmaP|0,1000);
target += cauchy_lpdf(b_Ctrl | 0, sqrt(2)/2);
target += cauchy_lpdf(b_Delta | 0, sqrt(2)/2);
target += cauchy_lpdf(y_Pts|mu_Pts,sigmaP);
target += cauchy_lpdf(y_Ctrl|mu_Ctrl,sigmaC);
}"
}else if(typeprior=="student"){
code.model = "
target += cauchy_lpdf(sigmaC|0,1000);
target += cauchy_lpdf(sigmaP|0,1000);
target += student_t_lpdf(b_Ctrl | 3, 0, 10);
target += student_t_lpdf(b_Delta | 3, 0, 10);
target += cauchy_lpdf(y_Pts|mu_Pts,sigmaP);
target += cauchy_lpdf(y_Ctrl|mu_Ctrl,sigmaC);
}"
}
code.generated.quantities = "generated quantities {
real y_pt_rep[Obs_Patients];
real y_ct_rep[Obs_Controls];
for(i in 1:Obs_Patients){
y_pt_rep[i] = cauchy_rng(mu_Pts[i], sigmaP);
}
for(i in 1:Obs_Controls){
y_ct_rep[i] = cauchy_rng(mu_Ctrl[i], sigmaC);
}
}"
if(!is.null(ran.matrices)){
ir = 1
for(ran in ran.matrices){
if(ncol(ran)>1){
code.data = paste(code.data,
paste0(" int<lower=1> Nrands",ir,"; //number of random coefficients for grouping factor ",ir),
paste0(" matrix[Obs_Controls,Nrands",ir,"] XR_Ctrl",ir,"; //the control random matrix for grouping factor ",ir),
paste0(" int grouping",ir,"[Obs_Controls]; //the index vector for the grouping factor ",ir),
paste0(" int<lower=1> Ngrouping",ir,"; // the total number of levels for grouping",ir),
sep ="\n"
)
code.parameter = paste(code.parameter,
paste0(" vector<lower=0>[Nrands",ir,"] sigma_u",ir,"; // random effects sd for grouping factor ",ir),
paste0(" cholesky_factor_corr[Nrands",ir,"] L_Omega",ir,";"),
paste0(" matrix[Nrands",ir,",Ngrouping",ir,"] z_u",ir,";"),
sep="\n"
)
code.transformed.parameter = paste(code.transformed.parameter,
paste0(" matrix[Nrands",ir,",Ngrouping",ir,"] u",ir,";"),
paste0(" u",ir," = (diag_pre_multiply(sigma_u",ir,", L_Omega",ir,") * z_u",ir,"); //random effects for grouping factor",ir),
sep="\n")
last.string.code.transformed.parameter = paste(last.string.code.transformed.parameter,
paste0("+ dot_product(u",ir,"[,grouping",ir,"[i]],XR_Ctrl",ir,"[i,])"))
code.model = paste(paste0("target += lkj_corr_cholesky_lpdf(L_Omega",ir," | 1);"),
paste0("target += normal_lpdf(to_vector(z_u",ir,") | 0, 1);"),
code.model,sep="\n")
}else if(dim(ran)[2]==1){
code.data = paste(code.data,
paste0(" int grouping",ir,"[Obs_Controls]; //the index vector for the grouping factor ",ir),
paste0(" int<lower=1> Ngrouping",ir,"; // the total number of levels for grouping",ir),
sep ="\n"
)
code.parameter = paste(code.parameter,
paste0(" real<lower=0> sigma_u",ir,"; // random effects sd for grouping factor ",ir),
paste0(" real[Ngrouping",ir,"] u",ir,";"),
sep="\n"
)
last.string.code.transformed.parameter = paste(last.string.code.transformed.parameter,
paste0("+ u",ir,"[grouping",ir,"[i]]"))
code.model = paste(paste0("target += normal_lpdf(u",ir," | 0, 10);"),
code.model,sep="\n")
}
ir = ir +1
}
}
code.transformed.parameter = paste(code.transformed.parameter,
paste0(last.string.code.transformed.parameter,";"),
" }",sep="\n")
code.model = paste(" model{
//priors",code.model,sep="\n")
out = paste(code.data," }",
code.parameter," }",
code.transformed.parameter," }",
code.model,
code.generated.quantities,sep="\n")
return(out)
}
.building.data.list = function(ran.matrices=NULL,grouping,
matrix.fix.ctrl,matrix.fix.pt,
data.ctrl,data.pt,formula){
data.list = list(
Nparameters = ncol(matrix.fix.ctrl),
y_Ctrl=data.ctrl[,as.character(formula[2])],
y_Pts =data.pt[,as.character(formula[2])],
XF_Ctrl=matrix.fix.ctrl,
XF_Pts =matrix.fix.pt,
Obs_Controls = nrow(matrix.fix.ctrl),
Obs_Patients = nrow(matrix.fix.pt)
)
if(!is.null(ran.matrices)){
ir = 1
for(ran in ran.matrices){
nn = length(data.list)
if(ncol(ran)>1){
data.list[[paste0("Nrands",ir)]] = ncol(ran)
data.list[[paste0("XR_Ctrl",ir)]] = ran
data.list[[paste0("grouping",ir)]] = as.numeric(data.ctrl[,grouping[[ir]]])
data.list[[paste0("Ngrouping",ir)]] = length(unique(data.ctrl[,grouping[[ir]]]))
}else if(dim(ran)[2]==1){
data.list[[paste0("grouping",ir)]] = as.numeric(data.ctrl[,grouping[[ir]]])
data.list[[paste0("Ngrouping",ir)]] = length(unique(data.ctrl[,grouping[[ir]]]))
}
ir = ir + 1
}
}
return(data.list)
}
#' @export
summary.BMSC = function(x){
if(class(x)[2]!="BMSC") stop("Not a valid BMSC object.")
se <- function(x){sd(x)/sqrt(length(x))}
if(x[[7]]=="normal"){
d0 <- dnorm(0,0,10)
}else if(x[[7]]=="cauchy"){
d0 <- dcauchy(0,0,sqrt(2)/2)
}else if(x[[7]]=="student"){
d0 <- LaplacesDemon::dst(0,10,3)
}
delta = extract(x[[2]],pars="b_Delta")
delta_logspl = apply(delta$b_Delta,2,logspline)
BF10_delta = lapply(delta_logspl,FUN=function(x){d0/dlogspline(0,x)})
beta = extract(x[[2]],pars="b_Ctrl")
beta_logspl = apply(beta$b_Ctrl,2,logspline)
BF10_beta = lapply(beta_logspl,FUN=function(x){d0/dlogspline(0,x)})
pts = beta$b_Ctrl+delta$b_Delta
pts_logspl = apply(pts,2,logspline)
BF10_pts = lapply(pts_logspl,FUN=function(x){d0/dlogspline(0,x)})
sum01 = as.data.frame(summary(x[[2]],pars="b_Ctrl")[[1]])
sum02 = as.data.frame(summary(x[[2]],pars="sigmaC")[[1]])
sum03 = as.data.frame(summary(x[[2]],pars="b_Delta")[[1]])
sum04 = as.data.frame(summary(x[[2]],pars="sigmaP")[[1]])
sum05 = as.data.frame(cbind(apply(pts,2,mean),apply(pts,2,se),apply(pts,2,sd),
apply(pts,2,quantile,probs=2.5/100),
apply(pts,2,quantile,probs=25/100),
apply(pts,2,quantile,probs=50/100),
apply(pts,2,quantile,probs=75/100),
apply(pts,2,quantile,probs=97.5/100)))
colnames(sum05) = c("mean","se_mean","sd","2.5%","25%","50%","75%","97.5%")
rownames(sum01) <- colnames(x[[5]]$XF_Ctrl)
rownames(sum03) <- colnames(x[[5]]$XF_Pts)
rownames(sum05) <- colnames(x[[5]]$XF_Pts)
sum01$BF10 <- BF10_beta
sum03$BF10 <- BF10_delta
sum05$BF10 <- BF10_pts
out = list(sum01,sum02,sum03,sum04,x,sum05,x[[7]])
class(out) = append(class(out),"summary.BMSC")
return(out)
}
#' @export
print.summary.BMSC = function(x, ...){
cat("\nBayesian Multilevel Single Case model\n\n")
print(x[[5]][[1]], ...)
cat("\n")
print(paste("Prior:",x[[7]]))
cat("\n\n Fixed Effects for the Control Group\n\n")
print(x[[1]], ...)
cat("\n")
print(x[[2]], ...)
cat("\n\n Fixed Effects for the Patient\n\n")
print(x[[6]], ...)
cat("\n\n Fixed Effects for the difference between the Patient and the Control Group\n\n")
print(x[[3]], ...)
cat("\n")
print(x[[4]], ...)
}
#' Posterior predictive checks for BMSC objects
#'
#' \code{pp_check()} plots the posterior predictive check for BMSC objects.
#'
#' @param object a BMSC object
#' @param type
#' \describe{
#' \item{dens} density overlay plot
#' \item{hist} histogram plot
#' \item{mean} the distribution of the mean statistic, over the simulated datasets, compared to the mean of the real data
#' }
#' @return a ggplot2 object
#' @export
pp_check.BMSC = function(object,type="dens"){
if(class(x)[2]!="BMSC") stop("Not a valid BMSC object.")
y_ct_rep = extract(object[[2]], pars = "y_ct_rep")
y_pt_rep = extract(object[[2]], pars = "y_pt_rep")
y_ct = object[[4]][,as.character(object[[1]][2])]
y_pt = object[[3]][,as.character(object[[1]][2])]
ans=list()
if ((requireNamespace("bayesplot", quietly = TRUE))&&(requireNamespace("gridExtra", quietly = TRUE))) {
if(type=="hist"){
ct = ppc_hist(y_ct, y_ct_rep[[1]][1:8, ])+ggtitle("Control Group")+xlim(c(-sd(y_ct)*4,sd(y_ct)*4))
pt = ppc_hist(y_pt, y_pt_rep[[1]][1:8, ])+ggtitle("Patient")+xlim(c(-sd(y_pt)*4,sd(y_pt)*4))
ans = gridExtra::grid.arrange(ct,pt)
}else if(type=="mean"){
ct = ppc_stat(y_ct, y_ct_rep[[1]])+ggtitle("Control Group")+xlim(c(-sd(y_ct)*4,sd(y_ct)*4))
pt = ppc_stat(y_pt, y_pt_rep[[1]])+ggtitle("Patient")+xlim(c(-sd(y_pt)*4,sd(y_pt)*4))
ans = gridExtra::grid.arrange(ct,pt)
}else{
ct = ppc_dens_overlay(y_ct, y_ct_rep[[1]][1:200, ])+ggtitle("Control Group")+xlim(c(-sd(y_ct)*4,sd(y_ct)*4))
pt = ppc_dens_overlay(y_pt, y_pt_rep[[1]][1:200, ])+ggtitle("Patient")+xlim(c(-sd(y_pt)*4,sd(y_pt)*4))
ans = gridExtra::grid.arrange(ct,pt)
}
}
return(ans)
}
#' Plot estimates from a \code{BMSC} object.
#'
#' @usage plot(object, who = "both", type = "interval")
#'
#' @param object An object of class \code{BMSC}.
#' @param who parameter to choose the estimates to plot
#' \describe{
#' \item{both} plot in the same graph both controls and the patient
#' \item{control} only the controls
#' \item{patient} only the patient (\beta + \delta)
#' \item{delta} only the difference between the patient and controls
#' }
#' @param type a parameter to select the typology of graph
#' \describe{
#' \item{interval} the estimates will be represented by means of pointrange, with median and the boundaries of the credible interval
#' \item{area} a density plot
#' \item{hist} a density histogram
#' }
#' @param CI the dimension of the Credible Interval
#'
#' @return a plot, a ggplot2 object, or a bayesplot object
#'
#' @export
plot.BMSC = function(mdl, who = "both", type = "interval", CI = 0.95){
if(class(mdl)[2]!="BMSC") stop("Not a valid BMSC object.")
limits = c((1-CI)/2,1-((1-CI)/2))
if(type=="interval"){
if(who == "both"){
beta <- extract(mdl[[2]],pars="b_Ctrl")
controls <- apply(beta[[1]],2,quantile,probs=c(limits[1],0.5,limits[2]))
delta <- extract(mdl[[2]],pars="b_Delta")
patient <- apply(delta[[1]]+beta[[1]],2,quantile,probs=c(limits[1],0.5,limits[2]))
namC <- colnames(mdl[[5]]$XF_Ctrl)
namP <- colnames(mdl[[5]]$XF_Pts)
dat <- data.frame(
low = c(patient[1,],controls[1,]),
med = c(patient[2,],controls[2,]),
high = c(patient[3,],controls[3,]),
group = c(rep("Patient",ncol(patient)),rep("Controls",ncol(controls))),
i = ordered(c(namP,namC),
level=namC[length(namC):1])
)
ans <- ggplot(dat,aes(x=i,ymin=low,ymax=high,y=med,colour=group))+
geom_pointrange()+
coord_flip()+xlab("coefficients")+ylab("")
}else if(who=="control"){
beta <- extract(mdl[[2]],pars="b_Ctrl")
controls <- apply(beta[[1]],2,quantile,probs=c(limits[1],0.5,limits[2]))
namC <- colnames(mdl[[5]]$XF_Ctrl)
dat <- data.frame(
low = controls[1,],
med = controls[2,],
high = controls[3,],
i = ordered(namC,level=namC[length(namC):1])
)
ans <- ggplot(dat,aes(x=i,ymin=low,ymax=high,y=med))+
geom_pointrange()+
coord_flip()+xlab("coefficients")+ylab("")
}else if(who=="patient"){
beta <- extract(mdl[[2]],pars="b_Ctrl")
delta <- extract(mdl[[2]],pars="b_Delta")
patient <- apply(delta[[1]]+beta[[1]],2,quantile,probs=c(limits[1],0.5,limits[2]))
namP <- colnames(mdl[[5]]$XF_Pts)
dat <- data.frame(
low = patient[1,],
med = patient[2,],
high = patient[3,],
i = ordered(namP,level=namP[length(namP):1])
)
ans <- ggplot(dat,aes(x=i,ymin=low,ymax=high,y=med))+
geom_pointrange()+
coord_flip()+xlab("coefficients")+ylab("")
}else if(who=="delta"){
delta <- extract(mdl[[2]],pars="b_Delta")
patient <- apply(delta[[1]],2,quantile,probs=c(limits[1],0.5,limits[2]))
namP <- colnames(mdl[[5]]$XF_Pts)
dat <- data.frame(
low = patient[1,],
med = patient[2,],
high = patient[3,],
i = ordered(namP,level=namP[length(namP):1])
)
ans <- ggplot(dat,aes(x=i,ymin=low,ymax=high,y=med))+
geom_pointrange()+
coord_flip()+xlab("coefficients")+ylab("")
}
} else if (type == "area") {
if(requireNamespace("reshape2", quietly = TRUE)){
if(who == "both"){
beta <- extract(mdl[[2]],pars="b_Ctrl")
controls <- reshape2::melt(beta[[1]])
namC <- colnames(mdl[[5]]$XF_Ctrl)
controls$Var2 <- factor(controls$Var2)
levels(controls$Var2) <- namC
controls$Var2 <- ordered(controls$Var2,level=namC[length(namC):1])
controls$group = "Controls"
delta <- extract(mdl[[2]],pars="b_Delta")
tmp <- delta[[1]]+beta[[1]]
patient <- reshape2::melt(tmp)
namP <- colnames(mdl[[5]]$XF_Pts)
patient$Var2 <- factor(patient$Var2)
levels(patient$Var2) <- namP
patient$Var2 <- ordered(patient$Var2,level=namP[length(namP):1])
patient$group = "Patient"
dat <- rbind(
patient,controls
)
ans <- ggplot(dat,aes(x=value,colour=group))+
geom_density()+
facet_grid(Var2~.)+
xlab("coefficients")+ylab("")
}else if(who=="control"){
beta <- extract(mdl[[2]],pars="b_Ctrl")
controls <- reshape2::melt(beta[[1]])
namC <- colnames(mdl[[5]]$XF_Ctrl)
controls$Var2 <- factor(controls$Var2)
levels(controls$Var2) <- namC
controls$Var2 <- ordered(controls$Var2,level=namC[length(namC):1])
dat <- controls
ans <- ggplot(dat,aes(x=value))+
geom_density()+
facet_grid(Var2~.)+
xlab("coefficients")+ylab("")
}else if(who=="patient"){
beta <- extract(mdl[[2]],pars="b_Ctrl")
delta <- extract(mdl[[2]],pars="b_Delta")
tmp <- delta[[1]]+beta[[1]]
patient <- reshape2::melt(tmp)
namP <- colnames(mdl[[5]]$XF_Pts)
patient$Var2 <- factor(patient$Var2)
levels(patient$Var2) <- namP
patient$Var2 <- ordered(patient$Var2,level=namP[length(namP):1])
dat <- patient
ans <- ggplot(dat,aes(x=value))+
geom_density()+
facet_grid(Var2~.)+
xlab("coefficients")+ylab("")
}else if(who=="delta"){
delta <- extract(mdl[[2]],pars="b_Delta")
patient <- reshape2::melt(delta[[1]])
namP <- colnames(mdl[[5]]$XF_Pts)
patient$Var2 <- factor(patient$Var2)
levels(patient$Var2) <- namP
patient$Var2 <- ordered(patient$Var2,level=namP[length(namP):1])
dat <- patient
ans <- ggplot(dat,aes(x=value))+
geom_density()+
facet_grid(Var2~.)+
xlab("coefficients")+ylab("")
}
}
} else if (type == "hist") {
if(requireNamespace("reshape2", quietly = TRUE)){
if(who == "both"){
beta <- extract(mdl[[2]],pars="b_Ctrl")
controls <- reshape2::melt(beta[[1]])
namC <- colnames(mdl[[5]]$XF_Ctrl)
controls$Var2 <- factor(controls$Var2)
levels(controls$Var2) <- namC
controls$Var2 <- ordered(controls$Var2,level=namC[length(namC):1])
controls$group = "Controls"
delta <- extract(mdl[[2]],pars="b_Delta")
tmp <- delta[[1]]+beta[[1]]
patient <- reshape2::melt(tmp)
namP <- colnames(mdl[[5]]$XF_Pts)
patient$Var2 <- factor(patient$Var2)
levels(patient$Var2) <- namP
patient$Var2 <- ordered(patient$Var2,level=namP[length(namP):1])
patient$group = "Patient"
dat <- rbind(
patient,controls
)
ans <- ggplot(dat,aes(x=value,fill=group))+
geom_histogram()+
facet_grid(Var2~.)+
xlab("coefficients")+ylab("")
}else if(who=="control"){
beta <- extract(mdl[[2]],pars="b_Ctrl")
controls <- reshape2::melt(beta[[1]])
namC <- colnames(mdl[[5]]$XF_Ctrl)
controls$Var2 <- factor(controls$Var2)
levels(controls$Var2) <- namC
controls$Var2 <- ordered(controls$Var2,level=namC[length(namC):1])
dat <- controls
ans <- ggplot(dat,aes(x=value))+
geom_histogram()+
facet_grid(Var2~.)+
xlab("coefficients")+ylab("")
}else if(who=="patient"){
beta <- extract(mdl[[2]],pars="b_Ctrl")
delta <- extract(mdl[[2]],pars="b_Delta")
tmp <- delta[[1]]+beta[[1]]
patient <- reshape2::melt(tmp)
namP <- colnames(mdl[[5]]$XF_Pts)
patient$Var2 <- factor(patient$Var2)
levels(patient$Var2) <- namP
patient$Var2 <- ordered(patient$Var2,level=namP[length(namP):1])
dat <- patient
ans <- ggplot(dat,aes(x=value))+
geom_histogram()+
facet_grid(Var2~.)+
xlab("coefficients")+ylab("")
}else if(who=="delta"){
delta <- extract(mdl[[2]],pars="b_Delta")
patient <- reshape2::melt(delta[[1]])
namP <- colnames(mdl[[5]]$XF_Pts)
patient$Var2 <- factor(patient$Var2)
levels(patient$Var2) <- namP
patient$Var2 <- ordered(patient$Var2,level=namP[length(namP):1])
dat <- patient
ans <- ggplot(dat,aes(x=value))+
geom_histogram()+
facet_grid(Var2~.)+
xlab("coefficients")+ylab("")
}
}
}
return(ans)
}
#' Computes log marginal likelihood via bridge sampling.
#'
#' @param object a BMSC object
#' @return an "psis_loo" "loo" object
#' @export
bridge_sampler.BMSC = function(object, ...){
if (requireNamespace("bridgesampling", quietly = TRUE)) {
ans = bridgesampling::bridge_sampler(object[[2]],...)
}
return(ans)
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.