Bayesian Federated Inference
BFI.Rdbfi function can be used (on the central server) to combine inference results from separate datasets (without combining the data) to approximate what would have been inferred had the datasets been merged. This function can handle linear, logistic and survival regression models.
Usage
bfi(theta_hats = NULL,
A_hats,
Lambda,
family = c("gaussian", "binomial", "survival"),
stratified = FALSE,
strat_par = NULL,
center_spec = NULL,
basehaz = c("weibul", "exp", "gomp", "poly", "pwexp"),
theta_A_polys = NULL,
q_ls,
center_zero_sample = FALSE,
which_cent_zeros,
zero_sample_covs,
refer_cats,
zero_cats,
lev_no_ref_zeros)Arguments
- theta_hats
a list of \(L\) vectors of the maximum a posteriori (MAP) estimates of the model parameters in the \(L\) centers. These vectors must have equal dimensions. See ‘Details’.
- A_hats
a list of \(L\) minus curvature matrices for \(L\) centers. These matrices must have equal dimensions. See ‘Details’.
- family
a character string representing the family name used for the local centers. Can be abbreviated.
- Lambda
a list of \(L+1\) matrices. The \(k^{\th}\) matrix is the chosen inverse variance-covariance matrix of the Gaussian distribution that is used as prior distribution in center \(k\), where \(k=1,2,\ldots,L\). The last matrix is the chosen variance-covariance matrix for the Gaussian prior of the (fictive) combined data set. If
stratified = FALSE, all \(L+1\) matrices must have equal dimensions. While, ifstratified = TRUE, the first \(L\) matrices must have equal dimensions and the last matrix should have a different (greater) dimention than the others. See ‘Details’.- stratified
logical flag for performing the stratified analysis. If
stratified = TRUE, the parameter(s) selected in thestrat_parargument are allowed to be different across centers, except when the argumentcenter_specis notNULL. Default isstratified = FALSE. See ‘Details’ and ‘Examples’.- strat_par
a one- or two-element integer vector for indicating the stratification parameter(s). The values \(1\) and/or \(2\) are/is used to indicate that the ``intercept'' and/or ``sigma2'' are allowed to vary, respectively. This argument is used only when
stratified = TRUEandcenter_spec = NULL. Default isstrat_par = NULL, but ifstratified = TRUE,strat_parcan not beNULLunless there is a center specific variable. For thebinomialfamily the length of the vector should be at most one which refers to ``intercept'', and the value of this element should be \(1\). Forgaussianfamily this vector can be \(1\) for indicating the ``intercept'' only, \(2\) for indicating the ``sigma2'' only, and c(\(1\), \(2\)) for both ``intercept'' and ``sigma2''. See ‘Details’ and ‘Examples’.- center_spec
a vector of \(L\) elements for representing the center specific variable. This argument is used only when
stratified = TRUEandstrat_par = NULL. Each element represents a specific feature of the corresponding center. There must be only one specific value or attribute for each center. This vector could be a numeric, characteristic or factor vector. Note that, the order of the centers in the vectorcenter_specmust be the same as in the list of the argumenttheta_hats. The used data type in the argumentcenter_specmust be categorical. Default iscenter_spec = NULL. See also ‘Details’ and ‘Examples’.- basehaz
a character string representing one of the available baseline hazard functions;
exponential("exp"),Weibull("weibul", the default),Gompertz("gomp"),exponentiated polynomial("poly"), andpiecewise exponential("pwexp"). It is only used whenfamily = "survival". Can be abbreviated.- theta_A_polys
a list with \(L\) elements so that each element is the array
theta_A_ploy(the output of theMAP.estimationfunction,MAP.estimation()$theta_A_ploy) for the correcponding center. This argument,theta_A_polys, is only used iffamily = "survival"andbasehaz = "poly". See ‘Details’.- q_ls
a vector with \(L\) elements in which each element is the order (minus 1) of the exponentiated polynomial baseline hazard function for the corresponding center, i.e., each element is the value of
q_l(the output of theMAP.estimationfunction,MAP.estimation()$q_l). This argument,q_ls, is only used iffamily = "survival"andbasehaz = "poly". It can also be a scalar which represents the maximum value of theq_l's across the centers.- center_zero_sample
logical flag indicating whether the center has a categorical covariate with no observations/individuals in one of the categories. Default is
center_zero_sample = FALSE. For more detailes see ‘References’.- which_cent_zeros
an integer vector representing the center(s) which has one categorical covariate with no individuals in one of the categories. It is used if
center_zero_sample = TRUE.- zero_sample_covs
a vector in which each element is a character string representing the categorical covariate that has no samples/observations in one of its categories for the corresponding center. Each element of the vector can be obtained from the output of the
MAP.estimationfunction for the corresponding center,MAP.estimation()$zero_sample_cov. It is used whencenter_zero_sample = TRUE.- refer_cats
a vector in which each element is a character string representing the reference category for the corresponding center. Each element of the vector can be obtained from the output of the
MAP.estimationfunction for the corresponding center,MAP.estimation()$refer_cat. This vector is used whencenter_zero_sample = TRUE.- zero_cats
a vector in which each element is a character string representing the category with no samples/observations for the corresponding center. Each element of the vector can be obtained from the output of the
MAP.estimationfunction for the corresponding center, i.e.,MAP.estimation()$zero_cat. It is used whencenter_zero_sample = TRUE.- lev_no_ref_zeros
a list in which the number of elements equals the length of the
which_cent_zerosargument. Each element of the list is a vector containing the names of the levels of the categorical covariate that has no samples/observations in one of its categories for the corresponding center. However, the name of the category with no samples and the name of the reference category are excluded from this vector. Each element of the list can be obtained from the output of theMAP.estimationfunction, i.e.,MAP.estimation()$lev_no_ref_zero. This argument is used ifcenter_zero_sample = TRUE.
Value
bfi returns a list containing the following components:
- theta_hat
the vector of estimates obtained by combining the inference results from the \(L\) centers with the
'BFI'methodology. If an intercept was fitted in every center andstratified = FALSE, there is only one general ``intercept'' in this vector, while ifstratified = TRUEandstrat_par = 1, there are \(L\) different intercepts in the model, for each center one;- A_hat
minus the curvature (or Hessian) matrix obtained by the
'BFI'method for the combined model. Ifstratified = TRUE, the dimension of the matrix is always greater than whenstratified = FALSE;- sd
the vector of standard deviation of the estimates in
theta_hatobtained from the matrix inA_hat, i.e., the vector equalssqrt(diag(solve(A_hat)))which equals the square root of the elements at the diagonal of the inverse of theA_hatmatrix.- family
the
familyobject used;- basehaz
the baseline hazard function used;
- stratified
whether a stratified analysis was done or not.;
Details
bfi function implements the BFI approach described in the papers Jonker et. al. (2024a), Pazira et. al. (2024) and Jonker et. al. (2024b) given in the references.
The inference results gathered from different (\(L\)) centers are combined, and the BFI estimates of the model parameters and curvature matrix evaluated at that point are returned.
The inference result from each center must be obtained using the MAP.estimation function separately, and then all of these results (coming from different centers) should be compiled into a list to be used as an input of bfi().
The models in the different centers should be defined in exactly the same way; among others, exactly the same covariates should be included in the models. The parameter vectors should be defined exactly the same, so that the \(L\) vectors and matrices in the input lists theta_hat's and A_hat's are defined in the same way (e.g., the covariates need to be included in the models in the same order).
Note that the order of the elements in the lists theta_hats, A_hats and Lambda, must be the same with respect to the centers, so that in every list the element at the \(\ell^{\th}\) position is from the center \(\ell\). This should also be the case for the vector center_spec.
If for the locations intercept = FALSE, the stratified analysis is not possible anymore for the binomial family.
If stratified = FALSE, both strat_par and center_spec must be NULL (the defaults), while if stratified = TRUE only one of the two must be NULL.
If stratified = FALSE and all the \(L+1\) matrices in Lambda are equal, it is sufficient to give a (list of) one matrix only.
In both cases of the stratified argument (TRUE or FALSE), if only the first \(L\) matrices are equal, the argument Lambda can be a list of two matrices, so that the fist matrix represents the chosen variance-covariance matrix for local centers and the second one is the chosen matrix for the combined data set.
The last matrix of the list in the argument Lambda can be built by the function inv.prior.cov().
If the data type used in the argument center_spec is continuous, one can use stratified = TRUE and center_spec = NULL, and set strat_par not to NULL (i.e., to \(1\), \(2\) or both \((1, 2)\)). Indeed, in this case, the stratification parameter(s) given in the argument strat_par are assumed to be different across the centers.
When family = 'survival' and basehaz = 'poly', the arguments theta_hats and A_hats should not be provided. Instead, the theta_A_polys and q_ls arguments should be defined using the local information, specifically MAP.estimation()$theta_A_poly and MAP.estimation()$q_l, respectively. See the last example in ‘Examples’.
References
Jonker M.A., Pazira H. and Coolen A.C.C. (2024a). Bayesian federated inference for estimating statistical models based on non-shared multicenter data sets, Statistics in Medicine, 43(12): 2421-2438. <https://doi.org/10.1002/sim.10072>
Pazira H., Massa E., Weijers J.A.M., Coolen A.C.C. and Jonker M.A. (2024). Bayesian Federated Inference for Survival Models, arXiv. <https://arxiv.org/abs/2404.17464>
Jonker M.A., Pazira H. and Coolen A.C.C. (2024b). Bayesian Federated Inference for regression models with heterogeneous multi-center populations, arXiv. <https://arxiv.org/abs/2402.02898>
Author
Hassan Pazira and Marianne Jonker
Maintainer: Hassan Pazira hassan.pazira@radboudumc.nl
Examples
#################################################
## Example 1: y ~ Binomial (L = 2 centers) ##
#################################################
# Setting a seed for reproducibility
set.seed(112358)
#------------------------------------#
# Data Simulation for Local Center 1 #
#------------------------------------#
n1 <- 30 # sample size of center 1
X1 <- data.frame(x1=rnorm(n1), # continuous variable
x2=sample(0:2, n1, replace=TRUE)) # categorical variable
# make dummy variables
X1x2_1 <- ifelse(X1$x2 == '1', 1, 0)
X1x2_2 <- ifelse(X1$x2 == '2', 1, 0)
X1$x2 <- as.factor(X1$x2)
# regression coefficients
beta <- 1:4 # beta[1] is the intercept
# linear predictor:
eta1 <- beta[1] + X1$x1 * beta[2] + X1x2_1 * beta[3] + X1x2_2 * beta[4]
# inverse of the link function ( g^{-1}(\eta) = \mu ):
mu1 <- binomial()$linkinv(eta1)
y1 <- rbinom(n1, 1, mu1)
#------------------------------------#
# Data Simulation for Local Center 2 #
#------------------------------------#
n2 <- 50 # sample size of center 2
X2 <- data.frame(x1=rnorm(n2), # continuous variable
x2=sample(0:2, n2, replace=TRUE)) # categorical variable
# make dummy variables:
X2x2_1 <- ifelse(X2$x2 == '1', 1, 0)
X2x2_2 <- ifelse(X2$x2 == '2', 1, 0)
X2$x2 <- as.factor(X2$x2)
# linear predictor:
eta2 <- beta[1] + X2$x1 * beta[2] + X2x2_1 * beta[3] + X2x2_2 * beta[4]
# inverse of the link function:
mu2 <- binomial()$linkinv(eta2)
y2 <- rbinom(n2, 1, mu2)
#---------------------------#
# MAP Estimates at Center 1 #
#---------------------------#
# Assume the same inverse covariance matrix (Lambda) for both centers:
Lambda <- inv.prior.cov(X1, lambda = 0.01, family = 'binomial')
fit1 <- MAP.estimation(y1, X1, family = 'binomial', Lambda)
theta_hat1 <- fit1$theta_hat # intercept and coefficient estimates
A_hat1 <- fit1$A_hat # minus the curvature matrix
#---------------------------#
# MAP Estimates at Center 2 #
#---------------------------#
fit2 <- MAP.estimation(y2, X2, family='binomial', Lambda)
theta_hat2 <- fit2$theta_hat
A_hat2 <- fit2$A_hat
#-----------------------#
# BFI at Central Center #
#-----------------------#
A_hats <- list(A_hat1, A_hat2)
theta_hats <- list(theta_hat1, theta_hat2)
bfi <- bfi(theta_hats, A_hats, Lambda, family='binomial')
class(bfi)
#> [1] "bfi"
summary(bfi, cur_mat=TRUE)
#>
#> Summary of the BFI model:
#>
#> Family: ‘binomial’
#> Link: ‘Logit’
#>
#> Coefficients:
#>
#> Estimate Std.Dev CI 2.5% CI 97.5%
#> (Intercept) 1.4479 0.8145 -0.1485 3.0444
#> x1 2.4951 0.8978 0.7355 4.2547
#> x21 4.1023 1.6217 0.9239 7.2807
#> x22 2.2199 1.2637 -0.2570 4.6968
#>
#> Dispersion parameter (sigma2): 1
#>
#> Minus the Curvature Matrix:
#>
#> (Intercept) x1 x21 x22
#> (Intercept) 3.8045 -2.8450 0.7721 0.8694
#> x1 -2.8450 4.0346 -1.2184 -0.5469
#> x21 0.7721 -1.2184 0.7821 0.0000
#> x22 0.8694 -0.5469 0.0000 0.8794
###---------------------###
### Stratified Analysis ###
###---------------------###
# By running the following line an error appears because
# when stratified = TRUE, both 'strat_par' and 'center_spec' can not be NULL:
Just4check1 <- try(bfi(theta_hats, A_hats, Lambda, family = 'binomial',
stratified = TRUE), TRUE)
class(Just4check1) # By default, both 'strat_par' and 'center_spec' are NULL!
#> [1] "try-error"
# By running the following line an error appears because when stratified = TRUE,
# last matrix in 'Lambda' should not have the same dim. as the other local matrices:
Just4check2 <- try(bfi(theta_hats, A_hats, Lambda, stratified = TRUE,
strat_par = 1), TRUE)
class(Just4check2) # All matices in Lambda have the same dimension!
#> [1] "try-error"
# Stratified analysis when 'intercept' varies across two centers:
newLam <- inv.prior.cov(X1, lambda=c(0.1, 0.3), family = 'binomial',
stratified = TRUE, strat_par = 1)
bfi <- bfi(theta_hats, A_hats, list(Lambda, newLam), family = 'binomial',
stratified=TRUE, strat_par=1)
summary(bfi, cur_mat=TRUE)
#>
#> Summary of the BFI model:
#>
#> Family: ‘binomial’
#> Link: ‘Logit’
#>
#> Coefficients:
#>
#> Estimate Std.Dev CI 2.5% CI 97.5%
#> (Intercept)_loc1 2.3167 2.5315 -2.6448 7.2783
#> (Intercept)_loc2 1.1865 0.7586 -0.3003 2.6733
#> x1 1.4500 0.7562 -0.0321 2.9322
#> x21 1.9848 1.1901 -0.3479 4.3174
#> x22 1.3621 1.0287 -0.6542 3.3784
#>
#> Dispersion parameter (sigma2): 1
#>
#> Minus the Curvature Matrix:
#>
#> (Intercept)_loc1 (Intercept)_loc2 x1 x21 x22
#> (Intercept)_loc1 0.1564 0.0000 0.0028 0.0082 0.0134
#> (Intercept)_loc2 0.0000 3.8380 -2.8478 0.7640 0.8561
#> x1 0.0028 -2.8478 4.3246 -1.2184 -0.5469
#> x21 0.0082 0.7640 -1.2184 1.0721 0.0000
#> x22 0.0134 0.8561 -0.5469 0.0000 1.1694
#################################################
## Example 2: y ~ Gaussian (L = 3 centers) ##
#################################################
# Setting a seed for reproducibility
set.seed(112358)
p <- 3 # number of coefficients without 'intercept'
theta <- c(1, rep(2, p), 1.5) # reg. coef.s (theta[1] is 'intercept') & 'sigma2' = 1.5
#------------------------------------#
# Data Simulation for Local Center 1 #
#------------------------------------#
n1 <- 30 # sample size of center 1
X1 <- data.frame(matrix(rnorm(n1 * p), n1, p)) # continuous variables
# linear predictor:
eta1 <- theta[1] + as.matrix(X1)
# inverse of the link function ( g^{-1}(\eta) = \mu ):
mu1 <- gaussian()$linkinv(eta1)
y1 <- rnorm(n1, mu1, sd = sqrt(theta[5]))
#------------------------------------#
# Data Simulation for Local Center 2 #
#------------------------------------#
n2 <- 40 # sample size of center 2
X2 <- data.frame(matrix(rnorm(n2 * p), n2, p)) # continuous variables
# linear predictor:
eta2 <- theta[1] + as.matrix(X2)
# inverse of the link function:
mu2 <- gaussian()$linkinv(eta2)
y2 <- rnorm(n2, mu2, sd = sqrt(theta[5]))
#------------------------------------#
# Data Simulation for Local Center 3 #
#------------------------------------#
n3 <- 50 # sample size of center 3
X3 <- data.frame(matrix(rnorm(n3 * p), n3, p)) # continuous variables
# linear predictor:
eta3 <- theta[1] + as.matrix(X3)
# inverse of the link function:
mu3 <- gaussian()$linkinv(eta3)
y3 <- rnorm(n3, mu3, sd = sqrt(theta[5]))
#---------------------------#
# Inverse Covariance Matrix #
#---------------------------#
# Creating the inverse covariance matrix for the Gaussian prior distribution:
Lambda <- inv.prior.cov(X1, lambda = 0.05, family='gaussian') # the same for both centers
#---------------------------#
# MAP Estimates at Center 1 #
#---------------------------#
fit1 <- MAP.estimation(y1, X1, family = 'gaussian', Lambda)
theta_hat1 <- fit1$theta_hat # intercept and coefficient estimates
A_hat1 <- fit1$A_hat # minus the curvature matrix
#---------------------------#
# MAP Estimates at Center 2 #
#---------------------------#
fit2 <- MAP.estimation(y2, X2, family = 'gaussian', Lambda)
theta_hat2 <- fit2$theta_hat
A_hat2 <- fit2$A_hat
#---------------------------#
# MAP Estimates at Center 3 #
#---------------------------#
fit3 <- MAP.estimation(y3, X3, family = 'gaussian', Lambda)
theta_hat3 <- fit3$theta_hat
A_hat3 <- fit3$A_hat
#-----------------------#
# BFI at Central Center #
#-----------------------#
A_hats <- list(A_hat1, A_hat2, A_hat3)
theta_hats <- list(theta_hat1, theta_hat2, theta_hat3)
bfi <- bfi(theta_hats, A_hats, Lambda, family = 'gaussian')
summary(bfi, cur_mat=TRUE)
#>
#> Summary of the BFI model:
#>
#> Family: ‘gaussian’
#> Link: ‘identity’
#>
#> Coefficients:
#>
#> Estimate Std.Dev CI 2.5% CI 97.5%
#> (Intercept) 0.8674 0.1008 0.6698 1.0649
#> X1 0.9386 0.1010 0.7407 1.1366
#> X2 -0.0077 0.0969 -0.1977 0.1823
#> X3 0.0258 0.1001 -0.1704 0.2220
#>
#> Dispersion parameter (sigma2): 1.177
#>
#> Minus the Curvature Matrix:
#>
#> (Intercept) X1 X2 X3 sigma2
#> (Intercept) 103.1207 -1.6070 17.4661 -13.7946 -0.2544
#> X1 -1.6070 98.6611 7.0539 -2.9114 -0.2818
#> X2 17.4661 7.0539 109.9461 -1.7567 0.0086
#> X3 -13.7946 -2.9114 -1.7567 101.7592 1.5926
#> sigma2 -0.2544 -0.2818 0.0086 1.5926 240.6152
###---------------------###
### Stratified Analysis ###
###---------------------###
# Stratified analysis when 'intercept' varies across two centers:
newLam1 <- inv.prior.cov(X1, lambda = c(0.1,0.3), family = 'gaussian',
stratified = TRUE, strat_par = 1, L=3)
# 'newLam1' is used the prior for combined data and
# 'Lambda' is used the prior for locals
list_newLam1 <- list(Lambda, newLam1)
bfi1 <- bfi(theta_hats, A_hats, list_newLam1, family = 'gaussian',
stratified = TRUE, strat_par = 1)
summary(bfi1, cur_mat = TRUE)
#>
#> Summary of the BFI model:
#>
#> Family: ‘gaussian’
#> Link: ‘identity’
#>
#> Coefficients:
#>
#> Estimate Std.Dev CI 2.5% CI 97.5%
#> (Intercept)_loc1 0.7071 0.2132 0.2892 1.1249
#> (Intercept)_loc2 0.9440 0.1681 0.6145 1.2735
#> (Intercept)_loc3 0.8817 0.1539 0.5801 1.1834
#> X1 0.9496 0.1022 0.7493 1.1498
#> X2 -0.0142 0.0973 -0.2049 0.1764
#> X3 0.0336 0.1016 -0.1655 0.2327
#>
#> Dispersion parameter (sigma2): 1.176
#>
#> Minus the Curvature Matrix:
#>
#> (Intercept)_loc1 (Intercept)_loc2 (Intercept)_loc3 X1
#> (Intercept)_loc1 22.1365 0.0000 0.0000 3.3627
#> (Intercept)_loc2 0.0000 38.6052 0.0000 -7.1754
#> (Intercept)_loc3 0.0000 0.0000 42.6289 2.2058
#> X1 3.3627 -7.1754 2.2058 98.7111
#> X2 1.1877 9.8332 6.4452 7.0539
#> X3 -1.3802 -13.1611 0.7468 -2.9114
#> sigma2 -0.0722 -0.0941 -0.0881 -0.2818
#> X2 X3 sigma2
#> (Intercept)_loc1 1.1877 -1.3802 -0.0722
#> (Intercept)_loc2 9.8332 -13.1611 -0.0941
#> (Intercept)_loc3 6.4452 0.7468 -0.0881
#> X1 7.0539 -2.9114 -0.2818
#> X2 109.9961 -1.7567 0.0086
#> X3 -1.7567 101.8092 1.5926
#> sigma2 0.0086 1.5926 240.8652
# Stratified analysis when 'sigma2' varies across two centers:
newLam2 <- inv.prior.cov(X1, lambda = c(0.1,0.3), family = 'gaussian',
stratified = TRUE, strat_par = 2, L = 3)
# 'newLam2' is used the prior for combined data and 'Lambda' is used the prior for locals
list_newLam2 <- list(Lambda, newLam2)
bfi2 <- bfi(theta_hats, A_hats, list_newLam2, family = 'gaussian',
stratified = TRUE, strat_par=2)
summary(bfi2, cur_mat = TRUE)
#>
#> Summary of the BFI model:
#>
#> Family: ‘gaussian’
#> Link: ‘identity’
#>
#> Coefficients:
#>
#> Estimate Std.Dev CI 2.5% CI 97.5%
#> (Intercept) 0.8672 0.1008 0.6697 1.0647
#> X1 0.9382 0.1010 0.7403 1.1361
#> X2 -0.0076 0.0969 -0.1976 0.1823
#> X3 0.0280 0.1001 -0.1682 0.2241
#> sigma2_loc1 1.3559 0.1285 1.1040 1.6079
#> sigma2_loc2 1.0349 0.1115 0.8164 1.2534
#>
#> Dispersion parameter (sigma2): 1.173
#>
#> Minus the Curvature Matrix:
#>
#> (Intercept) X1 X2 X3 sigma2_loc1 sigma2_loc2
#> (Intercept) 103.1707 -1.6070 17.4661 -13.7946 -0.0722 -0.0941
#> X1 -1.6070 98.7111 7.0539 -2.9114 -0.0915 -0.0983
#> X2 17.4661 7.0539 109.9961 -1.7567 0.0081 0.0020
#> X3 -13.7946 -2.9114 -1.7567 101.8092 -0.0110 1.6065
#> sigma2_loc1 -0.0722 -0.0915 0.0081 -0.0110 60.5220 0.0000
#> sigma2_loc2 -0.0941 -0.0983 0.0020 1.6065 0.0000 80.4577
#> sigma2_loc3 -0.0881 -0.0920 -0.0015 -0.0030 0.0000 0.0000
#> sigma2_loc3
#> (Intercept) -0.0881
#> X1 -0.0920
#> X2 -0.0015
#> X3 -0.0030
#> sigma2_loc1 0.0000
#> sigma2_loc2 0.0000
#> sigma2_loc3 100.4855
# Stratified analysis when 'intercept' and 'sigma2' vary across 2 centers:
newLam3 <- inv.prior.cov(X1, lambda = c(0.1,0.2,0.3), family = 'gaussian',
stratified = TRUE, strat_par = c(1, 2), L = 3)
# 'newLam3' is used the prior for combined data and 'Lambda' is used the prior for locals
list_newLam3 <- list(Lambda, newLam3)
bfi3 <- bfi(theta_hats, A_hats, list_newLam3, family = 'gaussian',
stratified = TRUE, strat_par = 1:2)
summary(bfi3, cur_mat = TRUE)
#>
#> Summary of the BFI model:
#>
#> Family: ‘gaussian’
#> Link: ‘identity’
#>
#> Coefficients:
#>
#> Estimate Std.Dev CI 2.5% CI 97.5%
#> (Intercept)_loc1 0.7079 0.2130 0.2904 1.1255
#> (Intercept)_loc2 0.9443 0.1597 0.6313 1.2572
#> (Intercept)_loc3 0.8817 0.1533 0.5812 1.1822
#> X1 0.9487 0.1019 0.7490 1.1483
#> X2 -0.0142 0.0941 -0.1986 0.1703
#> X3 0.0358 0.0992 -0.1586 0.2302
#> sigma2_loc1 1.3558 0.1285 1.1038 1.6077
#> sigma2_loc2 1.0348 0.1115 0.8163 1.2534
#>
#> Dispersion parameter (sigma2): 1.173
#>
#> Minus the Curvature Matrix:
#>
#> (Intercept)_loc1 (Intercept)_loc2 (Intercept)_loc3 X1
#> (Intercept)_loc1 22.1365 0.0000 0.0000 3.3627
#> (Intercept)_loc2 0.0000 38.6052 0.0000 -7.1754
#> (Intercept)_loc3 0.0000 0.0000 42.6289 2.2058
#> X1 3.3627 -7.1754 2.2058 98.8111
#> X2 1.1877 9.8332 6.4452 7.0539
#> X3 -1.3802 -13.1611 0.7468 -2.9114
#> sigma2_loc1 -0.0722 0.0000 0.0000 1.1877
#> sigma2_loc2 0.0000 -0.0941 0.0000 9.8332
#> sigma2_loc3 0.0000 0.0000 -0.0881 6.4452
#> X2 X3 sigma2_loc1 sigma2_loc2 sigma2_loc3
#> (Intercept)_loc1 -1.3802 0.0081 -0.0722 0.0000 0.0000
#> (Intercept)_loc2 -13.1611 0.0020 0.0000 -0.0941 0.0000
#> (Intercept)_loc3 0.7468 -0.0015 0.0000 0.0000 -0.0881
#> X1 7.0539 -2.9114 -0.0915 -0.0983 -0.0920
#> X2 110.0961 -1.7567 0.0081 0.0020 -0.0015
#> X3 -1.7567 101.9092 -0.0110 1.6065 -0.0030
#> sigma2_loc1 -0.0915 -0.0110 60.5220 0.0000 0.0000
#> sigma2_loc2 -0.0983 1.6065 0.0000 80.4577 0.0000
#> sigma2_loc3 -0.0920 -0.0030 0.0000 0.0000 100.4855
###----------------------------###
### Center Specific Covariates ###
###----------------------------###
# Assume the first and third centers have the same center-specific covariate value
# of '3', while this value for the second center is '1', i.e., center_spec = c(3,1,3)
newLam4 <- inv.prior.cov(X1, lambda=c(0.1, 0.2, 0.3), family='gaussian',
stratified=TRUE, center_spec = c(3,1,3), L=3)
# 'newLam4' is used the prior for combined data and 'Lambda' is used the prior for locals
l_newLam4 <- list(Lambda, newLam4)
bfi4 <- bfi(theta_hats, A_hats, l_newLam4, family = 'gaussian',
stratified = TRUE, center_spec = c(3,1,3))
summary(bfi4, cur_mat = TRUE)
#>
#> Summary of the BFI model:
#>
#> Family: ‘gaussian’
#> Link: ‘identity’
#>
#> Coefficients:
#>
#> Estimate Std.Dev CI 2.5% CI 97.5%
#> (Intercept)_1 0.9433 0.1681 0.6138 1.2728
#> (Intercept)_3 0.8234 0.1251 0.5781 1.0686
#> X1 0.9458 0.1020 0.7458 1.1458
#> X2 -0.0117 0.0972 -0.2022 0.1787
#> X3 0.0354 0.1015 -0.1635 0.2343
#>
#> Dispersion parameter (sigma2): 1.176
#>
#> Minus the Curvature Matrix:
#>
#> (Intercept)_1 (Intercept)_3 X1 X2 X3 sigma2
#> (Intercept)_1 38.6052 0.0000 -7.1754 9.8332 -13.1611 -0.0941
#> (Intercept)_3 0.0000 64.6655 5.5685 7.6329 -0.6334 -0.1603
#> X1 -7.1754 5.5685 98.8111 7.0539 -2.9114 -0.2818
#> X2 9.8332 7.6329 7.0539 110.0961 -1.7567 0.0086
#> X3 -13.1611 -0.6334 -2.9114 -1.7567 101.9092 1.5926
#> sigma2 -0.0941 -0.1603 -0.2818 0.0086 1.5926 240.8652
####################################################
## Example 3: Survival family (L = 2 centers) ##
####################################################
# Setting a seed for reproducibility
set.seed(112358)
p <- 4
theta <- c(1:4, 5, 6) # regression coefficients (1:4) & omega's (5:6)
#---------------------------------------------#
# Simulating Survival data for Local Center 1 #
#---------------------------------------------#
n1 <- 50
X1 <- data.frame(matrix(rnorm(n1 * p), n1, p)) # continuous (normal) variables
# Simulating survival data ('time' and 'status') from 'Weibull' with
# a predefined censoring rate of 0.3:
y1 <- surv.simulate(Z = list(X1), beta = theta[1:p], a = theta[5], b = theta[6],
u1 = 0.1, cen_rate = 0.3, gen_data_from = "weibul")$D[[1]][, 1:2]
Lambda <- inv.prior.cov(X1, lambda = c(0.1, 1), family = "survival", basehaz ="poly")
fit1 <- MAP.estimation(y1, X1, family = 'survival', Lambda = Lambda, basehaz = "poly")
theta_hat1 <- fit1$theta_hat # coefficient estimates
A_hat1 <- fit1$A_hat # minus the curvature matrix
summary(fit1, cur_mat=TRUE)
#>
#> Summary of the local model:
#>
#> Formula: Survival(time, status) ~ X1 + X2 + X3 + X4
#> Family: ‘survival’
#> Baseline: ‘poly’
#>
#> Coefficients:
#>
#> Estimate Std.Dev CI 2.5% CI 97.5%
#> X1 0.5096 0.1980 0.1215 0.8978
#> X2 0.6952 0.1906 0.3217 1.0687
#> X3 1.1450 0.1782 0.7957 1.4942
#> X4 1.4295 0.1554 1.1249 1.7341
#> omega_0 -1.5109 0.1967 -1.8963 -1.1254
#> omega_1 2.0833 0.3130 1.4699 2.6967
#> omega_2 0.5643 0.3092 -0.0417 1.1703
#>
#> log Lik Posterior: 1.275
#> Convergence: 0
#>
#> Minus the Curvature Matrix:
#>
#> X1 X2 X3 X4 omega_0 omega_1 omega_2
#> X1 28.9154 0.0644 -0.0075 8.4374 11.7632 4.4057 2.9786
#> X2 0.0644 32.5920 -8.1674 5.8666 -2.2318 -6.2069 -7.9569
#> X3 -0.0075 -8.1674 36.8872 3.2223 3.5844 -4.0484 -5.1086
#> X4 8.4374 5.8666 3.2223 58.7733 15.8110 -6.0420 -11.0501
#> omega_0 11.7632 -2.2318 3.5844 15.8110 48.5109 17.0175 17.9298
#> omega_1 4.4057 -6.2069 -4.0484 -6.0420 17.0175 26.9298 22.5453
#> omega_2 2.9786 -7.9569 -5.1086 -11.0501 17.9298 22.5453 31.5303
#---------------------------------------------#
# Simulating Survival data for Local Center 2 #
#---------------------------------------------#
n2 <- 50
X2 <- data.frame(matrix(rnorm(n2 * p), n2, p)) # continuous (normal) variables
# Survival simulated data from 'Weibull' with a predefined censoring rate of 0.3:
y2 <- surv.simulate(Z = list(X2), beta = beta, a = theta[5], b = theta[6], u1 = 0.1,
cen_rate = 0.3, gen_data_from = "weibul")$D[[1]][, 1:2]
fit2 <- MAP.estimation(y2, X2, family = 'survival', Lambda = Lambda, basehaz = "poly")
theta_hat2 <- fit2$theta_hat
A_hat2 <- fit2$A_hat
summary(fit2, cur_mat=TRUE)
#>
#> Summary of the local model:
#>
#> Formula: Survival(time, status) ~ X1 + X2 + X3 + X4
#> Family: ‘survival’
#> Baseline: ‘poly’
#>
#> Coefficients:
#>
#> Estimate Std.Dev CI 2.5% CI 97.5%
#> X1 0.2311 0.1717 -0.1054 0.5677
#> X2 0.6308 0.1939 0.2509 1.0108
#> X3 1.2914 0.2092 0.8813 1.7015
#> X4 1.4830 0.1763 1.1375 1.8285
#> omega_0 -1.2579 0.1868 -1.6240 -0.8918
#> omega_1 2.5449 0.3092 1.9389 3.1509
#> omega_2 -0.0842 0.2475 -0.5693 0.4009
#>
#> log Lik Posterior: -1.517
#> Convergence: 0
#>
#> Minus the Curvature Matrix:
#>
#> X1 X2 X3 X4 omega_0 omega_1 omega_2
#> X1 36.4073 -1.5410 6.7110 3.8278 -1.8760 -5.6677 -8.6908
#> X2 -1.5410 35.1388 12.6690 -0.2949 17.1580 5.0885 3.5375
#> X3 6.7110 12.6690 34.8678 3.3796 11.0675 -5.1800 -12.5700
#> X4 3.8278 -0.2949 3.3796 41.5954 3.6979 -10.7514 -16.9708
#> omega_0 -1.8760 17.1580 11.0675 3.6979 52.2579 17.8297 20.2039
#> omega_1 -5.6677 5.0885 -5.1800 -10.7514 17.8297 29.2039 30.1189
#> omega_2 -8.6908 3.5375 -12.5700 -16.9708 20.2039 30.1189 50.9909
#-----------------------#
# BFI at Central Center #
#-----------------------#
# When family = 'survival' and basehaz = "poly", only 'theta_A_polys'
# should be defined instead of 'theta_hats' and 'A_hats':
theta_A_hats <- list(fit1$theta_A_poly, fit2$theta_A_poly)
qls <- c(fit1$q_l, fit2$q_l)
bfi <- bfi(Lambda = Lambda, family = 'survival', theta_A_polys = theta_A_hats,
basehaz = "poly", q_ls = qls)
summary(bfi, cur_mat=TRUE)
#>
#> Summary of the BFI model:
#>
#> Family: ‘survival’
#> Baseline: ‘poly’
#>
#> Coefficients:
#>
#> Estimate Std.Dev CI 2.5% CI 97.5%
#> X1 0.3504 0.1269 0.1016 0.5991
#> X2 0.6991 0.1255 0.4530 0.9451
#> X3 1.2232 0.1288 0.9707 1.4756
#> X4 1.4394 0.1151 1.2137 1.6650
#> omega_0 -1.3841 0.1346 -1.6479 -1.1203
#> omega_1 2.4828 0.2253 2.0413 2.9244
#> omega_2 0.0840 0.1950 -0.2981 0.4662
#>
#> Minus the Curvature Matrix:
#>
#> X1 X2 X3 X4 omega_0 omega_1 omega_2
#> X1 65.2227 -1.4767 6.7036 12.2651 9.8872 -1.2620 -5.7122
#> X2 -1.4767 67.6307 4.5016 5.5717 14.9263 -1.1183 -4.4194
#> X3 6.7036 4.5016 71.6550 6.6019 14.6519 -9.2284 -17.6786
#> X4 12.2651 5.5717 6.6019 100.2687 19.5089 -16.7934 -28.0209
#> omega_0 9.8872 14.9263 14.6519 19.5089 99.7688 34.8472 38.1337
#> omega_1 -1.2620 -1.1183 -9.2284 -16.7934 34.8472 55.1337 52.6642
#> omega_2 -5.7122 -4.4194 -17.6786 -28.0209 38.1337 52.6642 81.5212