H-likelihood based Hierarchical Joint Models
This package fits shared parameter models for the joint modeling of longitudinal data and survival data, where the longitudinal responses may be of mixed types, such as binary and continuous, and may be left censored by lower limit of quantification. For statistical inference, we consider a computationally efficient approximate likelihood method based on h-likelihood method [3]. There is an extensive literature on h-likelihood method (e.g. [1]-[5]). Essentially, the h-likelihood method uses Laplace approximations to the intractable integral in the likelihood. Moreover, it can produce approximate MLEs for the mean parameters and approximate restricted maximum likelihood estimates (REML) for the variance-covariance (dispersion) parameters.
We also implement the adaptive Gauss-Hermite method to compare with the h-likelihood method.
A detailed example is given in the example folder.
# fit joint model
JMfit ( glmeObject = list( ), survObject = list( ), long.data, surv.data, idVar, eventTime, survFit, method,
itertol=0.001, Ptol=0.01, epsilon=1e-6, iterMax=10, ghsize=4, Silent = T )
# estimate SEs of h-likelihood based parameter estimates by using the adaptive GH method;
# can be used if method='h-likelihood' in JMfit()
JMsd_aGH( JMobject, ghsize=4, srcpath=NULL, paralle=F )
# print coefficient table
JMsummary( JMobject, newSD=NULL, digits=3 )
| glmeObject | A list, indicating the GLME models to be fitted. Details |
| survObject | A list, indicating the survival model (either Cox PH or Weibull model) to be fitted. Details |
| long.data | longitudinal data containing the variables named in formulas in glmeObject |
| surv.data | survival data containing the variables named in formulas in survObject |
| idVar | subject id |
| eventTime | observed event time |
| survFit | an object returned by the coxph() or survreg() function to represent a fitted survival model from the two-step method. |
| method | a vector indicating which method to apply. If method='h-likelihood' (by default), the h-likelihood method is used; if method='aGH', the adaptive Gauss-Hermite method is used. |
| itertol | Convergence tolerance on the relative absolute change in log-likelihood function between successive iterations. Convergence is declared when the change is less than itertol. Default is itertol = 0.001. |
| Ptol | Convergence tolerance on the average relative absolute change in parameter estimates between successive iterations. Convergence is declared when the change is less than Ptol. Default is Ptol = 0.01. |
| epsilon | A small numerical value, used to calculate the numerical value of the derivative of score function. The default value is 1e-6. |
| iterMax | The maximum number of iterations. The default value is 10. |
| ghsize | The number of quadrature points used in the adaptive GH method. The default value is 4. |
| Silent | logical: indicating if messages about convergence success or failure should be suppressed. |
| JMobject | output of JMfit() |
| srcpath | a character vector of full path names indicating the location of the R code; needed if parallel=T; srcpath=NULL by default. |
| parallel | logical: indicating if compute the standard errors of parameter estimates in parallel. By default, paralle=F. |
| newSD | If newSD=NULL (by default), the p-values of the parameter estimates are calculated based on the SEs in the output of JMfit() i.e. fixedsd. Otherwise, the p-values are calculated based on the new SEs given by newSD, e.g. the output of JMsd_aGH(). |
| fixedest | a named vector of estimated coefficients |
| fixedsd | standard errors of estimated coefficients, calculated based on the given method. |
| Bi | estimated random effects, corresponding to each subject |
| B | estimated random effects, corresponding to each measurement |
| covBi | estimated covariance matrix of the random effects |
| sigma | estimates of dispersion parameters |
| convergence | An integer code indicating type of convergence: 0 indicates successful convergence, 1 indicates that the maximum limit for iterations 'iterMax' has been reached without convergence. |
| loglike_value | value of approximate log likelihood function |
| ... | ... |
[1] Do Ha, I., Lee, Y., & Song, J. K. (2002). Hierarchical-likelihood approach for mixed linear models with censored data. Lifetime data analysis, 8(2), 163-176.
[2] Ha, I. D., Park, T., & Lee, Y. (2003). Joint modelling of repeated measures and survival time data. Biometrical journal, 45(6), 647-658.
[3] Lee, Y., Nelder, J. A., & Noh, M. (2007). H-likelihood: problems and solutions. Statistics and Computing, 17(1), 49-55.
[4] Lee, Y., & Nelder, J. A. (1996). Hierarchical generalized linear models. Journal of the Royal Statistical Society. Series B (Methodological), 619-678.
[5] Noh, M., & Lee, Y. (2007). REML estimation for binary data in GLMMs. Journal of Multivariate Analysis, 98(5), 896-915.