Maximum A Posteriori estimation

MAP.estimation function is used (in local centers) to compute Maximum A Posterior (MAP) estimators of the parameters for Generalized Linear Models (GLM) and Survival models.

Usage

MAP.estimation(y,
               X,
               family = c("gaussian", "binomial", "survival"),
               Lambda,
               intercept = TRUE,
               basehaz = c("weibul", "exp", "gomp", "poly", "pwexp", "unspecified"),
               treatment = NULL,
               refer_treat,
               coef_bfi = NULL,
               RCT_propens = NULL,
               initial = NULL,
               alpha = 0.1,
               max_order = 2,
               n_intervals = 4,
               min_max_times,
               center_zero_sample = FALSE,
               zero_sample_cov,
               refer_cat,
               zero_cat,
               control = list())

Arguments

y: response vector. If the “binomial” family is used, this argument is a vector with entries 0 (failure) or 1 (success). Alternatively, for this family, the response can be a matrix where the first column is the number of “successes” and the second column is the number of “failures”. For the “survival” family, the response is a matrix where the first column is the survival time, named “time”, and the second column is the censoring indicator, named “status”, with 0 indicating censoring time and 1 indicating event time.
X: design matrix of dimension $n_{\ell} \times p$, where $p$ is the number of covariates or predictors and $n_{\ell}$ is the number of indeviduals in the local center. Note that the order of the covariates must be the same across the centers; otherwise, the output estimates of bfi() will be incorrect.
family: a description of the error distribution. This is a character string naming a family of the model. In the current version of the package, the family of model can be “gaussian” (with identity link function), dQuotebinomial (with logit link function), or “survival”. Can be abbreviated. By default the “gaussian” family is used. In case of a linear regression model, family = "gaussian", there is an extra model parameter for the variance of measurement error. While in the case of survival model, family = "survival", the number of the model parameters depend on the choice of baseline hazard functions, see ‘Details’ for more information.
Lambda: the inverse variance-covariance matrix of the Gaussian distribution that is used as prior distribution for the model parameters. The dimension of the matrix depends on the number of columns of X, type of the covariates (continuous / dichotomous or categorical), family, and whether an intercept is included (if applicable). However, Lambda can be easily created by inv.prior.cov(). See inv.prior.cov for more information.
intercept: logical flag for fitting an intercept. If intercept=TRUE (the default), the intercept is fitted, i.e., it is included in the model, and if intercept = FALSE it is set to zero, i.e., it's not in the model. This argument is not used if family = "survival".
basehaz: a character string representing one of the available baseline hazard functions; exponential (“exp”), Weibull (“weibul”, the default), Gompertz (“gomp”), exponentiated polynomial (“poly”), piecewise constant exponential (“pwexp”), and unspecified baseline hazard (“unspecified”). Can be abbreviated. It is used only when family = "survival". If local sample size is large and the shape of the baseline hazard function is completely unknown, the “exponentiated polynomial” and “piecewise exponential” hazard functions would be preferred above the lower dimensional alternatives. However, if the local samples size is low, one should be careful using the “piecewise exponential” hazard function with many intervals. If basehaz = "unspecified", it means that a (semi-parametric) Cox model is considered, and the parameters (regression coefficients) are estimated using the partial log-likelihood. If treatment is not NULL, then basehaz must be set to "unspecified", as the regression coefficients are estimated using the weighted partial log-likelihood.
treatment: a character string representing the name or place of the binary covariate, respectively. This covariate indicates whether the patient got the new treatment ($z_{\ell i}=1$) or the placebo/standard treatment ($z_{\ell i}=0$). The treatment effect is estimated when this argument is NOT 'NULL'. If it is set to 'NULL' (the default), the treatment effect will not be estimated. In fact, for the first round, this argument should be 'NULL', and in the second round, it should not be 'NULL'. See ‘Details’.
refer_treat: a character string representing the reference category of the treatment variable. The reference category is considered as $z_{\ell i}=0$. This argument is used when treatment is not 'NULL'. Default is refer_treat = levels(X$treatment)[1].
coef_bfi: a vector specifying the BFI estimates of the coefficients received from the central server in the first round. The length of coef_bfi equals the number of regression coefficients, including the intercept if intercept=TRUE, but excluding $\zeta$, which represents the treatment effect, as well as the nuisance parameter $\sigma$ in the gaussian family and any parameters of the baseline hazard ($\boldsymbol{\omega}$) for survival. This argument is used only when the argument treatment is not 'NULL'. If treatment is not 'NULL' but coef_bfi = NULL, then the argument RCT_propens must not be 'NULL', indicating an RCT study. See ‘Details’.
RCT_propens: a vector specifying the propensity scores, which represent the probability of receiving the treatment given the covariates, which are known in the randomized studies (RCTs). For example, in a 1:1 randomized trial, the propensity scores are, by definition, equal to 1/2 (or 0.5), whereas in an unbalanced randomized trial, e.g., a 2:1 trial, the propensity scores are now 2/3 and 1/3 for the two arms, respectively. The length of RCT_propens equals to the number of individuals in the local center denoted as $n_{\ell}$. This argument is used only when the study is a randomized control trial, i.e., the propensity scores are known for this local center. In this case, there is no ‘second round’, and the argument treatment must not be 'NULL', whereas coef_bfi = NULL. Indeed, when 'treatment' is not 'NULL', one of the arguments 'RCT_propens' or 'coef_bfi' could be 'NULL'. See ‘Details’.
initial: a vector specifying initial values for the parameters to be optimized over. The length of initial is equal to the number of model parameters and thus, is equal to the number of rows or columns of Lambda. Since the 'L-BFGS-B' method is used in the algorithm, these values should always be finite. Default is a vector of zeros, except for the survival family with the poly function, where it is a vector with the first $p$ elements as zeros for coefficients ($\boldsymbol{\beta}$) and -0.5 for the remaining parameters ($\boldsymbol{\omega}$). For the gaussian family, the last element of the initial vector could also be considered negative, because the Gaussian prior was applied to $log(\sigma^2)$.
alpha: a significance level used in the chi-squared distribution (with one degree of freedom and 1-$\alpha$ representing the upper quantile) to conduct a likelihood ratio test for obtaining the order of the exponentiated polynomial baseline hazard function. It is only used when family = "survival" and basehaz = "poly". Default is 0.1. See ‘Details’.
max_order: an integer representing the maximum value of q_l, which is the order/degree minus 1 of the exponentiated polynomial baseline hazard function. This argument is only used when family = "survival" and basehaz = "poly". Default is 2.
n_intervals: an integer representing the number of intervals in the piecewise exponential baseline hazard function. This argument is only used when family = "survival" and basehaz = "pwexp". Default is 4.
min_max_times: a scalar representing the minimum of the maximum event times observed in the centers. The value of this argument should be defined by the central server (which has access to the maximum event times of all the centers) and is only used when family = "survival" and basehaz = "pwexp".
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.
zero_sample_cov: either a character string or an integer representing the categorical covariate that has no samples/observations in one of its categories. This covariate should have at least two categories, one of which is the reference. It is used when center_zero_sample = TRUE.
refer_cat: a character string representing the reference category. The category with no observations (the argument zero_cat) cannot be used as the reference in the argument refer_cat. It is used when center_zero_sample = TRUE.
zero_cat: a character string representing the category with no samples/observations. It is used when center_zero_sample = TRUE.
control: a list of control parameters. See ‘Details’.

Value

MAP.estimation returns a list containing the following components:

theta_hat

the vector corresponding to the maximum a posteriori (MAP) estimates of the parameters. For the gaussian family, the last element of this vector is $\sigma^2$;

A_hat

minus the curvature (or Hessian) matrix around the point theta_hat. The dimension of the matrix is the same as the argument Lambda;

sd

the vector of (posterior) standard deviation of the MAP estimates in theta_hat, that is sqrt(diag(solve(A_hat)));

Lambda

the inverse variance-covariance matrix of the Gaussian distribution that is used as prior distribution for the parameters. It's exactly the same as the argument Lambda;

formula

the formula applied;

names

the names of the model parameters;

n

sample size, $n_{\ell}$;

np

the number of coefficients;

q_l

the order/degree minus 1 of the exponentiated polynomial baseline hazard function determined for the current center by the likelihood ratio test. This output argument, q_l, is only shown when family = "survival" and basehaz = "poly", and will be used in the function bfi();

theta_A_poly

an array where the first component is a matrix with columns representing the MAP estimates of the parameters for different q_l's, i.e., q_l, q_l+1, ..., max_order. The other components are minus the curvature matrices for different q_l's, i.e., q_l, q_l+1, ..., max_order. Therefore, the first non-NA curvature matrix is equal to the output argument A_hat. This output argument, theta_A_poly, is only shown if family = "survival" and basehaz = "poly", and will be used in the function bfi();

lev_no_ref_zer

a vector containing the names of the levels of the categorical covariate that has no samples/observations in one of its categories. The name of the category with no samples and the name of the reference category are excluded from this vector. This argument is shown when family = "survival" and basehaz = "poly", and will be used in the function bfi();

treatment

a character string representing the name or place of the binary covariate, respectively. If it is set to 'NULL', the treatment effect will not be estimated;

refer_treat

the reference category of the treatment. It is shown when treatment is not 'NULL', and can be used in the function bfi();

coef_bfi

a vector specifying the BFI estimates of the coefficients received from the central server in the first round. If treatment = NULL, then coef_bfi must also be 'NULL';

RCT_propens

a vector specifying the propensity scores, which represent the probability of receiving the treatment given the covariates, which are known in the randomized studies (RCTs). If treatment = NULL, then RCT_propens must also be 'NULL';

propensity

a vector specifying the propensity scores (the probability a patient gets the treatment given the characteristics measured at baseline) calculated by $Pr(Z_\ell = 1 | X_\ell)$;

for_ATE

a vector used in the central server to calculate the average treatment effect (ATE). For family of binomial and gaussian, its elements are:

first: $m_{\ell 1}$, the number of patients in the treatment group, where $n_{\ell} = m_{\ell 1} + m_{\ell 2}$,
second: $m_{\ell 2}$, the number of patients in the reference group, where $m_{\ell 2} = n_{\ell} - m_{\ell 1}$,
third: $\sum_{i=1}^{n_{\ell}} z_{\ell i} y_{\ell i}$,
fourth: $\sum_{i=1}^{n_{\ell}} (z_{\ell i} y_{\ell i})^{2} $,
fifth: $\sum_{i=1}^{n_{\ell}} z_{\ell i} / e_{\ell i}$,
sixth: $\sum_{i=1}^{n_{\ell}} z_{\ell i} y_{\ell i} / e_{\ell i}$,
seventh: $\sum_{i=1}^{n_{\ell}} (1 - z_{\ell i}) / (1 - e_{\ell i})$,
eighth: $\sum_{i=1}^{n_{\ell}} (1 - z_{\ell i}) y_{\ell i} / (1 - e_{\ell i})$,

but in survival, its elements are:

first: $\hat \zeta_\ell$, the weighted MAP estimates of the treatment effect $ \zeta_\ell$ in center $\ell$,
second: $\hat{\bold{A}}_\ell$, minus the curvature (1$\times$1)-dimensional matrix estimates (variance) around $ \zeta_\ell$ in center $\ell$,
third: standard deviation of $ \zeta_\ell$ in center $\ell$, i.e., , sqrt(solve(for_ATE[2]))
fourth: an integer value used to encode the warnings and the errors related to the algorithm used to maximaze the weighted partial log likelihood;

zero_sample_cov

the categorical covariate that has no samples/observations in one of its categories. It is shown when center_zero_sample = TRUE, and can be used in the function bfi();

refer_cat

the reference category. It is shown when center_zero_sample = TRUE, and can be used in the function bfi();

zero_cat

the category with no samples/observations. It is shown when center_zero_sample = TRUE, and can be used in the function bfi();

value

the value of minus the log-likelihood posterior density evaluated at theta_hat;

family

the family used;

basehaz

the baseline hazard function used;

intercept

logical flag used to fit an intercept if TRUE, or set to zero if FALSE;

convergence

an integer value used to encode the warnings and the errors related to the algorithm used to fit the model. The values returned are:

0: algorithm has converged,
1: maximum number of iterations ('maxit') has been reached,
2: Warning from the 'L-BFGS-B' method. See the message after this value;

control

the list of control parameters used to compute the MAP estimates.

Details

MAP.estimation function finds the Maximum A Posteriori (MAP) estimates of the model parameters by maximizing the log-posterior density with respect to the parameters, i.e., the estimates equal the values for which the log-posterior density is maximal (the posterior mode). In other words, MAP.estimation() optimizes the log-posterior density with respect to the parameter vector to obtain its MAP estimates. In addition to the model parameters (i.e., coefficients ($\boldsymbol{\beta}$) and variance error ($\sigma^2_e$) for gaussian or the parameters of the baseline hazard ($\boldsymbol{\omega}$) for survival), the curvature matrix (Hessian of the log-posterior) is estimated around the mode.

The MAP.estimation function returns an object of class `bfi`. Therefore, summary() can be used for the object returned by MAP.estimation().

For the case where family = "survival" and basehaz = "poly", we assume that in all centers the $q_\ell$'s are equal. However, the order of the estimated polynomials may vary across the centers so that each center can have different number of parameters, say $q_\ell$+1. After obtaining the estimates within the local centers (by using MAP.estimation()) and having all estimates in the central server, we choose the order of the polynomial approximation for the combined data to be the maximum of the orders of the local polynomial functions, i.e., $\max \{q_1, \ldots, q_L \}$, to approximate the global baseline hazard (exponentiated polynomial) function more accurately. This is because the higher-order polynomial approximation can capture more complex features and details in the combined data. Using the higher-order approximation ensures that we account for the higher-order moments and features present in the data while maintaining accuracy. As a result, all potential cases are stored in the theta_A_poly argument to be used in bfi() by the central server. For further information on the survival family, refer to the 'References' section.

The three arguments treatment, coef_bfi, and RCT_propens are related to the estimation of the treatment effect. For observational and non-randomized studies, the treatment effect is estimated in two rounds; In the first round, only $\hat{\boldsymbol{\beta}}_{BFI}$ is estimated (in this case treatment should not be 'NULL', coef_bfi = NULL, and RCT_propens = NULL); In the second round, propensity scores are estimated, and some summary statistics are sent to the central server to estimate the average treatment estimate, ATE (in this case treatment and coef_bfi should not be 'NULL', but RCT_propens = NULL). In contrast, for the randomized control trial (RCT), the treatment effect can be estimated by only one round as the propensity scores are known (in this case treatment and RCT_propens should not be 'NULL', but coef_bfi = NULL). NOTE: the argument coef_bfi should not include estimates of the nuisance parameter $\sigma$ in the gaussian family or any parameters of the baseline hazard ($\boldsymbol{\omega}$) for survival. Check this: identical(names(coef_bfi), X[, -which(colnames(X) == treatment)]) should be TRUE.

To solve unconstrained and bound-constrained optimization problems, the MAP.estimation function utilizes an optimization algorithm called Limited-memory Broyden-Fletcher-Goldfarb-Shanno with Bound Constraints (L-BFGS-B), Byrd et. al. (1995). The L-BFGS-B algorithm is a limited-memory “quasi-Newton” method that iteratively updates the parameter estimates by approximating the inverse Hessian matrix using gradient information from the history of previous iterations. This approach allows the algorithm to approximate the curvature of the posterior distribution and efficiently search for the optimal solution, which makes it computationally efficient for problems with a large number of variables.

By default, the algorithm uses a relative change in the objective function as the convergence criterion. When the change in the objective function between iterations falls below a certain threshold (`factr`) the algorithm is considered to have converged. The convergence can be checked with the argument convergence in the output. See ‘Value’.

In case of convergence issue, it may be necessary to investigate and adjust optimization parameters to facilitate convergence. It can be done using the initial and control arguments. By the argument initial the initial points of the interative optimization algorithm can be changed, and the argument control is a list that can supply any of the following components:

maxit:: is the maximum number of iterations. Default is 150;
factr:: controls the convergence of the 'L-BFGS-B' method. Convergence occurs when the reduction in the objective is within this factor of the machine tolerance. Default for factr is 1e7, which gives a tolerance of about 1e-9. The exact tolerance can be checked by multiplying this value by .Machine$double.eps;
pgtol:: helps to control the convergence of the 'L-BFGS-B' method. It is a tolerance on the projected gradient in the current search direction, i.e., the iteration will stop when the maximum component of the projected gradient is less than or equal to pgtol, where pgtol$\geq 0$. Default is zero, when the check is suppressed;
trace:: is a non-negative integer. If positive, tracing information on the progress of the optimization is produced. Higher values may produce more tracing information: for the method 'L-BFGS-B' there are six levels of tracing. To understand exactly what these do see the source code of optim function in the stats package;
REPORT:: is the frequency of reports for the 'L-BFGS-B' method if 'control$trace' is positive. Default is every 10 iterations;
lmm:: is an integer giving the number of BFGS updates retained in the 'L-BFGS-B' method. Default is 5.

References

Jonker M.A., Pazira H. and Coolen A.C.C. (2024). 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>

Byrd R.H., Lu P., Nocedal J. and Zhu C. (1995). A limited memory algorithm for bound constrained optimization. SIAM Journal on Scientific Computing, 16, 1190-1208. <https://doi.org/10.1137/0916069>

Author

Hassan Pazira and Marianne Jonker
Maintainer: Hassan Pazira hassan.pazira@radboudumc.nl

Examples


###--------------###
### y ~ Gaussian ###
###--------------###

# Setting a seed for reproducibility
set.seed(11235)

# model parameters: coefficients and sigma2 = 1.5
theta <- c(1, 2, 2, 2, 1.5)

#----------------
# Data Simulation
#----------------
n   <- 30   # sample size
p   <- 3    # number of coefficients without intercept
X   <- data.frame(matrix(rnorm(n * p), n, p)) # continuous variables
# linear predictor:
eta <- theta[1] + theta[2] * X$X1 + theta[3] * X$X2 + theta[4] * X$X3
# inverse of the link function ( g^{-1}(\eta) = \mu ):
mu  <- gaussian()$linkinv(eta)
y   <- rnorm(n, mu, sd = sqrt(theta[5]))

# Load the BFI package
library(BFI)

#-----------------------------------------------
# MAP estimations for theta and curvature matrix
#-----------------------------------------------
# MAP estimates with 'intercept'
Lambda <- inv.prior.cov(X, lambda = c(0.1, 1), family = "gaussian")
(fit <- MAP.estimation(y, X, family = "gaussian", Lambda))
#> $theta_hat
#> (Intercept)          X1          X2          X3      sigma2 
#>    1.341258    2.236391    2.071001    2.164002    1.054571 
#> 
#> $A_hat
#>             (Intercept)         X1         X2         X3     sigma2
#> (Intercept)  28.5475747  1.9315730 -7.7385763  6.2804549 -0.2682462
#> X1            1.9315730 20.4369223 -2.2127254  4.5673270 -0.4472059
#> X2           -7.7385763 -2.2127254 49.3206700  3.1359341 -0.4141474
#> X3            6.2804549  4.5673270  3.1359341 24.3657850 -0.4327445
#> sigma2       -0.2682462 -0.4472059 -0.4141474 -0.4327445 64.2183745
#> 
#> $sd
#> (Intercept)          X1          X2          X3      sigma2 
#>   0.1982976   0.2269483   0.1476575   0.2153269   0.1248065 
#> 
#> $Lambda
#>             (Intercept)  X1  X2  X3 sigma2
#> (Intercept)         0.1 0.0 0.0 0.0      0
#> X1                  0.0 0.1 0.0 0.0      0
#> X2                  0.0 0.0 0.1 0.0      0
#> X3                  0.0 0.0 0.0 0.1      0
#> sigma2              0.0 0.0 0.0 0.0      1
#> 
#> $formula
#> [1] y ~ X1 + X2 + X3
#> 
#> $names
#> [1] "(Intercept)" "X1"          "X2"          "X3"          "sigma2"     
#> 
#> $n
#> [1] 30
#> 
#> $np
#> [1] 4
#> 
#> $treatment
#> NULL
#> 
#> $zero_sample_cov
#> NULL
#> 
#> $refer_cat
#> NULL
#> 
#> $zero_cat
#> NULL
#> 
#> $value
#> [1] 35.28046
#> 
#> $family
#> [1] "gaussian"
#> 
#> $basehaz
#> [1] "weibul"      "exp"         "gomp"        "poly"        "pwexp"      
#> [6] "unspecified"
#> 
#> $intercept
#> [1] TRUE
#> 
#> $convergence
#> [1] 0
#> 
#> $control
#> $control$maxit
#> [1] 100
#> 
#> 
#> attr(,"class")
#> [1] "bfi"
class(fit)
#> [1] "bfi"
summary(fit, cur_mat = TRUE)
#> 
#> Summary of the local model:
#> 
#>    Formula: y ~ X1 + X2 + X3 
#>     Family: ‘gaussian’ 
#>       Link: ‘identity’
#> 
#> Coefficients:
#> 
#>             Estimate Std.Dev CI 2.5% CI 97.5%
#> (Intercept)   1.3413  0.1983  0.9526   1.7299
#> X1            2.2364  0.2269  1.7916   2.6812
#> X2            2.0710  0.1477  1.7816   2.3604
#> X3            2.1640  0.2153  1.7420   2.5860
#> 
#> Dispersion parameter (sigma2):  1.055 
#>             log Lik Posterior:  -35.28 
#>                   Convergence:  0 
#> 
#> Minus the Curvature Matrix: 
#> 
#>             (Intercept)      X1      X2      X3  sigma2
#> (Intercept)     28.5476  1.9316 -7.7386  6.2805 -0.2682
#> X1               1.9316 20.4369 -2.2127  4.5673 -0.4472
#> X2              -7.7386 -2.2127 49.3207  3.1359 -0.4141
#> X3               6.2805  4.5673  3.1359 24.3658 -0.4327
#> sigma2          -0.2682 -0.4472 -0.4141 -0.4327 64.2184

# MAP estimates without 'intercept'
Lambda <- inv.prior.cov(X, lambda = c(0.1, 1), family = 'gaussian', intercept = FALSE)
(fit1 <- MAP.estimation(y, X, family = 'gaussian', Lambda, intercept = FALSE))
#> $theta_hat
#>       X1       X2       X3   sigma2 
#> 2.241637 1.832740 2.525235 2.496099 
#> 
#> $A_hat
#>                X1         X2         X3     sigma2
#> X1      8.6921038 -0.9348497  1.9296405 -0.4483243
#> X2     -0.9348497 20.8951379  1.3248942 -0.3665339
#> X3      1.9296405  1.3248942 10.3520007 -0.5050368
#> sigma2 -0.4483243 -0.3665339 -0.5050368 69.9843284
#> 
#> $sd
#>        X1        X2        X3    sigma2 
#> 0.3478802 0.2205610 0.3192969 0.1195753 
#> 
#> $Lambda
#>         X1  X2  X3 sigma2
#> X1     0.1 0.0 0.0      0
#> X2     0.0 0.1 0.0      0
#> X3     0.0 0.0 0.1      0
#> sigma2 0.0 0.0 0.0      1
#> 
#> $formula
#> [1] y ~ X1 + X2 + X3
#> 
#> $names
#> [1] "X1"     "X2"     "X3"     "sigma2"
#> 
#> $n
#> [1] 30
#> 
#> $np
#> [1] 3
#> 
#> $treatment
#> NULL
#> 
#> $zero_sample_cov
#> NULL
#> 
#> $refer_cat
#> NULL
#> 
#> $zero_cat
#> NULL
#> 
#> $value
#> [1] 63.9101
#> 
#> $family
#> [1] "gaussian"
#> 
#> $basehaz
#> [1] "weibul"      "exp"         "gomp"        "poly"        "pwexp"      
#> [6] "unspecified"
#> 
#> $intercept
#> [1] FALSE
#> 
#> $convergence
#> [1] 0
#> 
#> $control
#> $control$maxit
#> [1] 100
#> 
#> 
#> attr(,"class")
#> [1] "bfi"
summary(fit1, cur_mat = TRUE)
#> 
#> Summary of the local model:
#> 
#>    Formula: y ~ X1 + X2 + X3 
#>     Family: ‘gaussian’ 
#>       Link: ‘identity’
#> 
#> Coefficients:
#> 
#>    Estimate Std.Dev CI 2.5% CI 97.5%
#> X1   2.2416  0.3479  1.5598   2.9235
#> X2   1.8327  0.2206  1.4004   2.2650
#> X3   2.5252  0.3193  1.8994   3.1510
#> 
#> Dispersion parameter (sigma2):  2.496 
#>             log Lik Posterior:  -63.91 
#>                   Convergence:  0 
#> 
#> Minus the Curvature Matrix: 
#> 
#>             X1      X2      X3  sigma2
#> X1      8.6921 -0.9348  1.9296 -0.4483
#> X2     -0.9348 20.8951  1.3249 -0.3665
#> X3      1.9296  1.3249 10.3520 -0.5050
#> sigma2 -0.4483 -0.3665 -0.5050 69.9843

#------------------
#  Treatment Effect
#------------------




###-----------------###
### Survival family ###
###-----------------###

# Setting a seed for reproducibility
set.seed(112358)

#-------------------------
# Simulating Survival data
#-------------------------
n    <- 100
beta <- 1:4
p    <- length(beta)
X    <- data.frame(matrix(rnorm(n * p), n, p)) # continuous (normal) variables

## Simulating survival data from Weibull with a predefined censoring rate of 0.3
y <- surv.simulate(Z = list(X), beta = beta, a = 5, b = exp(1.8), u1 = 0.1,
                   cen_rate = 0.3, gen_data_from = "weibul")$D[[1]][, 1:2]

#---------------------------------------
# MAP estimations with "weibul" function
#---------------------------------------
Lambda <- inv.prior.cov(X, lambda = c(0.1, 1), family = 'survival', basehaz = "weibul")
fit2 <- MAP.estimation(y, X, family = 'survival', Lambda = Lambda, basehaz = "weibul")
fit2
#> $theta_hat
#>       X1       X2       X3       X4  omega_1  omega_2 
#> 1.216502 2.433061 3.305429 4.895585 1.411903 1.905060 
#> 
#> $A_hat
#>                 X1          X2          X3          X4     omega_1   omega_2
#> X1       54.813643    8.232601   -3.813501   -5.643886    7.220449  -49.8473
#> X2        8.232601   85.409611   -5.039660  -16.200600   15.679748 -130.8181
#> X3       -3.813501   -5.039660   57.371520    4.560337   -2.868448 -186.1929
#> X4       -5.643886  -16.200600    4.560337   56.028421   23.120562 -262.4108
#> omega_1   7.220449   15.679748   -2.868448   23.120562   77.588076 -229.2848
#> omega_2 -49.847299 -130.818053 -186.192918 -262.410849 -229.284799 2656.6802
#> 
#> $sd
#>         X1         X2         X3         X4    omega_1    omega_2 
#> 0.15896950 0.19980633 0.26038338 0.36632933 0.14912596 0.06994841 
#> 
#> $Lambda
#>          X1  X2  X3  X4 omega_1 omega_2
#> X1      0.1 0.0 0.0 0.0       0       0
#> X2      0.0 0.1 0.0 0.0       0       0
#> X3      0.0 0.0 0.1 0.0       0       0
#> X4      0.0 0.0 0.0 0.1       0       0
#> omega_1 0.0 0.0 0.0 0.0       1       0
#> omega_2 0.0 0.0 0.0 0.0       0       1
#> 
#> $formula
#> [1] "Survival(time, status) ~ X1 + X2 + X3 + X4"
#> 
#> $names
#> [1] "X1"      "X2"      "X3"      "X4"      "omega_1" "omega_2"
#> 
#> $n
#> [1] 100
#> 
#> $np
#> [1] 4
#> 
#> $treatment
#> NULL
#> 
#> $zero_sample_cov
#> NULL
#> 
#> $refer_cat
#> NULL
#> 
#> $zero_cat
#> NULL
#> 
#> $value
#> [1] -61.88167
#> 
#> $family
#> [1] "survival"
#> 
#> $basehaz
#> [1] "weibul"
#> 
#> $intercept
#> [1] FALSE
#> 
#> $convergence
#> [1] 0
#> 
#> $control
#> $control$maxit
#> [1] 100
#> 
#> 
#> attr(,"class")
#> [1] "bfi"
summary(fit2, cur_mat = TRUE)
#> 
#> Summary of the local model:
#> 
#>    Formula: Survival(time, status) ~ X1 + X2 + X3 + X4 
#>     Family: ‘survival’ 
#>   Baseline: ‘weibul’
#> 
#> Coefficients:
#> 
#>         Estimate Std.Dev CI 2.5% CI 97.5%
#> X1        1.2165  0.1590  0.9049   1.5281
#> X2        2.4331  0.1998  2.0414   2.8247
#> X3        3.3054  0.2604  2.7951   3.8158
#> X4        4.8956  0.3663  4.1776   5.6136
#> omega_1   1.4119  0.1491  1.1196   1.7042
#> omega_2   1.9051  0.0699  1.7680   2.0422
#> 
#> log Lik Posterior:  61.88 
#>       Convergence:  0 
#> 
#> Minus the Curvature Matrix: 
#> 
#>               X1        X2        X3        X4   omega_1   omega_2
#> X1       54.8136    8.2326   -3.8135   -5.6439    7.2204  -49.8473
#> X2        8.2326   85.4096   -5.0397  -16.2006   15.6797 -130.8181
#> X3       -3.8135   -5.0397   57.3715    4.5603   -2.8684 -186.1929
#> X4       -5.6439  -16.2006    4.5603   56.0284   23.1206 -262.4108
#> omega_1   7.2204   15.6797   -2.8684   23.1206   77.5881 -229.2848
#> omega_2 -49.8473 -130.8181 -186.1929 -262.4108 -229.2848 2656.6802

#-------------------------------------
# MAP estimations with "poly" function
#-------------------------------------
Lambda <- inv.prior.cov(X, lambda = c(0.1, 1), family = 'survival', basehaz = 'poly')
fit3 <- MAP.estimation(y, X, family = "survival", Lambda = Lambda, basehaz = "poly")
# Degree of the exponentiated polynomial baseline hazard
fit3$q_l + 1
#> [1] 3
# theta_hat for (beta_1, ..., beta_p, omega_0, ..., omega_{q_l})
fit3$theta_A_poly[,,1][,fit3$q_l+1] # equal to fit3$theta_hat
#>         X1         X2         X3         X4    omega_0    omega_1    omega_2 
#>  0.4587444  0.9065738  1.2387261  1.8013492 -2.0743547  4.7662933 -1.0941826 
# A_hat
fit3$theta_A_poly[,,fit3$q_l+2] # equal to fit3$A_hat
#>               [,1]      [,2]       [,3]       [,4]      [,5]       [,6]
#> X1      67.2578299  3.698096  -6.328249  -0.830014  7.296302   1.005843
#> X2       3.6980959 85.034720   3.130474  -6.311998 15.832762  -1.055534
#> X3      -6.3282486  3.130474  73.449859  11.351621 -2.661869 -25.299425
#> X4      -0.8300140 -6.311998  11.351621  68.043358 23.430052  -8.113705
#> omega_0  7.2963020 15.832762  -2.661869  23.430052 90.074528  43.191997
#> omega_1  1.0058432 -1.055534 -25.299425  -8.113705 43.191997  57.713003
#> omega_2 -0.7626327 -7.194974 -42.325003 -19.174732 48.713003  71.822389
#>                [,7]
#> X1       -0.7626327
#> X2       -7.1949742
#> X3      -42.3250033
#> X4      -19.1747319
#> omega_0  48.7130029
#> omega_1  71.8223891
#> omega_2 122.9402820
summary(fit3, 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.4587  0.1242  0.2153   0.7022
#> X2        0.9066  0.1159  0.6793   1.1338
#> X3        1.2387  0.1364  0.9713   1.5061
#> X4        1.8013  0.1460  1.5153   2.0874
#> omega_0  -2.0744  0.1647 -2.3971  -1.7516
#> omega_1   4.7663  0.2867  4.2044   5.3282
#> omega_2  -1.0942  0.1826 -1.4520  -0.7364
#> 
#> log Lik Posterior:  -1.084 
#>       Convergence:  0 
#> 
#> Minus the Curvature Matrix: 
#> 
#>              X1      X2       X3       X4 omega_0  omega_1  omega_2
#> X1      67.2578  3.6981  -6.3282  -0.8300  7.2963   1.0058  -0.7626
#> X2       3.6981 85.0347   3.1305  -6.3120 15.8328  -1.0555  -7.1950
#> X3      -6.3282  3.1305  73.4499  11.3516 -2.6619 -25.2994 -42.3250
#> X4      -0.8300 -6.3120  11.3516  68.0434 23.4301  -8.1137 -19.1747
#> omega_0  7.2963 15.8328  -2.6619  23.4301 90.0745  43.1920  48.7130
#> omega_1  1.0058 -1.0555 -25.2994  -8.1137 43.1920  57.7130  71.8224
#> omega_2 -0.7626 -7.1950 -42.3250 -19.1747 48.7130  71.8224 122.9403

#------------------------------------------------------
# MAP estimations with "pwexp" function with 3 intervals
#-------------------------------------------------------
# Assume we have 4 centers
Lambda <- inv.prior.cov(X, lambda = c(0.1, 1), family = 'survival',
                        basehaz = 'pwexp', n_intervals = 3)
# For this baseline hazard ("pwexp"), we need to know
# maximum survival times of the 3 other centers:
max_times <- c(max(rexp(30)), max(rexp(50)), max(rexp(70)))
# Minimum of the maximum values of the survival times of all 4 centers is:
min_max_times <- min(max(y$time), max_times)
fit4 <- MAP.estimation(y, X, family = "survival", Lambda = Lambda, basehaz = "pwexp",
                       n_intervals = 3, min_max_times=max(y$time))
#> 
#>  No. observations in the intervals :  73 18 9 
#>  
summary(fit4, cur_mat = TRUE)
#> 
#> Summary of the local model:
#> 
#>    Formula: Survival(time, status) ~ X1 + X2 + X3 + X4 
#>     Family: ‘survival’ 
#>   Baseline: ‘pwexp’
#> 
#> Coefficients:
#> 
#>         Estimate Std.Dev CI 2.5% CI 97.5%
#> X1        0.3310  0.1353  0.0658   0.5961
#> X2        0.6246  0.1242  0.3812   0.8680
#> X3        0.7795  0.1377  0.5097   1.0493
#> X4        1.1811  0.1588  0.8699   1.4924
#> omega_1  -0.3968  0.1787 -0.7471  -0.0465
#> omega_2   1.4806  0.2910  0.9104   2.0509
#> omega_3   1.4282  0.6395  0.1747   2.6817
#> 
#> log Lik Posterior:  -27.06 
#>       Convergence:  0 
#> 
#> Minus the Curvature Matrix: 
#> 
#>              X1      X2       X3      X4 omega_1  omega_2 omega_3
#> X1      58.2386  2.0588  -8.1905 -0.7741  7.5573  -0.4890  0.2409
#> X2       2.0588 74.7939   1.2829 -3.3496 16.9213  -0.5449 -0.5155
#> X3      -8.1905  1.2829  65.1478  9.3205  9.5898 -10.3783 -1.8274
#> X4      -0.7741 -3.3496   9.3205 64.3901 31.2096  -6.2708 -1.4469
#> omega_1  7.5573 16.9213   9.5898 31.2096 54.3965   0.0000  0.0000
#> omega_2 -0.4890 -0.5449 -10.3783 -6.2708  0.0000  14.5194  0.0000
#> omega_3  0.2409 -0.5155  -1.8274 -1.4469  0.0000   0.0000  2.5716


#--------------------------
# Semi-parametric Cox model
#--------------------------
Lambda <- inv.prior.cov(X, lambda = c(0.1), family = 'survival', basehaz = "unspecified")
fit5 <- MAP.estimation(y, X, family = 'survival', Lambda = Lambda, basehaz = "unspecified")
summary(fit5, cur_mat = TRUE)
#> 
#> Summary of the local model:
#> 
#>    Formula: Survival(time, status) ~ X1 + X2 + X3 + X4 
#>     Family: ‘survival’ 
#>   Baseline: ‘unspecified’
#> 
#> Coefficients:
#> 
#>    Estimate Std.Dev CI 2.5% CI 97.5%
#> X1   1.3040  0.1912  0.9293   1.6787
#> X2   2.4894  0.2755  1.9494   3.0293
#> X3   3.2086  0.3354  2.5512   3.8660
#> X4   5.0264  0.5126  4.0218   6.0311
#> 
#> log Lik Posterior:  -143.5 
#>       Convergence:  0 
#> 
#> Minus the Curvature Matrix: 
#> 
#>          X1       X2       X3       X4
#> X1  51.8496   3.0205  -6.6147 -10.8334
#> X2   3.0205  57.7207 -11.4154 -21.3849
#> X3  -6.6147 -11.4154  35.0165 -12.4757
#> X4 -10.8334 -21.3849 -12.4757  23.4777

#------------------
#  Treatment Effect
#------------------