BoneProfileR (version 4.0)

BP_FitMLPeriodicCompactness: Estimation of the likelihood of a bone section

Description

Estimation of the compactness of a bone section using radial model.
If the fitted.parameters and fixed.parameters are NULL and the analysis includes a BP_FitMLCompactness() result, the values of this result is used as a reference for fitted.parameters and fixed.parameters.
If no BP_FitMLCompactness() result is available, it will use:
fitted.parameters=c(P=0.5, S=0.05, Min=-2, Max=5); fixed.parameters=c(K1=1, K2=1).
The reference for radial estimation of compactness is the trigonometric circle for rotation.angle=0 in BP_EstimateCompactness():

• The top of the section is located at -pi/2.
  

• The left of the section is located at -pi and +pi.
  

• The bottom of the section is located at pi/2.
  

• The right of the section is 0.
  If rotation.angle is different from 0, the value of rotation.angle is added to the angle modulo 2.pi.
  The two-steps analysis performs first a quasi-Newton method, then a Bayesian MCMC and finally again a quasi-Newton method. It generally ensures that global minimum is found. On the other hand, it doubles the time to complete for each angle.
  To control the parallel computing, use:
  options(mc.cores = [put here the number of cores you want use])
  options(forking = FALSE) or options(forking = TRUE)
  The maximum number of cores is obtained by: parallel::detectCores()
  

Usage

BP_FitMLPeriodicCompactness(
  bone,
  fitted.parameters = NULL,
  fixed.parameters = NULL,
  analysis = 1,
  silent = FALSE,
  replicates.CI = NULL,
  twosteps = FALSE,
  amplitude.max = 0.1,
  control.optim = list(trace = 1)
)

Value

The -Ln L

Arguments

bone

The bone image to be used

fitted.parameters

Parameters of the model to be fitted

fixed.parameters

Fixed parameters of the model

analysis

Name or rank of analysis

silent

Should the function displays some information?

replicates.CI

Number of replicates to estimate confidence interval using Hessian

twosteps

Should a 2-steps analysis be performed? It can be sometimes useful.

amplitude.max

The maximum allowed amplitude for each parameter

control.optim

The list of options for optim.

Author

Marc Girondot

Details

BP_FitMLPeriodicCompactness estimates likelihood of global model of a bone section

See Also

Other BoneProfileR: BP_AutoFit(), BP_ChooseBackground(), BP_ChooseCenter(), BP_ChooseForeground(), BP_DetectBackground(), BP_DetectCenters(), BP_DetectForeground(), BP_DuplicateAnalysis(), BP_EstimateCompactness(), BP_FitBayesianCompactness(), BP_FitBayesianPeriodicCompactness(), BP_FitMLCompactness(), BP_FitMLRadialCompactness(), BP_GetFittedParameters(), BP_ListAnalyses(), BP_LnLCompactness(), BP_OpenImage(), BP_Report(), Erinaceus_europaeus, plot.BoneProfileR(), summary.BoneProfileR()

Examples

if (FALSE) {
# 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", cut.angle = 60)
 bone <- BP_FitMLCompactness(bone, analysis="logistic", twosteps=TRUE)
 plot(bone, type="observations+model", analysis="logistic")
 par <- BP_GetFittedParameters(bone, analysis="logistic", ML=TRUE, return.all=FALSE)[, "mean"]
 options(mc.cores=parallel::detectCores())
 
 #############################################
 # Periodic analysis
 #############################################
 bone <- BP_FitMLPeriodicCompactness(bone, analysis="logistic", control.optim=list(trace=2), 
                                     fitted.parameters=c(par, PSin=0.001, PCos=0.001, 
                                     SSin=0.001, SCos=0.001, MinSin=0.001, MinCos=0.001, 
                                     MaxSin=0.001, MaxCos=0.001), replicates.CI=2000)
 analysisP <- BP_GetFittedParameters(bone, analysis="logistic", type="periodic", 
                                     ML=TRUE, return.all=FALSE)[, "mean"]
 analysisP$par                                    
 plot(bone, type="periodic", parameter.name="compactness", col=rainbow(128))
 plot(bone, type="periodic", parameter.name="compactness", 
               col=hcl.colors(12, "YlOrRd", rev = TRUE))
 plot(bone, type="periodic", parameter.name="averagemodel")
 plot(bone, type="periodic", parameter.name="P", 
               rgb(red = 0.7, green = 0.7, blue = 0.7, alpha = 0.2))
 plot(bone, type="periodic", parameter.name="P", ylim=c(0, 1), 
               rgb(red = 0.7, green = 0.7, blue = 0.7, alpha = 0.2))
 boneNoPeriodic <- BP_FitMLPeriodicCompactness(bone, analysis="logistic", 
                                               control.optim=list(trace=2), 
                                     fitted.parameters=par, replicates.CI=2000)
 analysisNP <- BP_GetFittedParameters(boneNoPeriodic, analysis="logistic", ML=TRUE, 
                                      return.all=TRUE, type="periodic")
 analysisNP$par
 compare_AIC(PeriodicModel=analysisP, 
             NoPeriodicModel=analysisNP)
 
 #############################################
 
 # Note that the absolute likelihood is dependent on the number of angle cut
 # Only models analyzed with the same number of angle cuts can be compared
 
 dbinom(5, 10, prob=0.4, log=TRUE); 
       dbinom(2, 5, prob=0.4, log=TRUE)+dbinom(3, 5, prob=0.4, log=TRUE)
 # But the likelihood difference between two models are not:
 dbinom(5, 10, prob=0.4, log=TRUE)-dbinom(5, 10, prob=0.3, log=TRUE)
 dbinom(2, 5, prob=0.4, log=TRUE)+dbinom(3, 5, prob=0.4, log=TRUE)- 
      dbinom(2, 5, prob=0.3, log=TRUE)-dbinom(3, 5, prob=0.3, log=TRUE)
}