Goodness-of-fit tests in proportional hazards models with random effects

Received: 19 November 2020 Revised: 23 December 2021 Accepted: 20 February 2022 DOI: 10.1002/bimj.202000353 RESEARCH ARTICLE Goodness-of-fit tests in proportional hazards models with random effects Wenceslao González-Manteiga1María Dolores Martínez-Miranda2 Ingrid Van Keilegom3 1Department of Statistics, Mathematical Analysis and Operational Research, University of Santiago de Compostela, Santiago de Compostela, Spain 2Department of Statistics and Operations Research, University of Granada, Granada, Spain 3ORSTAT, KU Leuven, Leuven, Belgium Correspondence María Dolores Martínez-Miranda, Department of Statistics and Operations Research, University of Granada, Campus Fuentenueva s/n, 18071 Granada, Spain. Email: [email protected] Funding information Spanish Ministry of Economy and Competitiveness, Grant/Award Number: MTM2016-76969P; European Regional Development Fund; Ministerio de Ciencia e Innovación, Grant/Award Numbers: PID2020-116587GB-I00, MCIN/AEI/10.13039/501100011033; H2020 European Research Council, Grant/Award Number: 694409; University of Granada / CBUA This article has earned an open data badge “Reproducible Research”for making publicly available the code necessary to reproduce the reported results. The results reported in this article were partially reproduced due to data confidentiality issues and their computational complexity. Abstract This paper deals with testing the functional form of the covariate effects in a Cox proportional hazards model with random effects. We assume that the responses are clustered and incomplete due to right censoring. The estimation of the model under the null (parametric covariate effect) and the alternative (nonparametric effect) is performed using the full marginal likelihood. Under the alternative, the nonparametric covariate effects are estimated using orthogonal expansions. The test statistic is the likelihood ratio statistic, and its distribution is approximated using a bootstrap method. The performance of the proposed testing procedure is studied through simulations. The method is also applied on two real data sets one from biomedical research and one from veterinary medicine. KEYWORDS Cox regression, frailty, orthogonal polynomials, proportional hazards, random effects This is an open access article under the terms of the Creative Commons Attribution-NonCommercial-NoDerivs License, which permits use and distribution in any medium, provided the original work is properly cited, the use is non-commercial and no modifications or adaptations are made. © 2022 The Authors. Biometrical Journal published by Wiley-VCH GmbH. Biometrical Journal. 2023;65:2000353. www.biometrical-journal.com 1of22 https://doi.org/10.1002/bimj.202000353 2of22 MANTEIGA et al. 1 INTRODUCTION In survival analysis the Cox proportional hazards model (Cox, 1972) is without any doubt the most classical and wellknown semiparametric model to study the relationship between a survival time 𝑌and a 𝑑-dimensional covariate 𝑋. Under the Cox model the hazard function given 𝑋is described as: 𝜆(𝑡|𝑋)=𝜆 0(𝑡) exp(𝛽′𝑋), where 𝜆0(𝑡) is the unspecified baseline hazard and 𝛽=(𝛽 1,…,𝛽 𝑑)′is a 𝑑-dimensional vector of regression coefficients. The popularity of the Cox model can, among others, be explained by its easy interpretation and the fact that the nonparametric baseline hazard function cancels out in the likelihood, which makes the estimation of the parametric components of the model as simple as in a purely parametric model, like the accelerated failure time model. However, the simplicity of the Cox model cannot accommodate some common patterns of survival data such as correlated data due to the presence of clusters. In these cases the event time data in the same clusters may be correlated since they may be more similar than data from different clusters, although event time data from different clusters are still assumed to be independent. Some examples include multicenter and large-scale medical studies (where the patients’ survival rates may differ substantially across hospitals but may be similar within the same hospital), studies with repeated measurements (where the event data on the same individual may be correlated and therefore each individual may be viewed as a cluster), and recurrent event data (where each individual has several outcomes representing gap times between events). The latter situation is illustrated in this paper with events being recurring infections. In the analysis of clustered survival data, one approach is to introduce random effects to represent the cluster effects. This leads to survival models with random effects and the following extension of the Cox model: 𝜆(𝑡|𝑋𝑖𝑗,𝑏 𝑖)=𝜆 0(𝑡) exp {𝛽′𝑋𝑖𝑗 +𝑏 ′ 𝑖𝑍𝑖𝑗}.(1) Here, 𝜆(⋅|𝑋𝑖𝑗,𝑏 𝑖)is the conditional hazard function of the 𝑗th observation in the 𝑖th cluster (𝑗=1,…,𝑛 𝑖; 𝑖=1,…,𝑞), 𝜆0(⋅) ≥0is the unspecified baseline hazard, 𝑍𝑖𝑗 is a subvector of (1, 𝑋′ 𝑖𝑗)′of dimension 𝑠, and the 𝑏𝑖sare𝑠-dimensional random vectors of mean zero coming from a parametric distribution depending on a parameter vector 𝜃and corresponding to unknown and unobservableP random effects. The unknown random effects in the model play two roles: (1) They incorporate the correlation between the survival data within the same clusters, because the data in the same cluster share the same random effect; (2) they incorporate the variation in the survival data between different clusters since survival data may vary substantially across clusters. Model (1) is often called in the literature the multivariate frailty model. A special case is the so-called shared frailty model, which only involves one random effect and 𝑍𝑖𝑗 =1(i.e., a random intercept). This model is formulated as: 𝜆(𝑡|𝑋𝑖𝑗,𝑏 𝑖)=𝜆 0(𝑡) exp {𝛽′𝑋𝑖𝑗 +𝑏 𝑖}=𝑣 𝑖𝜆0(𝑡) exp {𝛽′𝑋𝑖𝑗},(2) where 𝑣𝑖=exp(𝑏 𝑖)is an unobserved random effect common to all individuals from cluster 𝑖(the shared frailty). Although model (1) is considerably more flexible than the basic Cox model without random effects, it might in some situations not be appropriate in practice. The assumed linearity of the regression function, 𝛽′𝑋𝑖𝑗, might be too strong for one or even several covariates in the model. Nevertheless, linearity is often assumed without any formal or even informal verification. Yu et al. (2012) recognize the importance of having a flexible, not necessarily linear structure and use penalized splines to capture the correct effect of the covariates on the conditional hazard. Previously Yu and Lin (2008)andYuand Lin (2010) used kernel smoothing. In this paper, we develop a test for the adequacy of the linear covariate effects in the shared frailty model (2) and the multivariate frailty model (1). We formulate the following testing problem: 𝐻0∶𝜆 (𝑡|𝑋𝑖𝑗,𝑏 𝑖)=𝜆 0(𝑡) exp {𝑑 ∑ 𝑘=1 𝛽𝑘𝑋𝑖𝑗𝑘 +𝑏 ′ 𝑖𝑍𝑖𝑗},(3) 𝐻1∶𝜆 (𝑡|𝑋𝑖𝑗,𝑏 𝑖)=𝜆 0(𝑡) exp {𝑝 ∑ 𝑘=1 𝛽𝑘𝑋𝑖𝑗𝑘 + 𝑑 ∑ 𝑘=𝑝+1 𝑚𝑘(𝑋𝑖𝑗𝑘)+𝑏 ′ 𝑖𝑍𝑖𝑗},(4) 15214036, 2023, 1, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/bimj.202000353 by Universidade de Santiago de Compostela, Wiley Online Library on [27/03/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License MANTEIGA et al. 3of22 for some 0≤𝑝≤𝑑−1given, where 𝑚𝑘(⋅) (𝑘=𝑝+1,…,𝑑) are nonparametric and nonlinear functions of the covariate 𝑋𝑖𝑗𝑘, which for identifiability reasons are supposed to have mean zero. We refer to model (4) as the nonlinear multivariate frailty model. As far as we know, no paper has considered so far the problem of studying the appropriateness of model (3) versus (4), whereas there is an abundant literature on testing linearity in other models (see, e.g., González-Manteiga & Crujeiras, 2013). For the Cox model, Gray (1994) suggests a likelihood-based test computed from a penalized partial likelihood (based on cubic splines with fixed knots). Lin et al. (2008) formulate a score test for the proportional hazards assumption and for the covariate effects (based on splines). Pan et al. (2015) formulate an additive hazards model with random effects and propose a class of diagnostics to assess the adequacy of the additive representation. For the Cox model with random effects, Xu et al. (2009) use a profile–Akaike information criterion (AIC) and a profile-likelihood ratio test for model selection in the multivariate frailty model, testing for the significance of a specified subset of random or fixed effects. The rest of paper is organized as follows. In the next section the estimation of the model under the null and the alternative hypothesis is explained. In Section 3we present the likelihood ratio test, whereas in Section 4simulations are carried out to study the finite sample performance of the proposed test. Two applications with real data are presented in Section 5. Finally, Section 6contains a discussion on future research directions in this area. 2 ESTIMATION Before attacking the problem of testing the null hypothesis (3) versus the alternative (4), we need to consider the estimation of the model under these two hypotheses. In the context of clustered survival data, we assume the situation where the available information is incomplete due to right censoring. This is, instead of observing the true survival time 𝑌𝑖𝑗, we observe 𝑇𝑖𝑗 =min(𝑌 𝑖𝑗,𝐶 𝑖𝑗), where 𝐶𝑖𝑗 denotes the censoring time. The noncensoring indicators are 𝛿𝑖𝑗 =𝐼(𝑇 𝑖𝑗 =𝑌 𝑖𝑗). The right censored sample is then described by triplets {(𝑇𝑖𝑗,𝑋 𝑖𝑗,𝛿 𝑖𝑗); 𝑗 = 1, … , 𝑛𝑖;𝑖=1,…,𝑞}. We assume that all 𝑋𝑖𝑗s are identically distributed, and 𝑋1,…,𝑋 𝑞are mutually independent, with 𝑋𝑖=(𝑋 𝑖1,…,𝑋 𝑖𝑛𝑖)′. We also assume that the random effects 𝑏1,…,𝑏 𝑞are independent and identically distributed. Finally, the random vectors (𝑇𝑖1,…,𝑇 𝑖𝑛𝑖,𝐶 𝑖1,…,𝐶 𝑖𝑛𝑖,𝑋 𝑖1,…,𝑋 𝑖𝑛𝑖,𝑏 𝑖)are independent (𝑖=1,…,𝑞), the variables 𝑇𝑖𝑗 and 𝐶𝑖𝑗 are independent conditional on 𝑋𝑖𝑗 and 𝑏𝑖,and𝑏𝑖is independent of 𝑋𝑖for 𝑖=1,…,𝑞. 2.1 Estimation under the null model In the (linear) Cox model with random effects, a popular procedure to estimate the model is to maximize the so-called penalized partial likelihood, see, for example, Ripatti and Palmgren (2000), among many other papers. However, although this procedure can be used for estimating the 𝛽-coefficients, the partial likelihood cannot be used for constructing likelihood ratio test statistics, as the partial likelihood is not equal to the full marginal or profile likelihood in the presence of random effects. In particular, when random effects are present, penalized partial likelihoods are not nested, and this means that a likelihood ratio test based on the partial likelihood can lead to negative values, which is clearly not acceptable. Since in this paper we aim at constructing a likelihood ratio test for the null hypothesis (3), we cannot use the partial likelihood and need to work with the full marginal likelihood instead. We follow Vaida and Xu (2000) (see also Abrahantes & Burzykowsky, 2005), who considered the full likelihood of the observed data and the EM algorithm to estimate the parameters in the model. This approach is briefly described below. Under model (1) the interest is to estimate the unknown quantities 𝜉 =(𝛽,𝜃,𝜆 0), and the inference is based on the full likelihood of the observed data, 𝐷={𝐷 𝑖}𝑞 𝑖=1 with 𝐷𝑖={𝑇 𝑖𝑗,𝑋 𝑖𝑗,𝛿 𝑖𝑗}𝑛𝑖 𝑗=1. For cluster 𝑖, conditional on the random effect 𝑏𝑖, the full (conditional) log-likelihood is: 𝑙𝑖0 =𝑙 𝑖(𝛽, 𝜆0;𝐷 𝑖|𝑏𝑖,𝐻 0)= 𝑛𝑖 ∑ 𝑗=1 {𝛿𝑖𝑗𝜆0(𝑇𝑖𝑗)+𝛿 𝑖𝑗(𝛽′𝑋𝑖𝑗 +𝑏 ′ 𝑖𝑍𝑖𝑗) −Λ0(𝑇𝑖𝑗)exp(𝛽′𝑋𝑖𝑗 +𝑏 ′ 𝑖𝑍𝑖𝑗)},(5) where Λ0(𝑡) = ∫𝑡 0𝜆0(𝑠)𝑑𝑠 is the cumulative baseline hazard. The full likelihood is then 𝐿(𝜉;𝐷|𝐻0)= 𝑞 ∏ 𝑖=1 ∫exp(𝑙𝑖0)𝑓(𝑏𝑖|𝜃)d𝑏𝑖,(6) 15214036, 2023, 1, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/bimj.202000353 by Universidade de Santiago de Compostela, Wiley Online Library on [27/03/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 4of22 MANTEIGA et al. where 𝑓(𝑏𝑖|𝜃) is the probability density function of the random effects. Let 𝑙(𝜉;𝐷,𝑏) = ∑𝑞 𝑖=1 𝑙𝑖0 +∑𝑞 𝑖=1 log 𝑓(𝑏𝑖|𝜃).The estimator of the vector (𝛽1,…,𝛽 𝑝,𝜃,𝜆 0)under 𝐻0is now defined as the maximizer of the marginal likelihood (6): (ˆ 𝛽𝐻0,ˆ 𝜃𝐻0,ˆ 𝜆0,𝐻0)=argmax𝛽,𝜃,𝜆0𝐿(𝛽1,…,𝛽 𝑝,𝜃,𝜆 0;𝐷|𝐻0). The EM algorithm provides a convenient method of understanding and computing the maximizers of the likelihood (6), by treating the random effects as “missing.” For the Cox model with random effects, the EM algorithm has the particular feature that the baseline hazard function is maximized nonparametrically, however it can be seen that the parametric EM algorithm applies also to this case. It may be shown that the maximum likelihood estimator of 𝜆0is concentrated at the observed failure times and therefore the maximization problem is equivalent to maximizing a parametric likelihood with parameter vector 𝜆=(𝜆 1,…,𝜆 ℎ)′corresponding to the baseline hazard, where 𝜆𝑙=𝜆 0(𝑡𝑙)(𝑙=1,…,ℎ)and𝑡1,…,𝑡 ℎare the distinct uncensored failures times. With the above definitions the EM algorithm starts with initial parameter values,  𝜉=( 𝛽,  𝜃,  𝜆0), and at each iteration it performs the following two steps until convergence. E-step. This step computes the expectation of the log-likelihood of the augmented data, (𝐷𝑖,𝑏 𝑖), conditional on the observed data, 𝐷𝑖, and the current parameter value,  𝜉=( 𝛽,  𝜃,  𝜆0). This expectation is denoted by 𝑄(𝜉) and given by 𝑄(𝜉) = 𝐸[𝑙(𝜉; 𝐷, 𝑏)|𝐷,  𝜉] ≡𝑄1(𝛽, 𝜆0)+𝑄 2(𝜃), with 𝑄1(𝛽, 𝜆0)= 𝑞 ∑ 𝑖=1 𝑛𝑖 ∑ 𝑗=1 {𝛿𝑖𝑗(log 𝜆0(𝑇𝑖𝑗)+𝛽 ′𝑋𝑖𝑗 +𝐸[𝑏 𝑖| 𝜉]′𝑍𝑖𝑗) −Λ0(𝑇𝑖𝑗)exp(𝛽′𝑋𝑖𝑗 + log 𝐸[exp(𝑏′ 𝑖𝑍𝑖𝑗)|𝐷,  𝜉])},(7) 𝑄2(𝜃) = 𝑞 ∑ 𝑖=1 𝐸[log𝑓(𝑏𝑖|𝜃)|𝐷,  𝜉]. (8) Here, 𝜆0(𝑇𝑖𝑗)equals the corresponding element in the vector 𝜆(defined above) if 𝑇𝑖𝑗 is an observed failure and 0 otherwise, and Λ0(𝑇𝑖𝑗)is the corresponding cumulative baseline hazard. M-step. The estimation of the regression coefficients 𝛽and of the vector 𝜆is conveniently separated from the estimation of the variance components 𝜃, by maximizing the terms 𝑄1and 𝑄2from the E-step. To maximize 𝑄1a profile likelihood approach is used. Assuming for simplicity that there are no ties (otherwise they can be handled in the usual way) the value of 𝜆that maximizes 𝑄1(𝛽, 𝜆), for fixed 𝛽,isˆ 𝜆=( ˆ 𝜆1,…,ˆ 𝜆ℎ)′, with ˆ 𝜆𝑙=1 ∑𝑇ij≥𝑡𝑙exp(𝛽′𝑋ij +𝑢 ij)(𝑙=1,…,ℎ), (9) where 𝑢𝑖𝑗 = log 𝐸[exp(𝑏′ 𝑖𝑍𝑖𝑗)|𝐷,  𝜉]. Substituting (9)into𝑄1gives the profile log-likelihood 𝑞 ∑ 𝑖=1 𝑛𝑖 ∑ 𝑗=1 𝛿𝑖𝑗{𝛽′𝑋𝑖𝑗 −log ∑ 𝑇𝑘𝑙≥𝑇𝑖𝑗 exp(𝛽′𝑋𝑘𝑙 +𝑢 𝑘𝑙)}.(10) Notice that (10) is of the form of the usual log partial likelihood in the Cox model with known offsets 𝑢𝑖𝑗, therefore standard software can be used to obtain the maximizer ˆ 𝛽. Finally 𝑄2is maximized to derive the variance components estimator ˆ 𝜃. To this goal it can be noted that 𝑄2is the log-likelihood corresponding to 𝑞independent observations from the random effects distribution 15214036, 2023, 1, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/bimj.202000353 by Universidade de Santiago de Compostela, Wiley Online Library on [27/03/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License MANTEIGA et al. 5of22 exp 𝐸[log 𝑓(𝑏𝑖|𝜃)|𝐷,  𝜉] and the solution is readily available. For the case of multivariate, Normal random effects with diagonal variance–covariance matrix, 𝑉𝜃=diag(𝜃1,…,𝜃 𝑠), 𝑄2(𝜃)=− 1 2 𝑠 ∑ 𝑔=1 (𝑞 log 𝜃𝑔+1 𝜃𝑔 𝑞 ∑ 𝑖=1 𝐸[𝑏2 𝑖𝑔| 𝜉]), where 𝑏𝑖=(𝑏 𝑖1,…,𝑏 𝑖𝑠)′, and in this case the variance parameter estimates are ˆ 𝜃𝑔=1 𝑞 𝑞 ∑ 𝑖=1 𝐸[𝑏2 ig| 𝜉] (𝑔=1,…,𝑠). The computation of the conditional expectations (7)and(8) is not trivial as in general they are not available in closed form. This means that the E-step involves 𝑠-dimensional numerical integration. For the shared frailty model with Gamma frailty the likelihood in (6) takes a simpler form and the E-step can be performed using closed-form expressions. A simplified algorithm for this case can be found, for example, in Klein (1992) and Nielsen et al. (1992). However this is not the case for other frailty distributions and/or multidimensional random effects (𝑠>1). In practice, to estimate the multivariate frailty model (1) in the general case of 𝑠≥1, we follow Vaida and Xu (2000), where Markov chain Monte Carlo (MCMC) simulation is suggested to approximate the integrals in the E-step. For convenience the authors assume that the random effects are multivariate normally distributed with diagonal variance–covariance matrix. Other distributions for the random effects in the multivariate frailty model can be considered by deriving and properly approximating the corresponding likelihood; see, for example, Korsgaard and Andersen (1999) and Legrand et al. (2005). See also Rondeau et al. (2008)andHaetal.(2011) for other proposals for the estimation of the multivariate frailty model. For a shared frailty model (𝑠=1) with a number of possible frailty distributions, Gorfine et al. (2006) proposed a simpler pseudo full likelihood approach. This method differs from the EM algorithm described above in the way the baseline hazard is estimated. Instead of estimating the cumulative baseline hazard, Λ0(𝑡), maximizing the marginal likelihood (6) at each iteration, it is estimated using a Breslow-type estimator with a noniterative two-stage procedure. Since the standard Breslow formula for the baseline hazard estimator at time 𝑡involves values of Λ0(𝑠), at each iteration, for times beyond time 𝑡, the first stage is thus required to obtain the estimator at those times. The second-stage estimator is the standard Breslow-type estimator that plugs in the first-stage estimator where necessary. The authors summarize their method in the following steps: 1. Obtain initial estimates of 𝛽and 𝜃. 2. Use the current values of 𝛽and 𝜃to estimate Λ0(𝑡) by the two-stage estimator in Gorfine et al. (2006). 3. Substitute the current estimate of Λ0(𝑡) in the likelihood (6) and maximize it to estimate 𝛽and 𝜃. 4. Iterate between Steps 2 and 3 until convergence. Gorfine et al. (2006) show that for the Gamma frailty case their method has essentially the same efficiency as the EM-based maximum likelihood approach. The EM algorithm (including the simplified version above) requires initial values for the parameters  𝜉=( 𝛽,  𝜃,  𝜆0).A good starting point is to consider for  𝛽and  𝜆0the estimates derived from standard Cox regression without random effects. For the variance parameters  𝜃, one can choose the identity vector. 2.2 Estimation under the alternative model Let us consider now the estimation under model (4). A possible way to estimate the nonparametric functions 𝑚𝑘(⋅) (𝑘= 𝑝+1,…,𝑑)wouldbetousepenalizedsplines.Yuetal.(2012) considered this approach in a doubly penalized partial likelihood context (see also Saegusa et al., 2014). These authors showed that their semiparametric model can be written as model (1) with augmented random effects. Using the full likelihood, this approach becomes prohibitive because it requires to integrate out the random effects. This involves multidimensional numerical integration, as we have described in the previous section. Another possible option would be using kernel smoothing. Following ideas commonly used in mixed models, Yu and Lin (2008) suggest local pseudo-partial likelihood to estimate the nonparametric effect of a single covariate using local 15214036, 2023, 1, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/bimj.202000353 by Universidade de Santiago de Compostela, Wiley Online Library on [27/03/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 6of22 MANTEIGA et al. polynomial kernel smoothing. The approach is extended to semiparametric models with time-varying coefficients by Yu and Lin (2010). See also Jin and He (2016) for a related approach using local linear regression. A promising and flexible alternative to the use of penalized splines or kernels is the use of orthogonal expansions. This means that we approximate the functions 𝑚𝑘(𝑥) by an expansion of the form: 𝑟 ∑ 𝓁=0 𝛾𝓁𝑢𝓁(𝑥) for some known orthogonal basis functions 𝑢0,…,𝑢 𝑟. Examples of common basis functions include different types of orthogonal polynomials (like Hermite, Laguerre, Legendre, Chebyshev ,or Bernstein polynomials), orthogonal B-splines, or orthogonal trigonometric functions (see Powell, 1981). The rationale behind using these kinds of approximations is that they can potentially approximate arbitrarily well any continuous function with respect to a certain distance, as long as the number of basis functions 𝑟is taken sufficiently large (see Rice, 1964). They have been used in many contexts, including for goodness-of-fit tests for the random effect in a frailty model based on Laguerre polynomials (see Geerdens et al., 2013), and for approximating the density of a variable with compact support based on Bernstein polynomials (see Bertrand et al., 2019). An important advantage of using orthogonal expansions is that the same estimation approach as under the null model can be used, except that the model now contains 𝑝+∑𝑑 𝑘=𝑝+1 𝑟𝑘coefficients instead of just 𝑑, where 𝑟𝑘is the number of basis functions used to approximate the function 𝑚𝑘. In practice, the number 𝑟𝑘needs to be chosen in an appropriate way, and we choose to use the AIC criterion for this. This entails the fitting of a total of 𝑃𝑑−𝑝 models if we consider candidate values of 𝑟𝑝+1,…,𝑟 𝑑ranging from 1 to some finite 𝑃, and selecting the model with the lowest AIC among these 𝑃𝑑−𝑝 candidate models. In the implementation (see a note on software in Section 4.6) we have considered data-driven shifting and rescaling of the Legendre polynomials in (−1, 1) using the Gram–Schmidt process. This process (also referred to as the Stieltjes process) constructs iteratively the next degree polynomial by removing the components in the directions of the previous ones. For a smooth function 𝑚(𝑥)of a one-dimensional covariate, let us consider the vector of observations of the covariate, 𝑋=(𝑋 1,…,𝑋 𝑛), and the raw polynomials evaluated at the observations, 𝑋𝓁=(𝑋 𝓁 1,…,𝑋𝓁 𝑛),for𝓁 = 0, 1, … , 𝑟. The function 𝑚is then approximated at the observations by ∑𝑟 𝓁=0 𝛾𝓁𝑢𝓁(𝑋), where 𝑢0,…,𝑢 𝑟are orthogonal vectors defined from the following recursion formulas: 𝑢0(𝑋) = 𝑋0= (1, … , 1), 𝑢𝓁(𝑋)=𝑋 𝓁−∑𝓁−1 𝑗=0 <𝑋𝓁,𝑢𝑗(𝑋)> <𝑢𝑗(𝑋),𝑢𝑗(𝑋)> 𝑢𝑗(𝑋) (𝓁 = 1, … , 𝑟), and <⋅,⋅>denotes the inner product, that is, < 𝑢, 𝑣 >= 𝑢′𝑣. For given 𝑟𝑝+1,…,𝑟 𝑑, the conditional log-likelihood under the alternative (4) can now be approximated by 𝑙𝑖1 =𝑙 𝑖(𝜉, 𝛾; 𝐷𝑖|𝑏𝑖,𝐻 1)= 𝑛𝑖 ∑ 𝑗=1 {𝛿𝑖𝑗𝜆0(𝑇𝑖𝑗)+𝛿 𝑖𝑗( 𝑝 ∑ 𝑘=1 𝛽𝑘𝑋𝑖𝑗𝑘 + 𝑑 ∑ 𝑘=𝑝+1 𝑟𝑘 ∑ 𝓁=0 𝛾𝑘,𝓁𝑢𝓁(𝑋𝑖𝑗𝑘)+𝑏 ′ 𝑖𝑍𝑖𝑗) −Λ0(𝑇𝑖𝑗)exp( 𝑝 ∑ 𝑘=1 𝛽𝑘𝑋𝑖𝑗𝑘 + 𝑑 ∑ 𝑘=𝑝+1 𝑟𝑘 ∑ 𝓁=0 𝛾𝑘,𝓁𝑢𝓁(𝑋𝑖𝑗𝑘)+𝑏 ′ 𝑖𝑍𝑖𝑗)},(11) where 𝛾=(𝛾 𝑝+1,0,…,𝛾 𝑑,𝑟𝑑)and the full marginal likelihood is then 𝐿(𝜉,𝛾; 𝐷|𝐻1)= 𝑞 ∏ 𝑖=1 ∫exp(𝑙𝑖1)𝑓(𝑏𝑖|𝜃)d𝑏𝑖.(12) The estimator of the vector (𝛽1,…,𝛽 𝑝,𝛾 𝑝+1,0,…,𝛾 𝑑,𝑟𝑑,𝜃,𝜆 0)under 𝐻1is now defined as the maximizer of the marginal likelihood (12): (ˆ 𝛽𝐻1,ˆ 𝛾𝐻1,ˆ 𝜃𝐻1,ˆ 𝜆0,𝐻1)=argmax𝛽,𝛾,𝜃,𝜆0𝐿(𝛽1,…,𝛽 𝑝,𝜃,𝜆 0,𝛾 𝑝+1,0,…,𝛾 𝑑,𝑟𝑑;𝐷|𝐻1). These estimators are derived again using the EM algorithm described in the previous section. 15214036, 2023, 1, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/bimj.202000353 by Universidade de Santiago de Compostela, Wiley Online Library on [27/03/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License MANTEIGA et al. 7of22 3 THE LIKELIHOOD RATIO TEST We consider a likelihood ratio test for the testing problem formulated in (3)–(4). The test statistic is defined as 𝐿𝑅 = −2{log 𝐿(ˆ 𝛽𝐻0,ˆ 𝜃𝐻0,ˆ 𝜆0,𝐻0|𝐻0)−log𝐿(ˆ 𝛽𝐻1,ˆ 𝛾𝐻1,ˆ 𝜃𝐻1,ˆ 𝜆0,𝐻1|𝐻1)}.(13) To calibrate the test in practice, we use a model-based bootstrap procedure that creates bootstrap samples satisfying the null hypothesis. The bootstrap procedure is an extension of the resampling plans described in Massonnet et al. (2006). It works both in the case of a shared frailty model (where 𝑏𝑖is one-dimensional) and in the case of a multivariate frailty model with an 𝑠-dimensional vector 𝑏𝑖of random effects. It assumes that the censoring time follows a Cox model, but this can be replaced by any other regression model, as long as it can be estimated in a consistent way. From a sample of the form {(𝑇𝑖𝑗,𝑋 𝑖𝑗,𝛿 𝑖𝑗);𝑗=1,…,𝑛 𝑖,𝑖=1,…,𝑞}, bootstrap samples are generated using the following algorithm. Resampling algorithm: 1. Under the null hypothesis (3), fit the model and get the estimators ˆ 𝛽𝐻0,ˆ 𝜃𝐻0,andˆ 𝜆0,𝐻0. 2. Draw i.i.d. random effects 𝑏∗ 𝑖,𝑖=1,…,𝑞, from their distribution with 𝜃replaced by ˆ 𝜃𝐻0. 3. Generate survival times 𝑌∗ 𝑖𝑗 (𝑗=1,…,𝑛 𝑖,𝑖=1,…,𝑞) from the estimated survival function: ˆ 𝑆(⋅) = ˆ 𝑆0(⋅)exp(ˆ 𝛽′ 𝐻0𝑋ij+𝑏∗′ 𝑖𝑍ij), with ˆ 𝑆0(⋅) being the estimated baseline survival function calculated from ˆ 𝜆0,𝐻0in step 1. 4. Generate censoring times 𝐶∗ 𝑖𝑗 (𝑗=1,…,𝑛 𝑖,𝑖=1,…,𝑞) from the Cox-regression estimator of the censoring distribution: ˆ 𝐺(⋅|𝑋ij)= ˆ 𝐺0(⋅)exp(ˆ 𝛿′𝑋ij). 5. Set 𝑇∗ 𝑖𝑗 =min(𝑌 ∗ 𝑖𝑗,𝐶∗ 𝑖𝑗)and 𝛿∗ 𝑖𝑗 =𝐼(𝑇 ∗ 𝑖𝑗 =𝑌 ∗ 𝑖𝑗). The bootstrap sample is then {(𝑇∗ 𝑖𝑗,𝑋 𝑖𝑗,𝛿∗ 𝑖𝑗);𝑗=1,…,𝑛 𝑖,𝑖=1,…,𝑞}. For each bootstrap sample we recalculate the optimal number of basis functions 𝑟∗ 𝑝+1,…,𝑟∗ 𝑑using the AIC, which then finally leads to the bootstrap test statistic 𝐿𝑅∗defined similarly as in (13), except that it is based on the bootstrap data and on the bootstrap estimates of the optimal degrees. This procedure is repeated a large number of times (say 𝐵times) leading to bootstrapped test statistics 𝐿𝑅∗ 1,…,𝐿𝑅∗ 𝐵, and the critical value of the test at level 𝛼is then approximated by the [(1 − 𝛼)𝐵]-th order statistic of these 𝐵values. 4 SIMULATION STUDY 4.1 Design of the study Here we describe the setup of the simulations following the structured approach for planning and reporting simulation studies (ADEMP) proposed by Morris et al. (2019). Aims: (i) To evaluate the type I error and power of our likelihood ratio test (formulated above) for the testing problem (3)–(4). (ii) To compare our test with two possible competitors in terms of type I error and power. (iii) To evaluate the sensitivity of our test to: misspecification of the frailty distribution, varying cluster sizes, and the dimension of the parameters. (iv) To evaluate the performance of our estimator of the nonparametric covariate effect under the alternative, including a comparison with an estimator based on splines. 15214036, 2023, 1, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/bimj.202000353 by Universidade de Santiago de Compostela, Wiley Online Library on [27/03/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 8of22 MANTEIGA et al. 4.1.1 Data-generating mechanisms We consider two different models: a shared frailty model and a multivariate frailty model. In the first case, the (onedimensional) random effect is simulated using a gamma distribution. In the second case, we consider two random effects, which are simulated independently from a normal distribution. In both cases the model contains two (one-dimensional) covariates, which are simulated independently, one (𝑋1) from a Bernoulli distribution (binary) and the other (𝑋2)from a uniform distribution. Under the null hypothesis both covariates have linear effects, represented by two parameters (𝛽1, 𝛽2). Under the alternative hypothesis, the binary covariate has a linear effect (𝛽1) and the other one a nonparametric effect, 𝑚(𝑋2). We consider four types of alternatives: quadratic, sinusoidal, medium-, and high-frequency models. The baseline hazard is set to one in both the shared frailty model and the multivariate frailty model. We generate samples with three sizes (𝑛 = 300, 600, 1200) in two situations: small clusters of size 𝑛𝑖=5and bigger clusters of size 𝑛=20. For aim (iii) above, we consider as well the case of varying cluster sizes. In all cases we generate samples with right censoring, being the censoring percentage of about 40%–70%. For the shared frailty model we simulate censoring times depending on the covariates using a Cox model, while in the multivariate frailty model the censoring is simulated from an exponential distribution, independent of the covariates. 4.1.2 Estimands and targets For aims (i)–(iii) the target is the null hypothesis. For aim (iv) the target is a nonparametric covariate effect under the alternative hypothesis. 4.1.3 Methods We evaluate the likelihood ratio test defined in Section 3for the testing problem (3)–(4). To the extent of our knowledge there are no competitors in the literature for the same problem. Thus, for comparison purposes we define two likelihood ratio tests, which could be possible competitors. First a test that ignores the correlation using simple Cox regression, and second a likelihood ratio test with parametric alternative. For aim (iv) we evaluate our nonparametric estimator using orthogonal polynomials defined in Section 2.2, as well as a competitor proposed by Yu et al. (2012) based on splines. 4.1.4 Performing measures For aims (i)–(iii) we asses the type I error and power of the tests. For assessing the type I error we simulate 500 samples under the null hypothesis and compute the percentage of rejections of the null hypothesis. This defines the empirical level of the test that should be around 5%, which is the assumed nominal level. We compute also the p-value of the test for each sample and report the average from the 500 samples. Under the null hypothesis these p-values follow a uniform distribution with mean 0.5. To evaluate the power we simulate 250 samples under the alternative hypothesis and compute the percentage of rejections of the null hypothesis, which gives the empirical power of the test. Monte Carlo standard errors for the empirical power and level values (100 × ˆ 𝑝) are computed as 100 × √ˆ 𝑝(1 − ˆ 𝑝)∕𝑛𝑠𝑖𝑚, where 𝑛𝑠𝑖𝑚 denotes the number of simulated samples. For aim (iv) we compute the (Monte-Carlo) average of the nonparametric estimators and the corresponding 95% confidence bands, based on the (Monte Carlo) approximated distribution of the estimates. 4.2 Type I error and power of the proposed test 4.2.1 Shared frailty model We first simulate a shared frailty model with Gamma frailty. We evaluate the level of the proposed test by simulating 500 samples under the null hypothesis, and calculating the empirical level defined as the observed proportion of rejections of 15214036, 2023, 1, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/bimj.202000353 by Universidade de Santiago de Compostela, Wiley Online Library on [27/03/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License MANTEIGA et al. 9of22 TABLE 1 Empirical level (%) of the test and Monte Carlo standard errors (between brackets) under the shared frailty model. Our proposal is compared to a simpler test ignoring the correlation. The nominal level is 5% Proposed test Competitor 𝒏𝒏 𝒊𝜽 = 𝟎.𝟓 𝜽 = 𝟐 𝜽 = 𝟎.𝟓 𝜽 = 𝟐 300 5 5.6 5.8 5.2 3.8 (1.0) (1.0) (1.0) (0.9) 20 5.2 6.0 5.6 6.6 (1.0) (1.1) (1.0) (1.1) 600 54.4 4.8 4.4 6.2 (0.9) (1.0) (0.9) (1.1) 20 4.6 5.4 4.8 6.2 (0.9) (1.0) (1.0) (1.1) 1200 5 5.1 5.7 4.6 5.2 (0.9) (1.0) (0.9) (1.0) 20 5.4 5.5 8.2 7.2 (1.0) (1.1) (1.2) (1.2) the null hypothesis. The model under the null hypothesis consists of the following specification of the conditional hazard: 𝜆(𝑡|𝑋ij1,𝑋 ij2,𝑣 𝑖)=𝑣 𝑖𝜆0(𝑡) exp {𝛽1𝑋ij1+𝛽 2𝑋ij2}, 𝑖=1,…,𝑞; 𝑗=1,…,𝑛 𝑖.(14) Here 𝑋1and 𝑋2are independent one-dimensional covariates with linear effects 𝛽1and 𝛽2, respectively, 𝑣𝑖are the frailties with a Gamma distribution with mean 1 and variance 𝜃,and𝜆0(𝑡) = 1 for all 𝑡. We consider values 𝜃=0.5and 2. The covariate 𝑋1has been generated as a binary random variable with an equal probability to be 0 or 1, and the true value of 𝛽1is set to 0.5. The covariate 𝑋2has been generated from a uniform distribution on the interval (0, 1),and𝛽2=1.We simulate samples with 𝑞clusters and 𝑛𝑖observations per cluster, with a total sample size of 𝑛 = 300, 600 and 1200.For each sample size we construct the samples in two cases: small clusters of 𝑛𝑖=5individuals and 𝑞 = 60, 120 and 240; and large clusters of 𝑛𝑖=20individuals and 𝑞 = 15, 30, and 60. In both cases the samples have been generated with right censoring. The censoring time has been generated depending on the covariates with distribution specified by 𝜆𝑐𝑒𝑛(𝑐|𝑋𝑖𝑗1,𝑋 𝑖𝑗2)=0.4exp{0.2𝑋 𝑖𝑗1 +0.5𝑋 𝑖𝑗2}, and the maximum follow-up time has been set to 5. This gives a censoring percentage of about 40%–70%. The results are shown in Table 1for the cases 𝜃 = 0.5, 2,and𝑛𝑖=5and 20. The table shows the percentage of rejections of the null hypothesis for a nominal level of 5%. A comparison with a simple competitor is also shown, which we describe below in Section 4.3. The proposed test has been calibrated using the bootstrap algorithm described in the previous section, considering 𝐵 = 200 bootstrap samples. The table shows that the level of the test with bootstrap calibration is close to the nominal level of 5% for the three sample sizes and all the considered cases. This shows that the test is correctly calibrated for our testing purposes. Moreover the averaged p-values are between 0.47 and 0.52, which agrees with the fact that under the null hypothesis these p-values follow a uniform distribution with mean 0.5. To evaluate the power of the test we have simulated 250 samples from four types of alternatives. In all cases the linear effect of the second covariate in model (14) has been perturbed adding a nonlinear function. The first three alternatives are sinusoidal functions so the second covariate effect is 𝑚𝑏(𝑥2)=𝛽 2𝑥2+𝑎sin(𝑏𝜋𝑥2),for𝑎and 𝑏constants. We consider 𝑏=2, 10, and 20, with the two last ones being more difficult to detect since they are higher frequency models. The last alternative is a quadratic function, 𝑚𝑞(𝑥2)=𝛽 2𝑥2+𝑎𝑥 2 2. Under these alternatives, the test is again performed considering the bootstrap calibration. The percentage of rejections of the null hypothesis is shown in Table 2for the cases 𝜃 = 0.5, 2, and 𝑛𝑖=5and 20. The alternatives in these tables are a sinusoidal (𝑏=2)with𝑎=0.3and 0.5, a mediumand a highfrequency alternatives (𝑏 = 10, 20)with𝑎=1and 1.5; and a quadratic alternative with 𝑎=1and 1.5. Figure 1shows graphics of these alternatives. The results show in general a good performance of the test for all alternatives, with the high-frequency model and the quadratic model being the ones that seem to be harder to detect for our test. As expected the power of the test decreases for 15214036, 2023, 1, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/bimj.202000353 by Universidade de Santiago de Compostela, Wiley Online Library on [27/03/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 16 of 22 MANTEIGA et al. FIGURE 4 Average estimates (with 95% confidence bands) using the splines method of Yu et al. (2012) 4.6 A note on software The computation of our test can be easily performed using standard software for the Cox model with random effects. Our test involves the estimation of the covariate effects under the null and the alternative. Our orthogonal expansion allows to solve both problems using the same tools, noting that under the alternative we will have more coefficients coming from the basis functions. In the computations shown in this paper we have used a basis of orthogonal polynomials constructed from the data using the function poly in R. These orthogonal polynomials are a data-driven shifting and rescaling of the Legendre polynomials on (−1, 1), and they show to be quite stable and work well for all computations here reported. For a given number of basis functions we estimate the model under the alternative as well as the null using maximum likelihood. From the maximized likelihood we estimate the number of basis functions using the AIC criterion. There are several packages available in R to perform these computations (for a good summary see Monaco et al., 2018). In our computations we have used two of these packages: phmm (Donohue & Ronghui, 2017)andfrailtySurv (Monaco et al., 2018). The first one implements the full likelihood and the MCMC-EM algorithm by Vaida and Xu (2000) (see the appendix in that paper to see the details of the MCMC-EM implementation), and it applies to a multivariate frailty model with normal random effects. At the moment this package has been removed from the CRAN site but older versions are still available at https://cran.r-project.org/src/contrib/Archive/phmm/. To the extent of our knowledge it is the only available software to estimate the multivariate frailty model formulated in this paper, apart from frailtypack (Rondeau et al., 2012), which uses a penalized log-likelihood instead. However, the implementation in the package phmm is very time consuming. Moreover, when more than one random effect is considered, it shows also some convergence problems. On the other hand the package frailtySurv implements the pseudo-marginal likelihood approach by Gorfine et al. (2006). It allows for only one random effect with several possible parametric distributions with finite moments. The authors (Monaco et al., 2018) showed that the pseudo-marginal likelihood approach provides quite similar results to the full likelihood with the EM algorithm, and computations are considerably less demanding. For our testing proposals both phmm and frailtySurv provide very similar results when both apply (for a comparison see the results shown for the chronic granulotomous disease [CGD] data below). In the simulation results reported above we have used frailtySurv for the shared frailty model and phmm for the multivariate frailty model with two normal random effects. At early stages of this work we also used the package frailtyEM (Balan & Putter, 2017). This implements the full likelihood with the EM algorithm for the shared frailty 15214036, 2023, 1, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/bimj.202000353 by Universidade de Santiago de Compostela, Wiley Online Library on [27/03/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License MANTEIGA et al. 17 of 22 TABLE 8 Estimated coefficients from the CGD data under the null hypothesis with standard errors, and limits of 95% confidence intervals trtmt inherit cortico prophy sex hosp1 hosp2 hosp3 age ˆ 𝛽1.14 0.82 −1.97 0.95 −0.96 −0.27 −1.09 −0.94 −0.04 SD 0.34 0.37 0.96 0.46 0.51 0.40 0.61 0.59 0.02 Lower 0.47 0.09 −3.85 0.04 −1.96 −1.06 −2.30 −2.10 −0.08 Upper 1.80 1.56 −0.08 1.86 0.04 0.52 0.11 0.21 −0.01 model only and allows for several distributions. But the results were quite similar to those from frailtySurv with larger computational time. R-code using the package frailtySurv is available with the paper as Supporting Information. It allows to reproduce the results reported in our first data application below. The computations of the test take about 5 min (4.463 min with a processor Intel(R) Core(TM) i7-10710U CPU 1.10 GHz 1.61 GHz) using 𝑃=7and 𝐵 = 200 bootstrap samples, and above 10 min using 𝐵 = 500, which is recommended to have an accurate approximation of the 95% quantile of the test statistic. 5 APPLICATIONS WITH REAL DATA In this section we illustrate the proposals in the paper using two data examples. The first one comes from biomedical research and the second one from veterinary medicine. 5.1 Chronic granulotomous disease The data come from a placebo-controlled randomized trial of gamma interferon in patients with CGD. There were a total number of 𝑛 = 203 records collected from 𝑞 = 128 patients. For each patient the time between entry into the study and first infection, and the times between any recurrent infections are given. The data are provided in Appendix D.2 of Fleming andHarrington(1991), as well as in the frailtySurv package (Monaco et al., 2018). These recurrent events data have been analyzed by various authors, see Vaida and Xu (2000) for references and a discussion of previous results. Vaida and Xu (2000) modeled the data using the parametric shared frailty model in (3), with a log-normal frailty. In this model each patient represents a cluster, and for the 𝑖th patient (𝑖=1,…,𝑞), the times 𝑇𝑖𝑗 (𝑗=1,…,𝑛 𝑖) are the censored/uncensored (correlated) infection times for the patient since the last infection or since study entry, whichever occurred later. As Vaida and Xu (2000) argued, such a model may be appropriate if the baseline hazard rate of the backward recurrence time is unaffected by the number of previous events, and practical experience seems to support this assumption. We consider here the same model assumptions as Vaida and Xu (2000), where the frailty has a log-normal distribution, with the exception that we allow for nonparametric covariate effects as in model (4). We use the same nine covariates as these authors, namely, treatment, pattern of inheritance, age, whether or not using corticosteroids at study entry, whether or not using prophylactic antibiotics at study entry, sex, and hospital category (four categories from which three binary covariates are created). Only one of these covariates, age, is continuous, and the others are binary. Our aim is to test whether age has indeed a linear effect on the hazard as these authors assumed. Therefore we formulate the testing problem as follows: 𝐻0∶𝜆 (𝑡|𝑋𝑖𝑗,𝑏 𝑖)=𝜆 0(𝑡) exp {𝛽′𝑋𝑖𝑗 +𝑏 𝑖},(16) 𝐻1∶𝜆 (𝑡|𝑋𝑖𝑗,𝑏 𝑖)=𝜆 0(𝑡) exp {8 ∑ 𝑘=1 𝛽𝑘𝑋𝑖𝑗𝑘 +𝑚 9(𝑋𝑖𝑗9)+𝑏 𝑖},(17) where 𝑋𝑖𝑗 is the vector of covariate values for the 𝑗th event of the 𝑖th individual, and 𝑚9(⋅) is a nonparametric function of the covariate age. To estimate the parameters in the model under the null and the alternative, we use the phmm package (Donohue & Ronghui, 2017). The estimated effects of the covariates under the null hypothesis are shown in Table 8. The estimated 15214036, 2023, 1, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/bimj.202000353 by Universidade de Santiago de Compostela, Wiley Online Library on [27/03/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 18 of 22 MANTEIGA et al. TABLE 9 Estimated coefficients from the CGD data under the alternative hypothesis with standard errors and limits of 95% confidence intervals trtmt inherit cortico prophy sex hosp1 hosp2 hosp3 age𝟏age𝟐age𝟑 ˆ 𝛽1.03 1.00 −1.94 1.13 −1.08 −0.28 −1.18 −0.91 −12.56 −12.50 −10.81 SD 0.31 0.37 0.78 0.43 0.49 0.39 0.59 0.57 4.87 6.29 4.92 Lower 0.41 0.27 −3.48 0.28 −2.03 −1.04 −2.34 −2.02 −22.11 −24.83 −20.46 Upper 1.64 1.73 −0.40 1.98 −0.12 0.48 −0.03 0.21 −3.01 −0.17 −1.16 FIGURE 5 Estimated effects under the null and the alternative for the data applications: CGD data on the left and culling data on the right FIGURE 6 Estimated baseline survival function under the null and the alternative hypotheses (CGD data) variance of the random effect is 0.67. These results are very close to those reported in table II of Vaida and Xu (2000). Under the alternative hypothesis we use the AIC criterion to select the degree of the orthogonal representation of the covariate age. In this case the chosen degree is 3 and the estimated effects of age are shown in Table 9. Figure 5shows a plot of the two fits and Figure 6shows the estimated baseline survival function. The two fits show slight differences apart from the right-hand side (older ages), which corresponds to an area of very few observations. From the survival function plots it can be seen that the two estimators differ slightly, which is in line with the differences shown between the two fits. 15214036, 2023, 1, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/bimj.202000353 by Universidade de Santiago de Compostela, Wiley Online Library on [27/03/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License MANTEIGA et al. 19 of 22 FIGURE 7 Estimated random effects (with 95% confidence intervals) for log-baseline hazard (CGD data) For comparison we also estimated the models using the frailtySurv package (Monaco et al., 2018) and considering both a normal distribution for the random effects 𝑏𝑖and a gamma distribution for the frailty exp(𝑏𝑖). In both cases the results are quite similar for most of the coefficients. Using these estimated coefficients we perform now the test for problem (16)–(17). Under the normal assumption and using the frailtySurv packageweget𝐿𝑅 = 7.17 with p-value 0.10, and 𝐿𝑅 = 7.61 with p-value 0.12 under the gamma assumption. The conclusion is therefore that we do not have evidence to reject the null hypothesis at the 5% significance level, but the results are not conclusive at 10% significance. So the linear representation of age in the model used by Vaida and Xu (2000) could be inappropriate. For comparison we have also computed the same kind of test but ignoring the correlation (see Section 4.3). For that case 𝐿𝑅 = 19.62 with a p-value of 0, which leads to the rejection of the null hypothesis. Recall that we do not recommend such a test for correlated data. Finally we compare the estimated random effects for each individual under the two models. The random effects are estimated as the conditional expectations 𝐸[𝑏𝑖| 𝜉], which are derived from the last step of the likelihood maximization algorithm. Figure 7shows these estimates along with 95% confidence intervals, computed using the approach in Vaida and Xu (2000). The patients (clusters) have been ordered according to the number of records, that is, the total number of infections. This shows that the random effects of patients with different numbers of records are different, as Vaida and Xu (2000) pointed out (see fig. 4 in the paper and related discussion on p. 3320). This is an interesting observation since when the number of infections is included in the shared frailty model as a covariate, it turns out not to be significant. This indicates that previous infections do not increase the risk of future recurrence, but it makes patients different, and these differences are picked up by the frailty in the model. Note that this effect can be observed under both null and alternative hypotheses in Figure 7, however it seems to be a bit smoother under the alternative. 5.2 Culling of dairy heifer cows Here we consider the culling data set described in Example 1.7 in Duchateau and Janssen (2008). The time (in days) to culling is studied in heifers as a function of the somatic cell count (SCC) measured between 5 and 15 days after calving. Cows are followed up for an entire lactation period (roughly 300–350 days) and if they are still alive at the end of the 15214036, 2023, 1, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/bimj.202000353 by Universidade de Santiago de Compostela, Wiley Online Library on [27/03/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 20 of 22 MANTEIGA et al. lactation period they are censored at that time. Cows are further clustered within herds and this clustering needs to be taken into account as culling policy and also SCC in early lactation might differ substantially between the herds. As the authors pointed out, high SCC might be a surrogate marker for intramammary infections and heifers, which have intramammary infections, or which are expected to develop intramammary infections in the future, are quite expensive to keep due to the high costs for drugs and the loss in milk production. For these data we consider the two available covariates to describe the time to culling: the day on which the SCC is assessed (Timeassess) and the logarithm of the SCC on that day (logSCC). In Duchateau and Janssen (2008)afrailty model with gamma frailty has been considered for these data assuming that logSCC has a linear effect. Our goal is to test whether this linear specification of logSCC is correct, and so we formulate the testing problem as follows: 𝐻0∶𝜆 (𝑡|𝑋𝑖𝑗,𝑏 𝑖)=𝜆 0(𝑡) exp {𝛽′𝑋𝑖𝑗 +𝑏 𝑖},(18) 𝐻1∶𝜆 (𝑡|𝑋𝑖𝑗,𝑏 𝑖)=𝜆 0(𝑡) exp {𝛽1𝑋𝑖𝑗1 +𝑚(𝑋 𝑖𝑗2)+𝑏 𝑖},(19) where 𝑋𝑖𝑗 is the vector of the two covariate values for the 𝑗th cow of the 𝑖th herd, 𝑚(⋅) is a nonparametric function of the covariate logSCC, and the frailty exp(𝑏𝑖)has a gamma distribution. First we estimate the model under the null and the alternative hypotheses. Since the data set is quite big (14,246 cows) we consider a random subsample of size 𝑛 = 5000. This sample size is still quite big so computations are demanding, and hence we only use the frailtySurv package (Monaco et al., 2018) to ease the computations. Figure 5shows the estimated effect of logSCC under the null and the alternative hypotheses. Under the alternative, using the AIC criterion we choose 4 basis functions, indicating that the linearity assumption for this covariate might not be appropriate. We check this result computing our test, which gives the value of 𝐿𝑅 = 11.37 and a p-value of 0.08 based on 500 bootstrap samples. At the 10% significance level we reject the null hypothesis and conclude that the covariate logSCC has a nonlinear effect on the time to culling. For comparison we have also computed the same kind of test but ignoring the correlation (see Section 4.3). For this case the estimated model under the alternative chooses 1 basis function (degree 1 as the null) and the test concludes that there is no evidence against the null hypothesis. 6DISCUSSION In this paper we developed a goodness-of-fit test for the functional form of the regression function in a Cox model with random effects. We used an approach based on the marginal likelihood of the model, and under the alternative we used orthogonal expansions to approximate the regression function. Simulations showed that the proposed bootstrap calibration works well in practice. When the linearity of a covariate is rejected by our test, our nonparametric estimator based on orthogonal expansions allows for the required flexibility to properly describe its effect on the survival time. Our nonparametric approach is quite simple compared to other approaches in the literature based on splines (Yu et al., 2012) or kernels (Yu & Lin, 2008). Our orthogonal expansions make the problem of estimating the nonparametric effects as simple as the well-known linear effects, moreover the same likelihood approach can be used. The empirical studies described in this paper consider the simpler situation of testing the linearity of one single covariate. This is convenient for us to speed up the computations and carry out the simulations. Our proposal nevertheless is more general and allows for a semiparametric alternative with more than one nonparametric covariate effect, 𝑚𝑘,𝑘=𝑝+1,…,𝑑(see Section 2.2). As an example consider the situation where the model involves 𝑑=3covariates and we want to test a nonlinear effect for two of them (𝑝=1). In this case the alternative would be 𝜆(𝑡|𝑋𝑖𝑗,𝑏 𝑖)= 𝜆0(𝑡) exp{𝛽1𝑋𝑖𝑗1 +𝑚 2(𝑋𝑖𝑗2)+𝑚 3(𝑋𝑖𝑗3)+𝑏 ′ 𝑖𝑍𝑖𝑗}, and to compute the test we would need to estimate such model, which involves two nonparametric functions 𝑚2and 𝑚3. Our proposal is to use orthogonal representations for each function and determine the number of basis functions for the first one, 𝑟2, and for the second one, 𝑟3, using the AIC criterion. This means that we would fit a total of 𝑃2candidate models with candidate values of 𝑟2and 𝑟3ranging from 1 to some finite but large enough 𝑃(𝑃=7was enough in our computations in the paper, except for the high-frequency model in the simulations where we considered 𝑃=12), selecting the values, which lead to the lowest AIC. Once the model under the alternative is fitted the rest of the computations for the test work in the same way as the situation of only one nonparametric effect. Therefore, the only difference with respect to the latter case is the number of candidate models we need to evaluate (𝑃in the case of one nonparametric effect, 𝑃2for two effects, 𝑃3for three, etc.). 15214036, 2023, 1, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/bimj.202000353 by Universidade de Santiago de Compostela, Wiley Online Library on [27/03/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License MANTEIGA et al. 21 of 22 A number of issues could be studied in more detail in future research. For instance, the use of a nonadditive structure under the alternative by using multidimensional orthogonal expansions would be worth studying in the future. Our testing procedure depends on certain parametric assumptions on the distribution of the random effects, which could be misspecified in practice. Our simulation experiments show that the misspecification of the frailty distribution has a minor impact on the power. An additional issue is the sensitivity of the power to the number of parameters in the model and the dimension of the data. We have described a brief simulation experiment, which shows that increasing the dimension does not have much impact on the power. In our first data application with recurrent events the frailty is assumed to be constant for each patient. In some situations the frailty changes over time with one value for each gap time. Modeling the frailty as a time series with autorregresive correlation structures has been considered in the literature by Yau and McGilchrist (1998), Tawiah et al. (2020), Putter and Van Houwelingen (2015), among others. The extension of our goodness-of-fit test to this model is worth studying in the future. An alternative way of describing the correlation in the proportional hazards model consists of incorporating the random effects in the log relative risk. The resulting models, under a suitable transformation, can be estimated as mixed effects regression models (see for example López-de-Ullibarri et al., 2012, for a recent discussion and references). Under this formulation it would be interesting to explore an extension of the test of González-Manteiga et al. (2016). ACKNOWLEDGMENTS The authors thank three anonymous reviewers and the associate editor for many valuable comments and suggestions, which have helped to improve the quality of the article. We also thank Prof. Zhangsheng Yu for providing sample code for the approach in Yu et al. (2012), and the University of Granada for providing the computing resources. This research was supported by the Spanish Ministry of Science and Innovation through grants PID2020-116587GB-I00 (MCIN/AEI/10.13039/501100011033) and MTM2016-76969P, and the European Research Council (2016-2021, Horizon 2020 / ERC grant agreement No. 694409). Funding for open access charge: University of Granada / CBUA. CONFLICT OF INTEREST The authors have declared no conflict of interest. OPEN RESEARCH BADGES This article has earned an Open Data badge for making publicly available the digitally-shareable data necessary to reproduce the reported results. The data is available in the Supporting Information section. This article has earned an open data badge “Reproducible Research” for making publicly available the code necessary to reproduce the reported results. The results reported in this article were partially reproduced due to data confidentiality issues and their computational complexity. ORCID Wenceslao González-Manteiga https://orcid.org/0000-0002-3555-4623 MaríaDolores Martínez-Miranda https://orcid.org/0000-0001-5468-627X IngridVan Keilegom https://orcid.org/0000-0001-8827-7642 REFERENCES Abrahantes, J. C., & Burzykowski, T. (2005). A version of the EM algorithm for proportional hazard model with random effects. Biometrical Journal,47,847–862. Balan, T. A., & Putter, H. (2017). frailtyEM: An R package for estimating semiparametric shared frailty models. Vignette https://cran.r-project. org/web/packages/frailtyEM/ Bertrand, A., Van Keilegom, I., & Legrand, C. (2019). Flexible parametric approach to classical measurement error variance estimation without auxiliary data. Biometrics,75,297–307. Cox, D. R. (1972). Regression models and life-tables (with discussion). Journal of the Royal Statistical Society - Series B,4,187–220. Donohue, M. C., & Ronghui, X. (2017). Proportional hazards mixed-effects models.R package version 0.7-10.https://cran.r-project.org/src/ contrib/Archive/phmm/ Duchateau, L., & Janssen, P. (2008). The frailty model. Springer-Verlag. Fleming, T. R., & Harrington, D. P. (1991). Counting processes and survival analysis. Wiley. 15214036, 2023, 1, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/bimj.202000353 by Universidade de Santiago de Compostela, Wiley Online Library on [27/03/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 22 of 22 MANTEIGA et al. Geerdens, C., Claeskens, G., & Janssen, P. (2013). Goodness-of-fit tests for the frailty distribution in proportional hazards models with shared frailty. Biostatistics,14,433–446. González-Manteiga, W., & Crujeiras, R. M. (2013). An updated review of goodness-of-fit tests for regression models. Test,22,361–411. González-Manteiga, W., Martínez-Miranda, M. D., & Van Keilegom, I. (2016). Goodness-of-fit test in parametric mixed effects models based on the estimation of the error distribution. Biometrika,103,133–146. Gorfine, M., Zucker, D. M., & Hsu, L. (2006). Prospective survival analysis with a general semiparametric shared frailty model: A pseudo full likelihood approach. Biometrika,93,735–741. Gray, R. J. (1994). Spline based tests in survival analysis. Biometrics,50,640–652. Ha, I. D., Sylvester, R., Legrand, C., & MacKenzie, G. (2011). Frailty modelling for survival data from multi-centre clinical trials. Statistics in Medicine,30, 2144–2159. Jin, Z., & He, W. (2016). Local linear regression on correlated survival data. Journal of Multivariate Analysis,147,285–294. Klein, J. P. (1992). Semiparametric estimation of random effects using the Cox model based on the EM algorithm. Biometrics,48, 795–806. Korsgaard, I. R., & Andersen, A. H. (1998). The additive genetic gamma frailty model. Scandinavian Journal of Statistics,25,255–269. Legrand, C., Ducrocq, V., Janssen, P., Sylvester, R., & Duchateau, L. (2005). A Bayesian approach to jointly estimate centre and treatment by centre heterogeneity in a proportional hazards model. Statistics in Medicine,24, 3789–3804. Lin, J., Zhang, D., & Davidian, M. (2008). Smoothing spline-based score tests for proportional hazards models. Biometrics,62,803–812. López-de-Ullibarri, I., Janssen, P., & Cao, R. (2012). Continuous covariate frailty models for censored and truncated clustered data. Journal of Statistical Planning and Inference,142,1864–1877. Massonnet, G., Burzykowski, T., & Janssen, P. (2006). Resampling plans for frailty models. Communications in Statistics-Simulation and Computation,35,497–514. Monaco, J. V., Gorfine, M., & Hsu, L. (2018). General semiparametric shared frailty model: Estimation and simulation with frailtySurv. Journal of Statistical Software,86(4), 1–42. Morris, T. P., White, I. R., & Crowther, M. J. (2019). Using simulation studies to evaluate statistical methods. Statistics in Medicine,38, 2074–2102. Nielsen, G. G., Gill, R. D., Andersen, P. K., & Sorensen, T. I. (1992). A counting process approach to maximum likelihood estimation in frailty models. Scandinavian Journal of Statistics,19,25–43. Pan, D., Liu, Y. Y., & Wu, Y. S. (2015). Additive hazards regression with random effects for clustered failure times. Acta Mathematica Sinica, English Series,31,511–525. Powell, M. J. D. (1981). Approximation theory and methods. Cambridge University Press. Putter, H., & Van Houwelingen, H. C. (2015). Dynamic frailty models based on compound birth-death processes. Biostatistics,16,550–564. Rice, J. J. (1964). The approximation of functions. Addison - Wesley. Ripatti, S., & Palmgren, J. (2000). Estimation of multivariate frailty models using penalized partial likelihood. Biometrics,56, 1016–1022. Rondeau, V., Michiels, S., Liquet, B., & Pignon, J. P. (2008). Investigating trial and treatment heterogeneity in an individual patient data metaanalysis of survival data by means of the penalized maximum likelihood approach. Statistics in Medicine,27, 1894–1910. Rondeau, R., Mazrouni, Y., & Gonzalez, J. R. (2012). frailtypack: An R package for the analysis of correlated survival data with frailty models using penalized likelihood estimation or parametrical estimation. Journal of Statistical Software,47(4), 1–28. Saegusa, T., Di, C., & Chen, Y. K. (2014). Hypothesis testing for an extended Cox model with time-varying coefficients. Biometrics,70,619–628. Tawiah, R., McLachlan, G. J., & Ng, S. K. (2020). Mixture cure models with time-varying and multilevel frailties for recurrent event data. Statistical Methods in Medical Research,29,1368–1385. Vaida, F., & Xu, R. (2000). Proportional hazards model with random effects. Statistics in Medicine,19, 3309–3324. Xu, R., Vaida, F., & Harrington, D. P. (2009). Using profile likelihood for semiparametric model selection with application to proportional hazards mixed models. Statistica Sinica,19,819–842. Yau, K. K. W., & McGilchrist, C. A. (1998). ML and REML estimation in survival analysis with time dependent correlated frailty. Statistics in Medicine,17, 1201–1213. Yu, Z., & Lin, X. (2008). Nonparametric regression using local kernel estimating equations for correlated failure time data. Biometrika,95, 123–137. Yu, Z., & Lin, X. (2010). Semiparametric regression with time-dependent coefficients for failure time data analysis. Statistica Sinica,10,853–869. Yu, Z., Lin, X., & Tu, W. (2012). Semiparametric frailty models for clustered failure time data. Biometrics,68,429–436. SUPPORTING INFORMATION Additional supporting information can be found online in the Supporting Information section at the end of this article. How to cite this article: González-Manteiga, W., Martínez-Miranda, M. D., & Van Keilegom, I. (2023). Goodness-of-fit tests in proportional hazards models with random effects. Biometrical Journal,65, 2000353. https://doi.org/10.1002/bimj.202000353 15214036, 2023, 1, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/bimj.202000353 by Universidade de Santiago de Compostela, Wiley Online Library on [27/03/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License