R/BP_FitBayesianCompactness.R

In BoneProfileR: Tools to Study Bone Compactness

#' BP_FitBayesianCompactness estimates Bayesian model of a bone section
#' @title Estimation of Bayesian model of a bone section
#' @author Marc Girondot \email{marc.girondot@@gmail.com}
#' @return The -Ln L
#' @param bone The bone image to be used
#' @param priors Priors 
#' @param n.iter Number of iterations
#' @param n.chains Number of chains
#' @param n.adapt Number of iteration to adapt
#' @param thin Thin parameter for analysis
#' @param analysis Name or rank of analysis
#' @param silent Should some information must me shown ?
#' @description Estimation of Bayesian model of a bone section.
#' @family BoneProfileR
#' @examples
#' \dontrun{
#' # Not run:
#'  library(BoneProfileR)
#'  path_Hedgehog <- system.file("extdata", "Erinaceus_europaeus_fem_2-1_small.png", 
#'                              package = "BoneProfileR")
#'  bone <- BP_OpenImage(file=path_Hedgehog)
#'  bone <- BP_DetectBackground(bone=bone, analysis="logistic")
#'  bone <- BP_DetectForeground(bone=bone, analysis="logistic")
#'  bone <- BP_DetectCenters(bone=bone, analysis="logistic")
#'  bone <- BP_EstimateCompactness(bone, analysis="logistic")
#'  bone <- BP_FitMLCompactness(bone, analysis="logistic")
#'  plot(bone)
#'  plot(bone, type="observations")
#'  plot(bone, type="observations+model", analysis=1)
#'  fittedpar <- BP_GetFittedParameters(bone, analysis="logistic")
#'  bone <- BP_DuplicateAnalysis(bone, from="logistic", to="flexit")
#'  bone <- BP_FitMLCompactness(bone, 
#'                 fitted.parameters=c(fittedpar, K1=1, K2=1), 
#'                 fixed.parameters=NULL, analysis="flexit")
#'  compare_AIC(Logistic=BP_GetFittedParameters(bone, analysis="logistic", alloptim=TRUE), 
#'              Flexit=BP_GetFittedParameters(bone, analysis="flexit", alloptim=TRUE))
#'  out4p <- plot(bone, type="observations+model", analysis="logistic")
#'  out6p <- plot(bone, type="observations+model", analysis="flexit")
#'  bone <- BP_FitBayesianCompactness(bone, analysis="logistic")
#'  plot(bone, type="observations+model", CI="MCMC", analysis="logistic")
#'  bone <- BP_FitBayesianCompactness(bone, analysis="flexit")
#'  plot(bone, type="observations+model", CI="MCMC", analysis="flexit")
#' }
#' @export


BP_FitBayesianCompactness <- function(bone=stop("A result from BP_FitMLCompactness() must be provided"),
                                      priors=NULL, 
                                      n.iter = 10000,
                                      n.chains = 1,
                                      n.adapt = 100,
                                      thin = 1, analysis=1, silent=TRUE) {
  
  # priors=NULL; n.iter = 10000; n.chains = 1; n.adapt = 100; thin = 1, analysis=1
  
  if (is.null(RM_get(x=bone, RMname=analysis, valuename = "optim"))) stop("The model must be first fitted with BP_FitMLCompactness()")
  
  if (is.null(priors)) {
    priors <- data.frame(Density=character(), 
                         Prior1=numeric(), 
                         Prior2=numeric(), 
                         SDProp=numeric(), 
                         Min=numeric(), 
                         Max=numeric(), 
                         Init=numeric(), stringsAsFactors = FALSE)
    p <- BP_GetFittedParameters(bone, analysis = analysis, alloptim = FALSE) 

    if (!is.na(p["P"])) {
      priors <- rbind(priors, data.frame(Density="dunif", 
                                         Prior1=0, 
                                         Prior2=1, 
                                         SDProp=0.2, 
                                         Min=0, 
                                         Max=1, 
                                         Init=unname(p["P"]), stringsAsFactors = FALSE, 
                                         row.names = "P"))
    }
    if (!is.na(p["S"])) {
      priors <- rbind(priors, data.frame(Density="dunif", 
                                         Prior1=0, 
                                         Prior2=10, 
                                         SDProp=0.2, 
                                         Min=0, 
                                         Max=10, 
                                         Init=unname(p["S"]), stringsAsFactors = FALSE, 
                                         row.names = "S"))
    }
    if (!is.na(p["Min"])) {
      priors <- rbind(priors, data.frame(Density="dunif", 
                                         Prior1=0, 
                                         Prior2=0.8, 
                                         SDProp=0.2, 
                                         Min=0, 
                                         Max=0.8, 
                                         Init=unname(p["Min"]), stringsAsFactors = FALSE, 
                                         row.names = "Min"))
    }
    if (!is.na(p["Max"])) {
      priors <- rbind(priors, data.frame(Density="dunif", 
                                         Prior1=0.2, 
                                         Prior2=1, 
                                         SDProp=0.2, 
                                         Min=0.2, 
                                         Max=1, 
                                         Init=unname(p["Max"]), stringsAsFactors = FALSE, 
                                         row.names = "Max"))
    }
    if (!is.na(p["K1"])) {
      priors <- rbind(priors, data.frame(Density="dunif", 
                                         Prior1=-10, 
                                         Prior2=10, 
                                         SDProp=0.2, 
                                         Min=-10, 
                                         Max=10, 
                                         Init=unname(p["K1"]), stringsAsFactors = FALSE, 
                                         row.names = "K1"))
    }
    if (!is.na(p["K2"])) {
      priors <- rbind(priors, data.frame(Density="dunif", 
                                         Prior1=-10, 
                                         Prior2=10, 
                                         SDProp=0.2, 
                                         Min=-10, 
                                         Max=10, 
                                         Init=unname(p["K2"]), stringsAsFactors = FALSE, 
                                         row.names = "K2"))
    }
    if (!silent) priors
  }
  
  fixedpar <- BP_GetFittedParameters(bone, analysis = analysis, alloptim = TRUE)$fixed.parameters
  
  mcmc <- MHalgoGen(
    likelihood = BP_LnLCompactness,
    bone=bone, 
    fixed.parameters=fixedpar, 
    parameters_name = "par", 
    parameters = priors, n.iter = n.iter,
    n.chains = n.chains,
    n.adapt = n.adapt,
    thin = thin, adaptive = TRUE)
  
  data <- RM_get(x=bone, RMname=analysis, valuename = "compactness.synthesis")
  outmcmc <- matrix(NA, ncol=nrow(data), nrow=n.iter)
  
  for (iter in 1:n.iter) {
    p <- c(mcmc$resultMCMC[[1]][iter, ], 
           fixedpar)
    
    
    Min <- p["Min"]
    Max <- p["Max"]
    
    # 21/2/2020
    p["S"] <- 1/(4*p["S"])
    
    c <- flexit(x = data$distance.center, 
                par = p) * (Max - Min) + Min
    
    outmcmc[iter, ] <- c
  }
  
  qmcmc <- apply(X = outmcmc, MARGIN = 2, FUN = function(x) quantile(x, probs = c(0.025, 0.5, 0.975)))
  colnames(qmcmc) <- data$distance.center
  
  mcmc <- modifyList(mcmc, list(quantiles=qmcmc))
  mcmc$timestamp <- date()
  
  mcmc$resultMCMC[[1]][, "Min"] <- mcmc$resultMCMC[[1]][, "Min"]
  mcmc$resultMCMC[[1]][, "Max"] <- mcmc$resultMCMC[[1]][, "Max"]
  
  summary.table <- data.frame(mean=apply(mcmc$resultMCMC[[1]], MARGIN=2, FUN = mean), 
  se=apply(mcmc$resultMCMC[[1]], MARGIN=2, FUN = sd))
  mcmc <- modifyList(mcmc, list(summary.table=summary.table))
  
  # Je retire bone car ça prend trop de place
  mcmc$parametersMCMC$control$bone <- NULL
  
  bone <- RM_add(x=bone, RMname = analysis, valuename = "mcmc", value=mcmc)
  
  return(bone)
}