Example 3: Multivariate generalized linear mixed models for continuous bounded outcomes
Henrique A. Laureano, Ricardo R. Petterle and Wagner H. Bonat - LEG/UFPR
1 Loading extra packages
library(knitr)
library(TMB)
library(betareg)
library(ggplot2)
library(tidyverse)
library(cowplot)
library(ggcorrplot)2 Body fat percentage data set
The motivating example analyzed in this article uses data of the body fat percentage of subjects assisted at the Endocrinology and Metabology Service of the Clinic’s Hospital of the Paraná Federal University, Curitiba, Brazil. The data set contains 298 observations and was previously analyzed by Petterle et al. (2020) by using the multivariate quasi-beta regression model. The data set is composed by two continuous and two categorical covariates, which correspond to: - subject age (in years, ranging from 18 to 87 years), - body mass index (in kg/m\(^2\), ranging from 18.5 to 29.9 kg/m\(^2\)), - gender of subjects (F: female; M: Male) and - level of physical activity practiced routinely (S: sedentary; IA: insufficiently active; A: active). The response variables refer to the body fat percentage (BFP) measured at five regions of the body.
data3 <- read.table("Data_Set.csv", header = TRUE, sep = ",", na.strings = "")tail(data3, n = 10)## AGE BMI SEX IPAQ ARMS LEGS TRUNK ANDROID GYNOID
## 289 79 29.4 Male Sedentary 0.367 0.444 0.452 0.473 0.479
## 290 80 28.3 Male Active 0.279 0.353 0.392 0.432 0.407
## 291 80 25.6 Male Insufficiently active 0.229 0.308 0.378 0.437 0.353
## 292 82 28.8 Female Sedentary 0.471 0.470 0.516 0.548 0.538
## 293 82 22.6 Female Sedentary 0.361 0.473 0.390 0.408 0.491
## 294 84 29.3 Male Insufficiently active 0.227 0.248 0.405 0.464 0.334
## 295 86 27.2 Male Sedentary 0.295 0.302 0.416 0.481 0.383
## 296 86 22.6 Male Sedentary 0.276 0.263 0.370 0.436 0.290
## 297 86 25.6 Male Sedentary 0.236 0.300 0.392 0.475 0.338
## 298 87 23.4 Male Sedentary 0.267 0.310 0.373 0.410 0.3943 Descriptive analysis
# Preparing
par(mfrow = c(4,5), mar=c(2.8, 2.6, 1.2, 0.5), mgp = c(1.6, 0.6, 0))
# cex (theme_cowplot)
fig_cex = 12
# Age
p1 <- data3 %>%
ggplot(aes(x = AGE, y = ARMS)) +
geom_point(size = 1.1, alpha = 0.8, color = "gray80") +
xlab("Age (years)") +
ylab("Arms") +
geom_smooth(method = lm, color = "gray70",
formula = y ~ splines::bs(x, 3),
se = T) +
theme_cowplot(fig_cex)
p2 <- data3 %>%
ggplot(aes(x = AGE, y = LEGS)) +
geom_point(size = 1.1, alpha = 0.8, color = "gray80") +
xlab("Age (years)") +
ylab("Legs") +
geom_smooth(method = lm, color = "gray70",
formula = y ~ splines::bs(x, 3),
se = T) +
theme_cowplot(fig_cex)
p3 <- data3 %>%
ggplot(aes(x = AGE, y = TRUNK)) +
geom_point(size = 1.1, alpha = 0.8, color = "gray80") +
xlab("Age (years)") +
ylab("Trunk") +
geom_smooth(method = lm, color = "gray70",
formula = y ~ splines::bs(x, 3),
se = T) +
theme_cowplot(fig_cex)
p4 <- data3 %>%
ggplot(aes(x = AGE, y = ANDROID)) +
geom_point(size = 1.1, alpha = 0.8, color = "gray80") +
xlab("Age (years)") +
ylab("Android") +
geom_smooth(method = lm, color = "gray70",
formula = y ~ splines::bs(x, 3),
se = T) +
theme_cowplot(fig_cex)
p5 <- data3 %>%
ggplot(aes(x = AGE, y = GYNOID)) +
geom_point(size = 1.1, alpha = 0.8, color = "gray80") +
xlab("Age (years)") +
ylab("Gynoid") +
geom_smooth(method = lm, color = "gray70",
formula = y ~ splines::bs(x, 3),
se = T) +
theme_cowplot(fig_cex)
# BMI
p6 <- data3 %>%
ggplot(aes(x = BMI, y = ARMS)) +
geom_point(size = 1.1, alpha = 0.8, color = "gray80") +
xlab(expression(paste("BMI (kg/m"^2,")"))) +
ylab("Arms") +
geom_smooth(method = lm, color = "gray70",
formula = y ~ splines::bs(x, 3),
se = T) +
theme_cowplot(fig_cex)
p7 <- data3 %>%
ggplot(aes(x = BMI, y = LEGS)) +
geom_point(size = 1.1, alpha = 0.8, color = "gray80") +
xlab(expression(paste("BMI (kg/m"^2,")"))) +
ylab("Legs") +
geom_smooth(method = lm, color = "gray70",
formula = y ~ splines::bs(x, 3),
se = T) +
theme_cowplot(fig_cex)
p8 <- data3 %>%
ggplot(aes(x = BMI, y = TRUNK)) +
geom_point(size = 1.1, alpha = 0.8, color = "gray80") +
xlab(expression(paste("BMI (kg/m"^2,")"))) +
ylab("Trunk") +
geom_smooth(method = lm, color = "gray70",
formula = y ~ splines::bs(x, 3),
se = T) +
theme_cowplot(fig_cex)
p9 <- data3 %>%
ggplot(aes(x = BMI, y = ANDROID)) +
geom_point(size = 1.1, alpha = 0.8, color = "gray80") +
xlab(expression(paste("BMI (kg/m"^2,")"))) +
ylab("Android") +
geom_smooth(method = lm, color = "gray70",
formula = y ~ splines::bs(x, 3),
se = T) +
theme_cowplot(fig_cex)
p10 <- data3 %>%
ggplot(aes(x = BMI, y = GYNOID)) +
geom_point(size = 1.1, alpha = 0.8, color = "gray80") +
xlab(expression(paste("BMI (kg/m"^2,")"))) +
ylab("Gynoid") +
geom_smooth(method = lm, color = "gray70",
formula = y ~ splines::bs(x, 3),
se = T) +
theme_cowplot(fig_cex)
# Gender
p11 <- data3 %>%
ggplot(aes(x = SEX, y = ARMS)) +
geom_boxplot(fill = "gray80") +
xlab("Gender") +
ylab("Arms") +
theme_cowplot(fig_cex)
p12 <- data3 %>%
ggplot(aes(x = SEX, y = LEGS)) +
geom_boxplot(fill = "gray80") +
xlab("Gender") +
ylab("Legs") +
theme_cowplot(fig_cex)
p13 <- data3 %>%
ggplot(aes(x = SEX, y = TRUNK)) +
geom_boxplot(fill = "gray80") +
xlab("Gender") +
ylab("Trunk") +
theme_cowplot(fig_cex)
p14 <- data3 %>%
ggplot(aes(x = SEX, y = ANDROID)) +
geom_boxplot(fill = "gray80") +
xlab("Gender") +
ylab("Android") +
theme_cowplot(fig_cex)
p15 <- data3 %>%
ggplot(aes(x = SEX, y = GYNOID)) +
geom_boxplot(fill = "gray80") +
xlab("Gender") +
ylab("Gynoid") +
theme_cowplot(fig_cex)
# IPAQ
data3$IPAQ <- factor(data3$IPAQ)
levels(data3$IPAQ) <- c("S", "IA", "A")
p16 <- data3 %>%
ggplot(aes(x = IPAQ, y = ARMS)) +
geom_boxplot(fill = "gray80") +
xlab("IPAQ") +
ylab("Arms") +
theme_cowplot(fig_cex)
p17 <- data3 %>%
ggplot(aes(x = IPAQ, y = LEGS)) +
geom_boxplot(fill = "gray80") +
xlab("IPAQ") +
ylab("Legs") +
theme_cowplot(fig_cex)
p18 <- data3 %>%
ggplot(aes(x = IPAQ, y = TRUNK)) +
geom_boxplot(fill = "gray80") +
xlab("IPAQ") +
ylab("Trunk") +
theme_cowplot(fig_cex)
p19 <- data3 %>%
ggplot(aes(x = IPAQ, y = ANDROID)) +
geom_boxplot(fill = "gray80") +
xlab("IPAQ") +
ylab("Android") +
theme_cowplot(fig_cex)
p20 <- data3 %>%
ggplot(aes(x = IPAQ, y = GYNOID)) +
geom_boxplot(fill = "gray80") +
xlab("IPAQ") +
ylab("Gynoid") +
theme_cowplot(fig_cex)
# Figure
plot_grid(p1, p2, p3, p4, p5,
p6, p7, p8, p9, p10,
p11, p12, p13, p14, p15,
p16, p17, p18, p19, p20,
label_size = 13,
labels = c("(A)", "(B)", "(C)", "(D)", "(E)",
"(F)", "(G)", "(H)", "(I)", "(J)",
"(K)", "(L)", "(M)", "(N)", "(O)",
"(P)", "(Q)", "(R)", "(S)", "(T)"),
hjust = -0.05, vjust = 1,
ncol = 5, nrow = 4)
4 Multivariate generalized linear mixed model
4.1 Beta distribution
Let \(Y_1, \ldots, Y_n\) be independent random variables with probability density function given by \[ f(y;\mu,\phi) = \dfrac{\Gamma(\phi)}{\Gamma(\mu \phi) \Gamma[(1 - \mu)\phi]} y^{\mu \phi - 1} (1 - y)^{(1 - \mu) \phi - 1}, \] where \(y, \mu \in (0,1)\) and \(\phi > 0\) denote the precision parameter. According to Ferrari and Cribari-Neto (2004) the expectation and variance of the beta distribution are given, respectively, by
\[\begin{equation} \label{betamoments} \mathrm{E}(Y) = \mu \quad~\mathrm{and}~\quad \mathrm{Var}(Y) = \dfrac{\mu(1-\mu)}{1 + \phi}. \end{equation}\]
4.2 Model specification
Let \(\mathbf{Y}\) be a response variable matrix \(N \times R\), where rows \(i = 1, \ldots, N\) are independent subjects and columns \(r = 1, \ldots, R\) refers to the \(R\) response variables. Let \(y_{ir} \in (0,1)\) be the observed value of a random variable \(Y_{ir}\), and let \(\mathbf{u}_i = (u_{i1}, \ldots, u_{iR})\) be the vector of random intercepts for the subject \(i\). For specification of the MGLMM, consider the following hierarchical structure \[\begin{align} \nonumber \label{eq:dist} Y_{ir} | u_{ir} & \sim \mathcal{B}(\mu_{ir}, \phi_{r}), \nonumber \\ g_r(\mu_{ir}) & = \mathbf{x}_{ir}^{\top}\boldsymbol{\beta}_r + u_{ir}, \end{align}\] where \(\phi_{r}\) is assumed to be constant for the \(r^{th}\) response variable. The function \(g(\cdot)\) is a suitable link function for binary data. In this notation, \(\mathbf{x}_{ir}\) and \(\boldsymbol{\beta}_r\) are \(Q_r \times 1\) vectors of known covariates and unknown regression coefficients associated to the expectation of the \(r^{th}\) response variable for the \(i^{th}\) subject. Note that, the index \(r\) was introduced to ensure different set of covariates for each response variable. Moreover, we assume that the random intercepts follow a multivariate Gaussian distribution, \(u_{ir} \sim \mathcal{MVN}_r(\mathbf{0}, \boldsymbol{\Sigma})\), where \(\boldsymbol{\Sigma}\) is the \(R \times R\) unstructured covariance matrix. Thus, the random intercepts distribution is given by \[\begin{equation} \label{eq:gaussmult} \begin{pmatrix} u_{i1} \\ \vdots \\ u_{iR} \\ \end{pmatrix} \sim \mathcal{MVN}_R \left[ \begin{pmatrix} 0 \\ \vdots \\ 0 \\ \end{pmatrix}, \mathop{\boldsymbol{\Sigma}}_{R \times R} = \begin{pmatrix} \sigma_1^2 \hspace{-0.2cm} & \hspace{-0.2cm} \ldots \hspace{-0.2cm} & \hspace{-0.2cm}~\rho_{1R} \sigma_1 \sigma_R \\ \vdots \hspace{-0.1cm} & \hspace{-0.2cm} \ddots \hspace{-0.2cm} & \hspace{-0.2cm} \vdots \\ \rho_{1R} \sigma_1 \sigma_R \hspace{-0.1cm} & \hspace{-0.1cm} \ldots \hspace{-0.2cm} & \hspace{-0.2cm} \sigma_R^2 \\ \end{pmatrix} \right], \end{equation}\] where \(\rho_{RR'} \in [-1, 1]\) and \(\sigma_R^2 > 0\), for \(r, r' = 1,, \ldots, R\), are the correlation and variance of the random intercepts, respectively. We highlight that, in our implementation the covariance matrix of the multivariate Gaussian distribution is specified as follows \[\begin{equation} \label{eq:covTMB} \boldsymbol{\Sigma} = \mathrm{diag}(\sigma_1, \ldots, \sigma_R) \boldsymbol{\Sigma}_{\rho} \mathrm{diag}(\sigma_1, \ldots, \sigma_R), \end{equation}\] where \(\mathrm{diag}(\sigma_1, \ldots, \sigma_R)\) for \(r = 1, \ldots, R\) is the diagonal matrix of the standard deviations of the random intercepts and \(\boldsymbol{\Sigma}_{\rho}\) is the \(R \times R\) unstructured correlation matrix given by \[\begin{equation} \label{eq:corrTMB} \boldsymbol{\Sigma}_{\rho} = \mathbf{D}^{-\tfrac{1}{2}} \mathbf{L} \mathbf{L}^{\top} \mathbf{D}^{-\tfrac{1}{2}}, \end{equation}\] where \(\mathbf{L}\) is a lower-triangular matrix, composed by 1’s on the diagonal and correlation parameters \(\theta_1, \ldots, \theta_{R(R-1)/2}\) in the lower triangle. For \(R = 5\), the lower-triangular matrix is given by \[\begin{equation} \nonumber \mathbf{L} = \begin{bmatrix} 1 & & & & \\ \theta_1 & 1 & & & \\ \theta_2 & \theta_5 & 1 & & \\ \theta_3 & \theta_6 & \theta_8 & 1 & \\ \theta_4 & \theta_7 & \theta_9 & \theta_{10} & 1 \\ \end{bmatrix}. \end{equation}\] Note that, in the correlation matrix \(\mathbf{D} = \mathrm{diag}(\mathbf{L} \mathbf{L}^{\top})\) and \(\mathbf{D}^{-\tfrac{1}{2}}\) is the corresponding inverse square root matrix. This parameterization allows us to estimate the \(R(R-1)/2\) correlation parameters respecting the constraint of positive definite matrix. Furthermore, such parameterization facilitates the estimation of the correlation parameters, since \(\theta_r \in \mathrm{R}\) for \(r = 1, \ldots, R(R-1)/2\).
4.3 Estimation and Inference
Let \(\boldsymbol{\theta} = (\boldsymbol{\beta}^{\top}, \boldsymbol{\phi}^{\top}, \boldsymbol{\tau}^{\top})^{\top}\) be a vector composed of three sets of parameters, where \(\boldsymbol{\beta} = (\boldsymbol{\beta}_1^{\top}, \ldots, \boldsymbol{\beta}_R^{\top})^{\top}\), \(\boldsymbol{\phi} = (\phi_1, \ldots, \phi_R)^{\top}\) and \(\boldsymbol{\tau} = (\sigma_1, \ldots, \sigma_R, \rho_{1R'}, \ldots, \rho_{RR'})^{\top}\) are vectors \(Q \times 1\), \(R \times 1\) and \(J \times 1\) of parameters associated with the regression, precision and covariance structures, respectively. In this notation, \(Q = \sum_r Q_r\) correspond to all regression parameters and \(\boldsymbol{\tau}\) has dimension \(J = R + R(R-1)/2\), where \(R\) and \(R(R - 1)/2\) are the number of standard deviations and correlation parameters of the covariance matrix, respectively.
Let \(\mathbf{y}_i = (y_{i1}, \ldots, y_{iR})\) denote the \(i^{th}\) row of the response variable matrix and let \(\mathbf{u}_i = (u_{i1}, \ldots, u_{iR})\) denote the corresponding vector of correlated random intercepts, for \(i = 1, \ldots, N\) and \(r = 1, \ldots, R\). For each independent subject, the marginal likelihood is given by \[\begin{equation} \label{eq:log-Lik} \mathcal{L}(\boldsymbol{\theta} | \mathbf{y}_i) = \int_{\mathrm{R}^r} \bigg\{ \prod_{r = 1}^{R} f_{r}(y_{ir} | \mathbf{u}_i, \boldsymbol{\beta}_r, \phi_r) \bigg\} f(\mathbf{u}_i | \boldsymbol{\Sigma}) d\mathbf{u}_i, \end{equation}\] where \(f_{r}(y_{ir} | \mathbf{u}_{i}, \boldsymbol{\beta}_r, \phi_r)\) denote the \(r^{th}\) conditional distribution following the hierarchical specification and \(f(\mathbf{u}_{i} | \boldsymbol{\Sigma})\) denote the multivariate Gaussian distribution of the random intercepts.
We assume independence between the \(N\) subjects. Thus, the full marginal likelihood function for the beta probability density function is given by \[\begin{align} \label{eq:log-LikBeta} \mathcal{L}(\boldsymbol{\theta} | \mathbf{y}) = \prod_{i = 1}^{N} \int_{\mathrm{R}^r} & \bigg\{ \prod_{r = 1}^{R} \dfrac{ \Gamma(\phi_r) }{\Gamma(\mu_{ir}\phi_r) \Gamma((1-\mu_{ir})\phi_r) } y_{ir}^{\mu_{ir} \phi_r-1} (1 - y_{ir})^{\phi_r (1 - \mu_{ir}) - 1 } \bigg\} \nonumber \\ & \times (2\pi)^{-R/2} |\boldsymbol{\Sigma}|^{-1/2} \exp \bigg( -\dfrac{1}{2} \mathbf{u}_i^{\top}\boldsymbol{\Sigma}^{-1} \mathbf{u}_i \bigg) d\mathbf{u}_i \end{align}\]
For details, see Petterle et al. (2021)
5 Template Model Builder (TMB) package
For details about TMB package see Kristensen et al. (2016).
5.1 Code written in C++ template
// MGLMM Beta
#include <TMB.hpp>
template<class Type>
Type objective_function<Type>::operator() ()
{
using namespace density;
DATA_VECTOR(Y1);
DATA_VECTOR(Y2);
DATA_VECTOR(Y3);
DATA_VECTOR(Y4);
DATA_VECTOR(Y5);
DATA_MATRIX(X);
PARAMETER_VECTOR(beta1);
PARAMETER_VECTOR(beta2);
PARAMETER_VECTOR(beta3);
PARAMETER_VECTOR(beta4);
PARAMETER_VECTOR(beta5);
PARAMETER_MATRIX(U);
PARAMETER_VECTOR(rho);
PARAMETER_VECTOR(sigma);
PARAMETER_VECTOR(phi);
// Preparing
vector<Type> rho_temp(10);
rho_temp = rho;
vector<Type> sigma_temp(5);
sigma_temp = sigma;
vector<Type> mu1(Y1.size());
vector<Type> mu2(Y2.size());
vector<Type> mu3(Y3.size());
vector<Type> mu4(Y4.size());
vector<Type> mu5(Y5.size());
vector<Type> phi_temp(5);
phi_temp = phi;
vector<Type> shapeA1(Y1.size());
vector<Type> shapeA2(Y2.size());
vector<Type> shapeA3(Y3.size());
vector<Type> shapeA4(Y4.size());
vector<Type> shapeA5(Y5.size());
vector<Type> shapeB1(Y1.size());
vector<Type> shapeB2(Y2.size());
vector<Type> shapeB3(Y3.size());
vector<Type> shapeB4(Y4.size());
vector<Type> shapeB5(Y5.size());
// Linear predictor for mean
mu1 = exp(X*beta1 + U.col(0).array())/(1 + exp(X*beta1 + U.col(0).array()));
mu2 = exp(X*beta2 + U.col(1).array())/(1 + exp(X*beta2 + U.col(1).array()));
mu3 = exp(X*beta3 + U.col(2).array())/(1 + exp(X*beta3 + U.col(2).array()));
mu4 = exp(X*beta4 + U.col(3).array())/(1 + exp(X*beta4 + U.col(3).array()));
mu5 = exp(X*beta5 + U.col(4).array())/(1 + exp(X*beta5 + U.col(4).array()));
// Shape A (Beta)
shapeA1 = mu1*exp(phi(0));
shapeA2 = mu2*exp(phi(1));
shapeA3 = mu3*exp(phi(2));
shapeA4 = mu4*exp(phi(3));
shapeA5 = mu5*exp(phi(4));
// Shape B (Beta)
shapeB1 = (1-mu1)*exp(phi(0));
shapeB2 = (1-mu2)*exp(phi(1));
shapeB3 = (1-mu3)*exp(phi(2));
shapeB4 = (1-mu4)*exp(phi(3));
shapeB5 = (1-mu5)*exp(phi(4));
// Full Log-likelihood
Type nll1 = 0;
for(int i=0; i<Y1.size(); i++)
nll1 -= dbeta(Y1(i), shapeA1(i), shapeB1(i), true);
Type nll2 = 0;
for(int i=0; i<Y2.size(); i++)
nll2 -= dbeta(Y2(i), shapeA2(i), shapeB2(i), true);
Type nll3 = 0;
for(int i=0; i<Y3.size(); i++)
nll3 -= dbeta(Y3(i), shapeA3(i), shapeB3(i), true);
Type nll4 = 0;
for(int i=0; i<Y4.size(); i++)
nll4 -= dbeta(Y4(i), shapeA4(i), shapeB4(i), true);
Type nll5 = 0;
for(int i=0; i<Y5.size(); i++)
nll5 -= dbeta(Y5(i), shapeA5(i), shapeB5(i), true);
Type nll = 0;
for(int i = 0; i < Y1.size(); i++)
nll += VECSCALE(UNSTRUCTURED_CORR(rho), sigma)(U.row(i));
matrix<Type> Cor(5,5);
Cor = UNSTRUCTURED_CORR(rho).cov();
REPORT(Cor);
return nll1 + nll2 + nll3 + nll4 + nll5 + nll;
}5.2 Univariate beta regression model
## Response 1
fit1 <- betareg::betareg(ARMS ~ AGE + BMI + SEX + IPAQ,
data = data3)
## Response 2
fit2 <- betareg::betareg(LEGS ~ AGE + BMI + SEX + IPAQ,
data = data3)
## Response 3
fit3 <- betareg::betareg(TRUNK ~ AGE + BMI + SEX + IPAQ,
data = data3)
## Response 4
fit4 <- betareg::betareg(ANDROID ~ AGE + BMI + SEX + IPAQ,
data = data3)
## Response 5
fit5 <- betareg::betareg(GYNOID ~ AGE + BMI + SEX + IPAQ,
data = data3)5.3 Compiling and loading
compile("MGLMM_Beta.cpp")
dyn.load(dynlib("MGLMM_Beta"))5.4 Data set
X <- model.matrix(~ AGE + BMI + SEX + IPAQ, data3)
data <- list(Y1 = data3$ARMS, Y2 = data3$LEGS, Y3 = data3$TRUNK,
Y4 = data3$ANDROID, Y5 = data3$GYNOID, X = X)5.5 Initial values from betareg package
beta1 = as.numeric(coef(fit1)[-7])
beta2 = as.numeric(coef(fit2)[-7])
beta3 = as.numeric(coef(fit3)[-7])
beta4 = as.numeric(coef(fit4)[-7])
beta5 = as.numeric(coef(fit5)[-7])
# Precision
phi1 <- log(as.numeric(coef(fit1)[7]))
phi2 <- log(as.numeric(coef(fit2)[7]))
phi3 <- log(as.numeric(coef(fit3)[7]))
phi4 <- log(as.numeric(coef(fit4)[7]))
phi5 <- log(as.numeric(coef(fit5)[7]))
parameters <- list(beta1 = beta1, beta2 = beta2, beta3 = beta3,
beta4 = beta4, beta5 = beta5,
phi = c(phi1, phi2, phi3, phi4, phi5),
U = matrix(0, ncol = 5, nrow = 298),
rho = rep(0,10), sigma = rep(0.11, 5))5.6 Log-likelihood function
MGLMM_Beta <- MakeADFun(data, parameters, DLL = "MGLMM_Beta",
random = "U",
hessian = TRUE, silent = T)5.7 Optimizing via nlminb()
fit_mglmm <- nlminb(start = MGLMM_Beta$par, objective = MGLMM_Beta$fn,
gradient = MGLMM_Beta$gr,
control = list(eval.max = 1000, iter.max = 1000,
abs.tol = 1e-04, rel.tol = 1e-04))## Warning in nlminb(start = MGLMM_Beta$par, objective = MGLMM_Beta$fn, gradient =
## MGLMM_Beta$gr, : NA/NaN function evaluationfit_mglmm## $par
## beta1 beta1 beta1 beta1 beta1
## -2.9537491927 0.0046383771 0.0918552429 -0.9389207160 -0.1166535159
## beta1 beta2 beta2 beta2 beta2
## -0.2619628418 -1.9672071507 0.0001473365 0.0692177283 -0.8835971770
## beta2 beta2 beta3 beta3 beta3
## -0.0467723463 -0.1577925196 -3.1611875259 0.0040957943 0.1047017779
## beta3 beta3 beta3 beta4 beta4
## -0.4356444155 -0.1016530631 -0.1865957786 -3.0316100097 0.0049453633
## beta4 beta4 beta4 beta4 beta5
## 0.1047544125 -0.3767117598 -0.1117457458 -0.1891741565 -1.5116539010
## beta5 beta5 beta5 beta5 beta5
## -0.0009244796 0.0623411756 -0.7485710014 -0.0656823897 -0.1526492326
## rho rho rho rho rho
## 1.2885883549 1.2852405657 0.1192039273 9.5178060962 -0.4353073592
## rho rho rho rho rho
## 7.7765089152 6.3719306256 4.3162482303 1.3488597561 0.0594100007
## sigma sigma sigma sigma sigma
## 0.3074049218 0.2888021111 0.2744367528 0.3022942515 0.2424567260
## phi phi phi phi phi
## 8.6044702209 8.1910687573 10.8583572000 6.3380327317 7.4774895770
##
## $objective
## [1] -3155.856
##
## $convergence
## [1] 0
##
## $iterations
## [1] 317
##
## $evaluations
## function gradient
## 470 318
##
## $message
## [1] "relative convergence (4)"5.8 Getting parameter estimates
report <- sdreport(MGLMM_Beta)5.9 Summary for the fixed effects
summary(report, "fixed", p.value = TRUE)## Estimate Std. Error z value Pr(>|z^2|)
## beta1 -2.9537491927 0.1514013061 -19.5094036 9.134143e-85
## beta1 0.0046383771 0.0011655556 3.9795417 6.904825e-05
## beta1 0.0918552429 0.0067704794 13.5670219 6.282308e-42
## beta1 -0.9389207160 0.0373109021 -25.1647820 9.737430e-140
## beta1 -0.1166535159 0.0545352174 -2.1390492 3.243168e-02
## beta1 -0.2619628418 0.0523640232 -5.0027256 5.652536e-07
## beta2 -1.9672071507 0.1426065524 -13.7946477 2.744905e-43
## beta2 0.0001473365 0.0010973312 0.1342680 8.931906e-01
## beta2 0.0692177283 0.0063713807 10.8638507 1.713673e-27
## beta2 -0.8835971770 0.0351284309 -25.1533346 1.299332e-139
## beta2 -0.0467723463 0.0513979967 -0.9100033 3.628208e-01
## beta2 -0.1577925196 0.0493346757 -3.1984100 1.381877e-03
## beta3 -3.1611875259 0.1344413063 -23.5135139 2.966984e-122
## beta3 0.0040957943 0.0010348929 3.9576987 7.567534e-05
## beta3 0.1047017779 0.0060076606 17.4280448 5.054482e-68
## beta3 -0.4356444155 0.0331198398 -13.1535786 1.623071e-39
## beta3 -0.1016530631 0.0484720614 -2.0971475 3.598052e-02
## beta3 -0.1865957786 0.0465310580 -4.0101340 6.068429e-05
## beta4 -3.0316100097 0.1545654598 -19.6137611 1.179764e-85
## beta4 0.0049453633 0.0011864367 4.1682487 3.069489e-05
## beta4 0.1047544125 0.0069007307 15.1801913 4.783848e-52
## beta4 -0.3767117598 0.0380200344 -9.9082435 3.833258e-23
## beta4 -0.1117457458 0.0555184668 -2.0127671 4.413913e-02
## beta4 -0.1891741565 0.0533191199 -3.5479610 3.882257e-04
## beta5 -1.5116539010 0.1214004619 -12.4517970 1.367316e-35
## beta5 -0.0009244796 0.0009328939 -0.9909804 3.216951e-01
## beta5 0.0623411756 0.0054209359 11.5000761 1.317991e-30
## beta5 -0.7485710014 0.0298810308 -25.0517128 1.672408e-138
## beta5 -0.0656823897 0.0436734978 -1.5039416 1.325964e-01
## beta5 -0.1526492326 0.0419176217 -3.6416482 2.708981e-04
## rho 1.2885883549 0.1287308027 10.0099458 1.378281e-23
## rho 1.2852405657 0.1124908383 11.4252910 3.126066e-30
## rho 0.1192039273 0.0877568424 1.3583434 1.743547e-01
## rho 9.5178060962 5.5194688501 1.7244062 8.463455e-02
## rho -0.4353073592 0.8111829307 -0.5366328 5.915213e-01
## rho 7.7765089152 4.4158232060 1.7610553 7.822904e-02
## rho 6.3719306256 3.5652211129 1.7872470 7.389758e-02
## rho 4.3162482303 2.4667196220 1.7497928 8.015407e-02
## rho 1.3488597561 0.7587222054 1.7778045 7.543597e-02
## rho 0.0594100007 0.5757128768 0.1031938 9.178092e-01
## sigma 0.3074049218 0.0141011794 21.7999440 2.323008e-105
## sigma 0.2888021111 0.0128919670 22.4017105 3.787544e-111
## sigma 0.2744367528 0.0113103508 24.2642123 4.681304e-130
## sigma 0.3022942515 0.0137363555 22.0068745 2.474754e-107
## sigma 0.2424567260 0.0110253496 21.9908424 3.523851e-107
## phi 8.6044702209 3.3615126898 2.5597018 1.047620e-02
## phi 8.1910687573 1.7501704994 4.6801547 2.866585e-06
## phi 10.8583572000 3.9578992824 2.7434648 6.079457e-03
## phi 6.3380327317 0.1321140068 47.9739649 0.000000e+00
## phi 7.4774895770 0.6311663196 11.8470985 2.227711e-325.10 Correlation matrix
Beta_report <- MGLMM_Beta$report()
COR_Beta <- Beta_report$Cor
colnames(COR_Beta) <- rownames(COR_Beta) <- c("Arms", "Legs", "Trunk",
"Android", "Gynoid")
COR_Beta## Arms Legs Trunk Android Gynoid
## Arms 1.0000000 0.7900159 0.7871365 0.7713563 0.8088809
## Legs 0.7900159 1.0000000 0.6666091 0.5877548 0.9749532
## Trunk 0.7871365 0.6666091 1.0000000 0.9905704 0.7815699
## Android 0.7713563 0.5877548 0.9905704 1.0000000 0.7131318
## Gynoid 0.8088809 0.9749532 0.7815699 0.7131318 1.00000005.11 Plot correlation matrix
ggcorrplot::ggcorrplot(COR_Beta, hc.order = F, type = "lower",
lab = TRUE, show.legend = F, lab_size = 6, tl.cex = 16,
ggtheme = ggplot2::theme_light(),
colors = c("gray60", "gray60", "gray60"))
References
Ferrari, Silvia, and Francisco Cribari-Neto. 2004. “Beta Regression for Modelling Rates and Proportions.” Journal of Applied Statistics 31 (7): 799–815.
Kristensen, Kasper, Anders Nielsen, Casper W Berg, Hans Skaug, and Bradley M Bell. 2016. “TMB: Automatic Differentiation and Laplace Approximation.” Journal of Statistical Software 70 (5).
Petterle, Ricardo Rasmussen, Henrique Aparecido Laureano, Guilherme Parreira da Silva, and Wagner Hugo Bonat. 2021. “Multivariate Generalized Linear Mixed Models for Continuous Bounded Outcomes: Analyzing the Body Fat Percentage Data.” Statistical Methods in Medical Research (Accepted): 1–26.
Petterle, Ricardo R, Wagner H Bonat, Cassius T Scarpin, Thaı́sa Jonasson, and Victória ZC Borba. 2020. “Multivariate Quasi-Beta Regression Models for Continuous Bounded Data.” The International Journal of Biostatistics 1 (ahead-of-print).