BayesSizeAndShape.jl

Bayesian regression models for size-and-shape data
Author GianlucaMastrantonio
Popularity
0 Stars
Updated Last
12 Months Ago
Started In
January 2023

BayesSizeAndShape

This package implements a Bayesian regression for size-and-shape data, based on "Di Noia, A., Mastrantonio, G., and Jona Lasinio, G., “Bayesian Size-and-Shape regression modelling”, arXiv e-prints, 2023". To install, simply do

julia> ]

pkg> add BayesSizeAndShape

at the julia prompt. At the present moment, the package implements a model for two-dimensional data with reflection information. The function to use is SizeAndShapeMCMC

Basic Usage

In the demo directory, there is a julia file Ex_p2withreflection.jl with an example of how to implement the model, which we will describe also here.

Let

  • n be the number of shapes;
  • k be the number of landmarks (for the pre-form matrix X);
  • p be the dimension of each landmark (only p=2 is implemented)
  • d be the number of covariates to use;

First we load the required packages

using Random, Distributions, LinearAlgebra, StatsBase
using Kronecker, DataFrames, StatsModels, CategoricalArrays
using BayesSizeAndShape

To be able to use the function, we have to simulate some data. We first simulate the regressive coefficients, which are in the reg object - reg must be of dimension (kd, p):

n::Int64 = 100;
p::Int64 = 2;
k::Int64 = 10;
d::Int64 = 3;

# regression
reg::Matrix{Float64} = zeros(Float64, k*d, p);
reg[:] = rand(Normal(0.0, 5.0), prod(size(reg)));

The regressive coefficients are standardized with the function standardize_reg, for identifiability purposes. The second argument is here used to specify the dimension of each landmark, Value2() is used for two-dimensional data, while Valuep3() (not yet implemented) is for three-dimensional data:

standardize_reg(reg::Matrix{Float64}, Valuep2())

A DataFrame named zmat, containing the covariates is simulated, and the function compute_designmatrix is used to compute the design matrix - zmat must be of dimension (n, d):

zmat = DataFrame(
    x1 = rand(Normal(10.0,1.0 ),n),
    x2 = sample(["A", "B"],n)
)
zmat[:,1] = (zmat[:,1] .- mean(zmat[:,1])) ./ std(zmat[:,1])
zmat.x2 = categorical(zmat.x2)

zmat_modmat_ModelFrame = ModelFrame(@formula(1 ~ 1+x1+x2), zmat);
zmat_modmat = ModelMatrix(zmat_modmat_ModelFrame).m
design_matrix = compute_designmatrix(zmat_modmat, k);

Please notice that, since in the MCMC algorithm the predictors are standardized, we standardize them before simulating the data.

The pre-form matrix X, here saved in the object dataset_complete, is simulated from a multivariate normal, and its size-and-shape version, contained in the object dataset, is obtained by using the function compute_ss_from_pre. The third argument of the function is used to specify if reflection must be kept (true is the only available option in this implementation)

sigma::Symmetric{Float64,Matrix{Float64}} = Symmetric(rand(InverseWishart(k + 2, 5.0 * Matrix{Float64}(I, k, k))));
dataset_complete = zeros(Float64,k,p,n);
dataset = zeros(Float64, k, p, n);
for i_n = 1:n
    for i_p = 1:p
        dataset_complete[:, i_p, i_n] = rand(MvNormal(design_matrix[:, :, i_n] * reg[:, i_p], sigma))
        
    end
end
rmat = compute_ss_from_pre(dataset_complete, dataset, true);

If instead of the preform you want to simulate a landmark matrix (or if you have landmark data), you can use the function compute_helmertized_configuration to obtain the preforms from the landmark.

Posterior samples are obtained with the function SizeAndShapeMCMC_p2withreflection. To specify the regressive formula, you can use @formula, where on the left-hand side there must be 1 and on the right-hand side is the actual regressive formula.

betaout, sigmaout = SizeAndShapeMCMC(;
    dataset = dataset,
    fm = @formula(1 ~ 1 + x1 + x2),
    covariates = zmat,
    iterations = (iter=1000, burnin=200, thin=2),
    betaprior = Normal(0.0, 10000.0),
    sigmaprior = InverseWishart(k + 2, 5.0 * Matrix{Float64}(I, k, k)),
    beta_init = zeros(Float64, k * d, p),
    sigma_init = Symmetric(Matrix{Float64}(I, k, k)),
    rmat_init = reshape(vcat([Matrix{Float64}(I, p, p)[:] for i = 1:n]...), (p, p, n)),
    reflection = true
);

The objects betaout and sigmaout contain the posterior samples of the regressive coefficients and covariance matrix, respectively.

Citing

See CITATION.bib