Some Genetic Regression Models for Multiple Quantitative End Points Data

Abstract

Multiple endpoints data are common in practice. There are various statistical methods for the analysis of this type of data, however, genetic models for familial observations with multiple endpoints data are relatively few, and the existing methods are basically variations of the Elston-Stewart algorithm. Here we consider several joint statistical models for such data with quantitative measurements with a new algorithm, which is computationally more efficient than the existing method. The proposed method is detailed in some commonly used parametric, semi parametric and nonparametric settings for this type of data. For un-genotyped data, the commonly used models are the mixture and variance components models. We elaborate how these genetic models can be extended for multiple endpoints data with the proposed method

Keywords: Censoring; Endpoints data; Familial structure; Genotype; Missing observation

Introduction

Endpoints data are observed responses from patients of some pre-specified clinical events of interests, such as death, loss of vision, occurrences of certain diseases, or other symptomatic events. In medical research, study participants are often followed for a long time, during which some participants may drop out early, so that random censorship may be present in the data. Such data have missing observations, which may be inhomogeneous across the patients. For example, in one patient we have observations on the lung cancer and kidney disease, and on another we have observations on lung cancer, diabetes and asthma. Analyses of such data largely fall into two categories: hypothesis testing (usually non-model based) and model inference. Here we concentrate on the model inference of such data.

For censored data modeling there are extensive literatures [1-8], just mention a few. For multiple endpoints data analyses, there are various statistical methods [9-13], for example. Wei and Glidden [14] provided an overview for some of the methods in this field.

Family genetic data differ from the ordinary data in that they are collected in familial units, often with varying structures and sizes, and with/without genotyping. These features make the models distinct. The key in the modeling is the familial dependence structure and the implementation of the genetic mechanism, existing methods are basically variations of the Elston-Stewart algorithm, which is a multi-level mixture model, and the computation is often challenging. Genetic models for multiple endpoints data are relatively limited. Here we consider some statistical joint models for such data with quantitative measurements with a new algorithm, which is a one-level mixture model, thus enhance the computation considerably. The parametric method is used when one has some confidence about the model specification. The semi parametric method can be used when there is not enough information about the full parametric model specification. The nonparametric method is used for the robustness of least model assumptions. We elaborate our methods for the parametric, semi parametric and nonparametric cases. The methods we describe below are valid for arbitrary pedigrees; however, in this article we focus on the simpler case of nuclear family for illustration.

In genetic analysis, the data contains genotypes, partially genotypes, or no genotypes. However, even if the data are genotyped, it is still of interest to know whether there are some other unknown gene(s) behind the response functioning. There are reports that with added unknown gene locus, the likelihood Akaike information reduced (e.g., [15], p.1091, [16], which makes sense, as correct parameter(s) added to the model will reduce its AIC), or the segregation analysis guided to some other gene(s) which deserve(s) further investigation. So even for the genotyped data, a segregation model is still of importance. It is also the general model including the genotyped data case. In the following we derive some commonly used regressive models for this case, including the parametric model, semiparametric proportional hazards model, nonparametric least squares model, variance components model and the competing risks model. Also, hypothesis testing on parameters of interest can be conducted using the likelihood ratio statistics based on the parametric models. Our aim here is to present several new parametric, semi parametric and nonparametric models for this familial data, and thus we focus derivations of the basic forms of these models. Implementations of these models and applications to real data will follow in our future work.

Methods

Suppose there are d responses observed with some clinical events of interests, along with r covariates, for each member in a family. We concentrate on nuclear family structure for simplicity. In practice we only observe a subset of the responses and covariates for each patient in a family. Let y_i=(y_if,y_im,y_is)(i=1,…,n) be the vector of responses for the i-th nuclear family, with corresponding covariates x_i=(x_if,x_im,x_is) (i=1,…,n). Here y_if a d_if(<d) dimensional vector of responses of the father in the i-th family, which belongs to a dif dimensional subspace of the d-dimensional space, with a response non-missing indicator vector I_if and covariate non-missing indicator vector J_if. Similarly, y_im denotes a vector of responses of the mother in the i-th family and y_is=(y_i1,…,y_ibi) is for offspring with each of the y_ij s has the same data structure as that for y_if. For example, there are three responses to be observed, we have the first and third on the father, then I_if=(1,0,1), and d_if=|I_if|=2 is its dimension or cardinality. If there are total of five covariates in the design, and we may only have the first, second, fourth covariates for the father, then J_if=(1,1,0,1,0), and r_if=|J_if|=3 is its dimension. We assume random censorship. Let δ_if=(δ_if1,…, δ_ifd_if) be the censoring indicator of y_if, i.e. δ_ifj=1 if y_ifj is uncensored, and δ_ifj=0 otherwise. Similar notations are used for the mother. For the off springs, y_is = (y_i1,...,y_ki) denote the response vector, with y_ij be d_ij dimensional observation for the j^th sib, with response configuration I_ij and covariate configuration j_if, (j=1,…,k_i). Let d_i=(d_if,d_im,d_is), I_i=(I_if,I_im,I_is), J_i=(J_if,J_im,J_is), δ_i=(δ_if,δ_im,δ_is) and δ_i=(C_if,C_im,C_is). The complete data information consists of Z_i=(y_i,x_i,I_i,J_i,δ_i), (i=1,…,n). Let χ(⋅) be the indicator function, i.e. χ(gir=s)=1 if the genotype of the rth individual is s and zero otherwise. Let πs be the population proportion of the S-th genotype, t(s|i,k) be the transmission probability of a off string’s genotype s given the parents’ genotype (i,k), and θ be the collection of all the parameters, including the as and βs in the mean and the parameters in the within and between individual covariance matrices and the genotype frequencies πk_s and transmission probabilities t(s|i,k)_s. With unobserved genotypes, the computation is a serious challenge because of the mixture nature of the model. Let f(y_r|θ), F(y_r|θ), and S(y_r|θ)=1-F(y_r|θ) be the density function, distribution function and survival function of Y_r, respectively. They may be described by its genotype g_r through the mean function specification. To simplify model specification, we assume random mating so that father and mother can be viewed as independent in most cases. The case of non-random mating or within parent’s dependence can be treated similarly with more involved notations.

Note that there is within family dependence but independence among different families. We assume the genotypes of each patient are unobserved; the case of observed genotypes is automatically covered and simpler. Now we describe the methods in some common settings below.

Parametric model

Let the genotypes at the locus of interest be coded as 1,…,k. We first consider the case of no missing record and censoring. The regressive model assumes y_ir=μ(g_ir)+∈_ir, (i,=1,…,n; r=f,m,1,…,n),

where μ(g_ir)=μ0+αχ(g_ir)+βx_ir is the mean phenotypic value, α=(α₁,…,α_k)', χ(g_ir)=(χ(g_ir=1),…,χ(g_ir=k), μ₀ is the intercept vector,. The residual error term ∈_ir is a d dimensional random vector where they are independent across i but dependent across r.

In the case of multivariate familial quantitative response data, under the commonly used Elston-Stewart [15] algorithm or its variants, the likelihood of the observation yi for the i-th nuclear family is

$L_i\mathrm{<mo\; stretchy="false">(</mo>}{y}_i\mathrm{<mo>|</mo>}\theta \mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{\displaystyle \sum _k}\mathrm{<mo\; stretchy="false">(</mo>}{\pi }_{k}f\mathrm{<mo\; stretchy="false">(</mo>}{y}_{if}\mathrm{<mo>|</mo>}\theta \mathrm{,}{g}_{if}\mathrm{<mo>=</mo>}k\mathrm{<mo\; stretchy="false">)</mo>}{\displaystyle \sum _j}\mathrm{<mo\; stretchy="false">(</mo>}{\pi }_{j}f\mathrm{<mo\; stretchy="false">(</mo>}{y}_{im}\mathrm{<mo>|</mo>}\theta \mathrm{,}{g}_{im}\mathrm{<mo>=</mo>}j\mathrm{<mo\; stretchy="false">)</mo>}K\mathrm{<mo\; stretchy="false">(</mo>}\theta \mathrm{,}k\mathrm{,}j\mathrm{<mo\; stretchy="false">)</mo>}$

$\times {\displaystyle \prod _{l\mathrm{<mo>=</mo>1}}^{b_i}}{\displaystyle \sum _r}t\mathrm{<mo\; stretchy="false">(</mo>}r\mathrm{<mo>|</mo>}k\mathrm{,}j\mathrm{<mo\; stretchy="false">)</mo>}f\mathrm{<mo\; stretchy="false">(</mo>}{y}_{il}\mathrm{<mo>|</mo>}\theta \mathrm{,}{g}_{if}\mathrm{<mo>=</mo>}k\mathrm{,}{g}_{im}\mathrm{<mo>=</mo>}j\mathrm{,}{g}_{il}\mathrm{<mo>=</mo>}r\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>,<mo\; stretchy="false">(</mo>1<mo\; stretchy="false">)</mo>}$

where typically the density is assumed multivariate normal with covariate matrix Σ, f(⋅|θ,k) is the density for residual ∈i with genotype g_i in the mean vector specification, f(y_l|θ,i,j,r) is the density for residual ∈_r with genotype g_r and with adjusted mean and variance given by

μ(g_il=r)-ΩΣ^-1[(y_if-μ(g_f=k))+(y_m-μ(g_m=j))] and Σ-ΩΣ^-1Ω

where Σ=Cov(Y_l,Y_l) and Ω= Cov(Y_l,Y_m); K(θ,I,j) is a quantity that depends on the parents’ genotypes and the mean [18]. It is well known that when the number of genotypes is relatively large, this model is computationally inefficient [19]. Proposed a computational more efficient model. In light of [19], let G_f=G_m=(π₁,…,π_I), then the joint likelihood for the i-th family is written as

$L_i\mathrm{<mo\; stretchy="false">(</mo>}{y}_i\mathrm{<mo>|</mo>}\theta \mathrm{<mo\; stretchy="false">)</mo><mo>=</mo><mo\; stretchy="false">(</mo>}{\displaystyle \sum _k}{\pi }_{k}f\mathrm{<mo\; stretchy="false">(</mo>}{y}_{if}\mathrm{<mo>|</mo>}\theta \mathrm{,}{g}_{if}\mathrm{<mo>=</mo>}k\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo><mo\; stretchy="false">(</mo>}{\displaystyle \sum _j}{\pi }_{j}f\mathrm{<mo\; stretchy="false">(</mo>}{y}_{im}\mathrm{<mo>|</mo>}\theta \mathrm{,}{g}_{im}\mathrm{<mo>=</mo>}j\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}$

$K\mathrm{<mo\; stretchy="false">(</mo>}\theta \mathrm{<mo\; stretchy="false">)</mo>}{\displaystyle \prod _{l\mathrm{<mo>=</mo>1}}^{b_i}}{\displaystyle \sum _r}T\mathrm{<mo\; stretchy="false">(</mo>}r\mathrm{<mo\; stretchy="false">)</mo>}f\mathrm{<mo\; stretchy="false">(</mo>}{y}_{il}\mathrm{<mo>|</mo>}G\mathrm{,}\theta \mathrm{,}{g}_{il}\mathrm{<mo>=</mo>}r\mathrm{<mo\; stretchy="false">)</mo>,<mo\; stretchy="false">(</mo>2<mo\; stretchy="false">)</mo>}$

where $K\mathrm{<mo\; stretchy="false">(</mo>}\theta \mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{\displaystyle {\sum }_{i\mathrm{<mo>=</mo>1}}^k}{\displaystyle {\sum }_{j\mathrm{<mo>=</mo>1}}^k}K\mathrm{<mo\; stretchy="false">(</mo>}\theta \mathrm{,}i\mathrm{,}j\mathrm{<mo\; stretchy="false">)</mo>}$

$T\mathrm{<mo\; stretchy="false">(</mo>}r\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}T\mathrm{<mo\; stretchy="false">(</mo>}r\mathrm{<mo>|</mo>}{G}_f\mathrm{,}{G}_m\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{\displaystyle \sum _{i\mathrm{<mo>=</mo>1}}^k}{\displaystyle \sum _{s\mathrm{<mo>=</mo>1}}^k}t\mathrm{<mo\; stretchy="false">(</mo>}r\mathrm{<mo>|</mo>}i\mathrm{,}s\mathrm{<mo\; stretchy="false">)</mo>}P\mathrm{<mo\; stretchy="false">(</mo>}{g}_f\mathrm{<mo>=</mo>}i\mathrm{<mo\; stretchy="false">)</mo>}P\mathrm{<mo\; stretchy="false">(</mo>}{g}_m\mathrm{<mo>=</mo>}s\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{\displaystyle \sum _{i\mathrm{<mo>=</mo>1}}^k}{\displaystyle \sum _{s\mathrm{<mo>=</mo>1}}^k}t\mathrm{<mo\; stretchy="false">(</mo>}r\mathrm{<mo>|</mo>}i\mathrm{,}s\mathrm{<mo\; stretchy="false">)</mo>}{\pi }_i{\pi }_s\mathrm{,}$

In the mean specification f(y_l|G,θ,r) is the density of residual ∈l with genotype gl=r with adjusted mean and variance given by

$\mu \mathrm{<mo\; stretchy="false">(</mo>}{g}_{il}\mathrm{<mo>=</mo>}r\mathrm{<mo\; stretchy="false">)</mo>}-{{\Omega }^{\prime }}_p{\Sigma }_p^{-1}\mathrm{<mo\; stretchy="false">(</mo>}{y}_p-{\mu }_p\mathrm{<mo\; stretchy="false">)</mo><mi>\; </mi><mi>\; </mi><mi>\; </mi><mi>\; </mi>}and\mathrm{<mi>\; </mi><mi>\; </mi><mi>\; </mi><mi>\; </mi>}\Sigma -{{\Omega }^{\prime }}_p{\Sigma }_p^{-1}{\Omega }_p\mathrm{,}$

where y_p=(y_f,y_m),μ_p=(μ_f,μ_m), ${\mu }_f\mathrm{<mo>=</mo>}{\displaystyle {\sum }_{i\mathrm{<mo>=</mo>1}}^k}{\pi }_i{\mu }_f\mathrm{<mo\; stretchy="false">(</mo>}{g}_f\mathrm{<mo>=</mo>}i\mathrm{<mo\; stretchy="false">)</mo>}$ similarly for ; μ_m

${\Omega }_p\mathrm{<mo>=</mo>}\left(\begin{array}{c}\Omega \\ \Omega \end{array}\right)\mathrm{<mi>\; </mi><mi>\; </mi><mi>\; </mi><mi>\; </mi>}{\Sigma }_p\mathrm{<mo>=</mo>}\left(\begin{array}{cc}\Sigma & \Omega \\ \Omega & \Sigma \end{array}\right)\mathrm{.}$

In comparison, model (1) has three layers of mixing (summation) corresponding to b_ik³ function evaluations that grow exponentially with the number of genotypes. On the other hand, model (2) has only one layer of mixing in three factors each, or (b_i+2)k function evaluations that are linearly proportional the number of genotypes. The reduction of computation will be more significant for multiple loci case.

Here we extend this model in the case of censoring and partial observation. In this case, the mean is modeled as

$\mu \mathrm{<mo\; stretchy="false">(</mo>}{g}_{ir}\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{I}_{ir}\odot \mathrm{<mo\; stretchy="false">(</mo>}{\mu }_0+\alpha \chi \mathrm{<mo\; stretchy="false">(</mo>}{g}_{ir}\mathrm{<mo\; stretchy="false">)</mo>}+\beta \mathrm{<mo\; stretchy="false">)</mo>}{J}_{ir}\odot x_{ir}\mathrm{,}-0.6cm\text{ (3)}$

where the operation $I_{ir}\odot {\mu }_0$ means the projection of μ₀ onto the subspace corresponding to the nonzero elements of Iir, similarly for $I_{ir}\odot \alpha \chi \mathrm{<mo\; stretchy="false">(</mo>}{g}_{ir}\mathrm{<mo\; stretchy="false">)</mo>}$ and $J_{ir}\odot x_{ir}$ . The corresponding error is now $I_{ir}\odot {\epsilon }_{ir}$ .

Recall that in the case of 1-dimensional observation without genetic implementation, the likelihood for an observation yi with a censoring indicator di is

$L_i\mathrm{<mo\; stretchy="false">(</mo>}{y}_i\mathrm{<mo>|</mo>}\theta \mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}f{\mathrm{<mo\; stretchy="false">(</mo>}{y}_i\mathrm{<mo>|</mo>}\theta \mathrm{<mo\; stretchy="false">)</mo>}}^{{\delta }_i}S{\mathrm{<mo\; stretchy="false">(</mo>}{y}_i\mathrm{<mo>|</mo>}\theta \mathrm{<mo\; stretchy="false">)</mo>}}^{1-{\delta }_i}\mathrm{.}$

To extend this to our situation, for any dimension indicator I_i, and any d-variable function v(⋅), let $I_i\odot v\mathrm{<mo\; stretchy="false">(</mo>}\cdot \mathrm{<mo\; stretchy="false">)</mo>}$ be the marginal version of v(⋅) with respect to the non-zero entry of I_i and ${\delta }_i{I}_i\odot v\mathrm{<mo\; stretchy="false">(</mo>}\cdot \mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{\delta }_i\odot \mathrm{<mo\; stretchy="false">(</mo>}{I}_i\odot v\mathrm{<mo\; stretchy="false">(</mo>}\cdot \mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>.}$ Let 1-d_i be the indicator with the same length of d_i but with 0 and 1 reversed. The full likelihood is

$L\mathrm{<mo\; stretchy="false">(</mo>}z\mathrm{<mo>|</mo>}\theta \mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{\displaystyle \prod _{i\mathrm{<mo>=</mo>1}}^n}\mathrm{<mo\; stretchy="false">(</mo>}{\displaystyle \sum _{j\mathrm{<mo>=</mo>1}}^k}{\pi }_j{\delta }_{if}{I}_{if}\odot f\mathrm{<mo\; stretchy="false">(</mo>}{y}_{if}\mathrm{<mo>|</mo>}\theta \mathrm{,}j\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo><mo\; stretchy="false">(</mo>}{\displaystyle \sum _{j\mathrm{<mo>=</mo>1}}^k}{\pi }_j\mathrm{<mo\; stretchy="false">(</mo>1}-{\delta }_{if}\mathrm{<mo\; stretchy="false">)</mo>}{I}_{if}\odot S\mathrm{<mo\; stretchy="false">(</mo>}{y}_{if}\mathrm{<mo>|</mo>}\theta \mathrm{,}j\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}$

$\times \mathrm{<mo\; stretchy="false">(</mo>}{\displaystyle \sum _{j\mathrm{<mo>=</mo>1}}^k}{\pi }_j{\delta }_{im}{I}_{im}\odot f\mathrm{<mo\; stretchy="false">(</mo>}{y}_{im}\mathrm{<mo>|</mo>}\theta \mathrm{,}j\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo><mo\; stretchy="false">(</mo>}{\displaystyle \sum _{j\mathrm{<mo>=</mo>1}}^k}{\pi }_j\mathrm{<mo\; stretchy="false">(</mo>1}-{\delta }_{im}\mathrm{<mo\; stretchy="false">)</mo>}{I}_{im}\odot S\mathrm{<mo\; stretchy="false">(</mo>}{y}_{im}\mathrm{<mo>|</mo>}\theta \mathrm{,}j\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}$

$\times {\displaystyle \prod _{j\mathrm{<mo>=</mo>1}}^{b_i}}\mathrm{<mo\; stretchy="false">(</mo>}{\displaystyle \sum _{r\mathrm{<mo>=</mo>1}}^k}T\mathrm{<mo\; stretchy="false">(</mo>}r\mathrm{<mo\; stretchy="false">)</mo>}{\delta }_{ij}{I}_{ij}\odot f\mathrm{<mo\; stretchy="false">(</mo>}{y}_j\mathrm{<mo>|</mo>}G\mathrm{,}\theta \mathrm{,}r\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo><mo\; stretchy="false">(</mo>}{\displaystyle \sum _{r\mathrm{<mo>=</mo>1}}^k}T\mathrm{<mo\; stretchy="false">(</mo>}r\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">(</mo>1}-{\delta }_{ij}\mathrm{<mo\; stretchy="false">)</mo>}{I}_{ij}\odot S\mathrm{<mo\; stretchy="false">(</mo>}{y}_j\mathrm{<mo>|</mo>}G\mathrm{,}\theta \mathrm{,}r\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>.<mo\; stretchy="false">(</mo>4<mo\; stretchy="false">)</mo>}$

Here extra caution should be taken since the observation vector from each individual may vary in dimensions and sub-spaces. For a d-dimensional vector v, let ${\delta }_{ij}{I}_{ij}\odot v$ be its margin with respect to ; and for a d-dimensional matrix A, ${\delta }_{ij}{I}_{ij}\odot A\odot I_{il}{\delta }_{il}$ denote the submatrix by the rows corresponding to the non-zero entry of d_ijI_ij and columns corresponding to the non-zero entry of d_ilI_il. In particular, ${\delta }_{ij}{I}_{ij}\odot A\odot I_{il}{\delta }_{il}$ has adjusted mean given by

${\delta }_{ij}{I}_{ij}\odot \mu \mathrm{<mo\; stretchy="false">(</mo>}{g}_{ij}\mathrm{<mo>=</mo>}r\mathrm{<mo\; stretchy="false">)</mo>}-\mathrm{<mo\; stretchy="false">(</mo>}{\delta }_{ij}{I}_{ij}\odot {{\Omega }^{\prime }}_p\odot I_{ip}{\delta }_{ip}\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">(</mo>}{\delta }_{ip}{I}_{ip}\odot {\Sigma }_p^{-1}\odot I_{ip}{\delta }_{ip}\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">(</mo>}{\delta }_{ip}{I}_{ip}\odot \mathrm{<mo\; stretchy="false">(</mo>}{y}_p-{\mu }_p\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>,}$

and adjusted variance matrix given by

${\delta }_{ij}{I}_{ij}\odot \Sigma \odot I_{ij}{\delta }_{ij}-\mathrm{<mo\; stretchy="false">(</mo>}{\delta }_{ij}{I}_{ij}\odot {{\Omega }^{\prime }}_p\odot I_{ip}{\delta }_{ip}\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">(</mo>}{\delta }_{ip}{I}_{ip}\odot {\Sigma }_p^{-1}\odot I_{ip}{\delta }_{ip}\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">(</mo>}{\delta }_{ip}{I}_{ip}\odot {\Omega }_p\odot I_{ij}{\delta }_{ij}\mathrm{<mo\; stretchy="false">)</mo>,}$

where I_ip=(I_if,I_im). The corresponding adjustment in $\mathrm{<mo\; stretchy="false">(</mo>1}-{\delta }_{ij}\mathrm{<mo\; stretchy="false">)</mo>}{I}_{ij}\odot S\mathrm{<mo\; stretchy="false">(</mo>}{y}_j\mathrm{<mo>|</mo>}G\mathrm{,}\theta \mathrm{,}r\mathrm{<mo\; stretchy="false">)</mo>}$is made. The parameter θ is estimated by its MLE $\widehat{\theta }$ under (3), along with the restriction ${\displaystyle {\sum }_{j\mathrm{<mo>=</mo>1}}^k}{\pi }_j\mathrm{<mo>=</mo>1}$

Semiparametric model

For censored data, a commonly used semi parametric regression model is Cox’s proportional hazards model [20,21]. In the univariate case, let y₍₃₎<y₍₄₎<…<y_(n) be the ordered observations of y₁,…,y_n (assume no ties for simplicity), x_(i) and d_(i) be the associated quantities, for y_(i), of the xi’s and di’s. Let R_(i) be the i-th risk set, the set of all individuals who are still under study at the ‘time’ just prior to y_(i), U be the set of all uncensored individuals, and

$\lambda \mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{<mo>|</mo>}x\mathrm{,}\theta \mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}\frac{f\mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{<mo>|</mo>}x\mathrm{,}\theta \mathrm{<mo\; stretchy="false">)</mo>}}{1-F\mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{<mo>|</mo>}x\mathrm{,}\theta \mathrm{<mo\; stretchy="false">)</mo>}}$

be the hazard function. The proportional hazards model has a form of λ(y/x,θ)=h(β'x)λ0(y)-0.2cm for some known positive function h(⋅), and unspecified baseline hazard rate λ0(⋅), which implies that the distribution belongs to the Lehmann family [4] 1-F=(1-F₀) for some F₀(⋅) and γ>₀. Under these assumptions, the conditional likelihood (partial likelihood [20,21]; marginal rank likelihood, [4]) is

$L_c\mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{<mo>|</mo>}\theta \mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{\displaystyle \prod _{i\in U}}\frac{h\mathrm{<mo\; stretchy="false">(</mo>}{\beta }^{\prime }{x}_{\mathrm{<mo\; stretchy="false">(</mo>}i\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">)</mo>}}{{\displaystyle \sum _{j\in R_{\mathrm{<mo\; stretchy="false">(</mo>}i\mathrm{<mo\; stretchy="false">)</mo>}}}}h\mathrm{<mo\; stretchy="false">(</mo>}{\beta }^{\prime }{x}_{\mathrm{<mo\; stretchy="false">(</mo>}j\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{,}$

where the estimate of θ is the MLE $\widehat{\theta }$under L_c(y|θ). The optimality property of $\widehat{\theta }$is studied extensively. In the case of multivariate observations, various extensions of this method have focused on each marginal distribution and Markov chain Monte Carlo on the margins [22]. Proposed a multivariate extension of the proportional hazards model, or frailty model, which is equivalent to an exponential specification of the joint survival function [23]. Proposed a class of multivariate failure time distributions, including a multivariate version of Cox’s proportional hazards model, in which the within family dependence is modeled by a common latent variable with a known parametric distribution given that all the family members are independent. Then the joint distribution is obtained by taking expectation of the conditional one. All these frailty models assume that there is a shared common dependent latent variable. This assumption basically requires that the distribution be interchangeable among the involved individuals. This is reasonable for some familial data but not generally true. Other existing multivariate proportional hazards models [24-26] are similar in nature. Here we model the within family dependence in a manifest way to be desirable for our genetic analysis. We adopt a successive conditional version of the proportional hazards model where we assume a special semi parametric form of the survival function in order to evaluate the conditioning in closed form easily. More specifically, in our multivariate proportional hazards model, we assume h(⋅) and λ₀(⋅) are functions of d-variates each. Let y_i,(3),…,y_i,(ni) be the ordered observations on the i-th variable (i=1,…,d), define x_i,(j),d_i,(j),R_i,(j), and U_i accordingly. Note there are structures in h(⋅) through the dependent effects among the covariates. Recall β'x is a d-vector, let

$h\mathrm{<mo\; stretchy="false">(</mo>}{\beta }^{\prime }x\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{e}^{-\frac{1}{2}\mathrm{<mo\; stretchy="false">(</mo>}{\beta }^{\prime }x{)}^{\prime }\Omega {\beta }^{\prime }x}\mathrm{,}$

where Ω is the within individual covariance matrix. Then h(⋅) behaves as a d-variate normal density, and its marginal and conditional versions are well defined and in closed forms, although it is not a proper density function. We need the successive ‘conditioning’ form of h(⋅) to apply the proportional hazards method. Specifically, let w_ij be the j-th diagonal element of Ω where Ω_j be the upperleft j-dimensional sub-matrix of it, a_j be the first j elements in the j-th column; [β'x]_j be the first j components of β'x,h_j+1|j(⋅), be the conditional version of covariates [β'x]_j+1 given [β'x]j. Then h_j+1|j(β'x) is a univariate normal kernel with mean $-{a^{\prime }}_j{\Omega }_j^{-1}{\mathrm{<mo\; stretchy="false">(</mo>}{\beta }^{\prime }x\mathrm{<mo\; stretchy="false">)</mo>}}_j$ and variance${\omega }_{jj}-{a^{\prime }}_j{\Omega }_j^{-1}a$, and

$h\mathrm{<mo\; stretchy="false">(</mo>}{\beta }^{\prime }x\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{h}_{\mathrm{1<mo>|</mo>0}}\mathrm{<mo\; stretchy="false">(</mo><mo\; stretchy="false">[</mo>}{\beta }^{\prime }x{\mathrm{<mo\; stretchy="false">]</mo>}}_1\mathrm{<mo\; stretchy="false">)</mo>}{\displaystyle \prod _{j\mathrm{<mo>=</mo>1}}^{d-1}}{h}_{j+\mathrm{1<mo>|</mo>}j}\mathrm{<mo\; stretchy="false">(</mo><mo\; stretchy="false">[</mo>}{\beta }^{\prime }x{\mathrm{<mo\; stretchy="false">]</mo>}}_{j+1}\mathrm{<mo\; stretchy="false">)</mo>.<mo\; stretchy="false">(</mo>5<mo\; stretchy="false">)</mo>}$

Thus, without mixing over gene, for singleton multivariate observations, the joint conditional likelihood is

$L_c\mathrm{<mo\; stretchy="false">(</mo>}z\mathrm{<mo>|</mo>}\theta \mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{\displaystyle \prod _{i\in U_1}}\frac{h\mathrm{<mo\; stretchy="false">(</mo><mo\; stretchy="false">[</mo>}{\beta }^{\prime }{x}_{\mathrm{<mo\; stretchy="false">(</mo>}i\mathrm{<mo\; stretchy="false">)</mo>}}{\mathrm{<mo\; stretchy="false">]</mo>}}_1\mathrm{<mo\; stretchy="false">)</mo>}}{{\displaystyle \sum _{l\in R_{\mathrm{1,<mo\; stretchy="false">(</mo>}i\mathrm{<mo\; stretchy="false">)</mo>}}}}h\mathrm{<mo\; stretchy="false">(</mo><mo\; stretchy="false">[</mo>}{\beta }^{\prime }{x}_{\mathrm{<mo\; stretchy="false">(</mo>}l\mathrm{<mo\; stretchy="false">)</mo>}}{\mathrm{<mo\; stretchy="false">]</mo>}}_1\mathrm{<mo\; stretchy="false">)</mo>}}{\displaystyle \prod _{j\mathrm{<mo>=</mo>1}}^{d-1}}{\displaystyle \prod _{i\in U_j}}\frac{h_{j+\mathrm{1<mo>|</mo>}j}\mathrm{<mo\; stretchy="false">(</mo><mo\; stretchy="false">[</mo>}{\beta }^{\prime }{x}_{\mathrm{<mo\; stretchy="false">(</mo>}i\mathrm{<mo\; stretchy="false">)</mo>}}{\mathrm{<mo\; stretchy="false">]</mo>}}_j\mathrm{<mo\; stretchy="false">)</mo>}}{{\displaystyle \sum _{l\in R_{j\mathrm{,<mo\; stretchy="false">(</mo>}i\mathrm{<mo\; stretchy="false">)</mo>}}}}{h}_{j+\mathrm{1<mo>|</mo>}j}\mathrm{<mo\; stretchy="false">(</mo><mo\; stretchy="false">[</mo>}{\beta }^{\prime }{x}_{\mathrm{<mo\; stretchy="false">(</mo>}l\mathrm{<mo\; stretchy="false">)</mo>}}{\mathrm{<mo\; stretchy="false">]</mo>}}_j\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{.}$

Now for the case of nuclear family, inspired by (4), we assume h(⋅) has the form

$h\mathrm{<mo\; stretchy="false">(</mo>}{\mu }_i\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{\displaystyle \sum _{j\mathrm{<mo>=</mo>1}}^k}{\pi }_j{h}_{I_{if}}\mathrm{<mo\; stretchy="false">(</mo>}\mu \mathrm{<mo\; stretchy="false">(</mo>}{g}_j\mathrm{,}{I}_{if}\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}{\displaystyle \sum _{j\mathrm{<mo>=</mo>1}}^k}{\pi }_j{h}_{I_{im}}\mathrm{<mo\; stretchy="false">(</mo>}\mu \mathrm{<mo\; stretchy="false">(</mo>}{g}_j\mathrm{,}{I}_{im}\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}{\displaystyle \prod _{l\mathrm{<mo>=</mo>1}}^{b_i}}{\displaystyle \sum _{j\mathrm{<mo>=</mo>1}}^k}T\mathrm{<mo\; stretchy="false">(</mo>}j\mathrm{<mo\; stretchy="false">)</mo>}{h}_{I_{ij}}\mathrm{<mo\; stretchy="false">(</mo>}\mu \mathrm{<mo\; stretchy="false">(</mo>}{g}_j\mathrm{,}{I}_{ij}\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>,<mo\; stretchy="false">(</mo>6<mo\; stretchy="false">)</mo>}$

Treat h(⋅) as a ‘density’. Recall μ_i=(μif,μ_im,μ_i1,…,μ_ibi). The conditioning [μ_i]_j+1|[μ_i]_j can be applied component-wise, i.e.

[μ_i]_j+1|[μ_i]j=([μ_if]_j+1|[μ_if]_j, [μ_im]_j+1|[μ_im]_j, [μ_i1]_j+1|[μ_i1]_j,…, [μ_ibi]_j+1|[μ_ibi]_j)

Now we have

$h_{j+\mathrm{1<mo>|</mo>}j}\mathrm{<mo\; stretchy="false">(</mo><mo\; stretchy="false">[</mo>}{\mu }_i{\mathrm{<mo\; stretchy="false">]</mo>}}_{j+1}\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo><mo\; stretchy="false">(</mo>}{\displaystyle \sum _{j\mathrm{<mo>=</mo>1}}^k}{\pi }_j{h}_{I_{if}\mathrm{,1<mo>|</mo>0}}\mathrm{<mo\; stretchy="false">(</mo><mo\; stretchy="false">[</mo>}\mu \mathrm{<mo\; stretchy="false">(</mo>}{g}_j\mathrm{,}{I}_{if}{\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">]</mo>}}_1\mathrm{<mo\; stretchy="false">)</mo>}{\displaystyle \prod _{j\mathrm{<mo>=</mo>1}}^{\mathrm{<mo>|</mo>}{I}_{if}\mathrm{<mo>|</mo>}-1}}{h}_{I_{if}\mathrm{,}j+\mathrm{1<mo>|</mo>}j}\mathrm{<mo\; stretchy="false">(</mo><mo\; stretchy="false">[</mo>}\mu \mathrm{<mo\; stretchy="false">(</mo>}{g}_j\mathrm{,}{I}_{if}{\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">]</mo>}}_j\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}$

$\times \mathrm{<mo\; stretchy="false">(</mo>}{\displaystyle \sum _{j\mathrm{<mo>=</mo>1}}^k}{\pi }_j{h}_{I_{im}\mathrm{,1<mo>|</mo>0}}\mathrm{<mo\; stretchy="false">(</mo><mo\; stretchy="false">[</mo>}\mu \mathrm{<mo\; stretchy="false">(</mo>}{g}_j\mathrm{,}{I}_{im}{\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">]</mo>}}_1\mathrm{<mo\; stretchy="false">)</mo>}{\displaystyle \prod _{j\mathrm{<mo>=</mo>1}}^{\mathrm{<mo>|</mo>}{I}_{im}\mathrm{<mo>|</mo>}-1}}{h}_{I_{if}\mathrm{,}j+\mathrm{1<mo>|</mo>}j}\mathrm{<mo\; stretchy="false">(</mo><mo\; stretchy="false">[</mo>}\mu \mathrm{<mo\; stretchy="false">(</mo>}{g}_j\mathrm{,}{I}_{im}{\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">]</mo>}}_j\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}$

$\times {\displaystyle \prod _{l\mathrm{<mo>=</mo>1}}^{b_i}}\mathrm{<mo\; stretchy="false">(</mo>}{\displaystyle \sum _{j\mathrm{<mo>=</mo>1}}^k}T\mathrm{<mo\; stretchy="false">(</mo>}j\mathrm{<mo\; stretchy="false">)</mo>}{h}_{I_{il}\mathrm{,1<mo>|</mo>0}}\mathrm{<mo\; stretchy="false">(</mo><mo\; stretchy="false">[</mo>}\mu \mathrm{<mo\; stretchy="false">(</mo>}{g}_j\mathrm{,}{I}_{il}{\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">]</mo>}}_1\mathrm{<mo\; stretchy="false">)</mo>}{\displaystyle \prod _{j\mathrm{<mo>=</mo>1}}^{\mathrm{<mo>|</mo>}{I}_{il}\mathrm{<mo>|</mo>}-1}}{h}_{I_{il}\mathrm{,}j+\mathrm{1<mo>|</mo>}j}\mathrm{<mo\; stretchy="false">(</mo><mo\; stretchy="false">[</mo>}\mu \mathrm{<mo\; stretchy="false">(</mo>}{g}_j\mathrm{,}{I}_{il}{\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">]</mo>}}_j\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>.<mo\; stretchy="false">(</mo>7<mo\; stretchy="false">)</mo>}$

In (7), [μ(g_j,I_ir)]1, means the first component of μ(g_j,I_ir) in I_ir, and |I_ir| denote its cardinality (r=f,m,1,…,b_i). Now, the conditional likelihood is

$L_c\mathrm{<mo\; stretchy="false">(</mo>}z\mathrm{<mo>|</mo>}\theta \mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{\displaystyle \prod _{i\in U_1}}\frac{h\mathrm{<mo\; stretchy="false">(</mo><mo\; stretchy="false">[</mo>}{\mu }_{\mathrm{<mo\; stretchy="false">(</mo>}i\mathrm{<mo\; stretchy="false">)</mo>}}{\mathrm{<mo\; stretchy="false">]</mo>}}_1\mathrm{<mo\; stretchy="false">)</mo>}}{{\displaystyle \sum _{l\in R_{\mathrm{1,<mo\; stretchy="false">(</mo>}i\mathrm{<mo\; stretchy="false">)</mo>}}}}h\mathrm{<mo\; stretchy="false">(</mo><mo\; stretchy="false">[</mo>}{\mu }_{\mathrm{<mo\; stretchy="false">(</mo>}i\mathrm{<mo\; stretchy="false">)</mo>}}{\mathrm{<mo\; stretchy="false">]</mo>}}_1\mathrm{<mo\; stretchy="false">)</mo>}}{\displaystyle \prod _{j\mathrm{<mo>=</mo>1}}^{d-1}}{\displaystyle \prod _{i\in U_j}}\frac{h_{j+\mathrm{1<mo>|</mo>}j}\mathrm{<mo\; stretchy="false">(</mo><mo\; stretchy="false">[</mo>}{\mu }_{\mathrm{<mo\; stretchy="false">(</mo>}i\mathrm{<mo\; stretchy="false">)</mo>}}{\mathrm{<mo\; stretchy="false">]</mo>}}_j\mathrm{<mo\; stretchy="false">)</mo>}}{{\displaystyle \sum _{l\in R_{j\mathrm{,<mo\; stretchy="false">(</mo>}i\mathrm{<mo\; stretchy="false">)</mo>}}}}{h}_{j+\mathrm{1<mo>|</mo>}j}\mathrm{<mo\; stretchy="false">(</mo><mo\; stretchy="false">[</mo>}{\mu }_{\mathrm{<mo\; stretchy="false">(</mo>}i\mathrm{<mo\; stretchy="false">)</mo>}}{\mathrm{<mo\; stretchy="false">]</mo>}}_j\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{,<mo\; stretchy="false">(</mo>8<mo\; stretchy="false">)</mo>}$

where h_j+1|j([μ_(i)]j is given by (7). The MLE $\widehat{\theta }$ of θ is obtained under (8).

Nonparametric model

For univariate censored data, [27,28] considered a class of estimators, including the weighted least squares estimators, for censored data. Here the weights are determined by the ordered statistics of the observations and the associated censoring indicators, and are derived from the empirical survival function, i.e., the Kaplan- Meier product limit estimator [29-32]. Formulated the multivariate Kaplan-Meier estimator. Using the product integral, the mathematical expressions are quite involved. So instead of choosing the weights according to the multivariate Kaplan-Meier estimator, we use the nonparametric locally weighted least squares method, also called locally linear regression smoothers [33,34]. Let Y and X be the d and J-dimensional random vectors corresponding to the full observation and the covariates for an individual. Let μ(x)=E(Y|X=x) denote the regression function. In the univariate observation case, the locally linear estimator $\widehat{\mu }\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo\; stretchy="false">)</mo>}$ of μ(x) is first to find $\widehat{a}$ and $\widehat{b}$ to minimize

${\displaystyle \sum _{i\mathrm{<mo>=</mo>1}}^n}{\mathrm{<mo\; stretchy="false">(</mo>}{y}_i-a-b\mathrm{<mo\; stretchy="false">(</mo>}x-x_i\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}}^{2}K\mathrm{<mo\; stretchy="false">(</mo>}\frac{x-x_i}{h_n}\mathrm{<mo\; stretchy="false">)</mo>,}$

Where K(⋅) is a kernel function, hn is the bandwidth, and $\widehat{\mu }\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}\widehat{a}$ . In our case, keep the notations in section 1. We choose the kernel to be the J-dimensional standard normal density f(⋅), and (⋅) to be its distribution function. To simplify the expression of the likelihood, let

${\tilde{y}}_{ifr}\mathrm{<mo>=</mo>}{y}_{if}-I_{if}\odot \mu \mathrm{<mo\; stretchy="false">(</mo>}r\mathrm{<mo\; stretchy="false">)</mo>,}{\tilde{x}}_{if}\mathrm{<mo>=</mo>}{J}_{if}\odot \mathrm{<mo\; stretchy="false">(</mo>}x-x_{if}\mathrm{<mo\; stretchy="false">)</mo>,}{\stackrel{⌣}{x}}_{if}\mathrm{<mo>=</mo><mo\; stretchy="false">(</mo>1}-J_{if}\mathrm{<mo\; stretchy="false">)</mo>}\odot \mathrm{<mo\; stretchy="false">(</mo>}x-x_{if}\mathrm{<mo\; stretchy="false">)</mo>,}$

and similarly for imr ${\tilde{y}}_{imr}$, ${\tilde{y}}_{ijr}$, ${\tilde{x}}_{im}$, ${\stackrel{⌣}{x}}_{im}$ , ${\tilde{x}}_{ij}$, and ${\tilde{x}}_{ij}$. To estimate$\widehat{\mu }\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo\; stretchy="false">)</mo>}$ , inspired from the univariate locally linear estimator and (4), we first find $\mathrm{<mo\; stretchy="false">(</mo>}{\widehat{\mu }}_0\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo\; stretchy="false">)</mo>,}\widehat{\alpha }\mathrm{,}\widehat{\beta }\mathrm{,}\widehat{\pi }\mathrm{<mo\; stretchy="false">)</mo>}$ to minimize

${\displaystyle \sum _{i\mathrm{<mo>=</mo>1}}^n}\mathrm{<mo\; stretchy="false">(</mo>}{\displaystyle \sum _{r\mathrm{<mo>=</mo>1}}^k}{\pi }_r^2\tilde{y}{\mathrm{<mo>\text{'}</mo>}}_{ifr}\mathrm{<mo\; stretchy="false">(</mo>}{I}_{if}\odot \Omega \odot I_{if}\mathrm{<mo\; stretchy="false">)</mo>}{\tilde{y}}_{ifr}{J}_{if}\odot \phi \mathrm{<mo\; stretchy="false">(</mo>}\frac{{\tilde{x}}_{if}}{h_n}\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">(</mo>1}-\mathrm{<mo\; stretchy="false">(</mo>1}-J_{if}\mathrm{<mo\; stretchy="false">)</mo>}\odot \Phi \mathrm{<mo\; stretchy="false">(</mo>}\frac{{\stackrel{⌣}{x}}_{if}}{h_n}\mathrm{<mo\; stretchy="false">)</mo>}$

$+{\displaystyle \sum _{r\mathrm{<mo>=</mo>1}}^k}{\pi }_r^2\tilde{y}{\mathrm{<mo>\text{'}</mo>}}_{imr}\mathrm{<mo\; stretchy="false">(</mo>}{I}_{im}\odot \Omega \odot I_{im}\mathrm{<mo\; stretchy="false">)</mo>}{\tilde{y}}_{imr}{J}_{im}\odot \phi \mathrm{<mo\; stretchy="false">(</mo>}\frac{{\tilde{x}}_{im}}{h_n}\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">(</mo>1}-\mathrm{<mo\; stretchy="false">(</mo>1}-J_{im}\mathrm{<mo\; stretchy="false">)</mo>}\odot \Phi \mathrm{<mo\; stretchy="false">(</mo>}\frac{{\stackrel{⌣}{x}}_{im}}{h_n}\mathrm{<mo\; stretchy="false">)</mo>}$

$+{\displaystyle \sum _{j\mathrm{<mo>=</mo>1}}^{b_i}}{\displaystyle \sum _{r\mathrm{<mo>=</mo>1}}^k}{T}^2\mathrm{<mo\; stretchy="false">(</mo>}r\mathrm{<mo\; stretchy="false">)</mo>}\tilde{y}{\mathrm{<mo>\text{'}</mo>}}_{ijr}\mathrm{<mo\; stretchy="false">(</mo>}{I}_{ij}\odot \Omega \odot I_{ij}\mathrm{<mo\; stretchy="false">)</mo>}{y}_{ijr}{J}_{ij}\odot \phi \mathrm{<mo\; stretchy="false">(</mo>}\frac{{\tilde{x}}_{ij}}{h_n}\mathrm{<mo>|</mo>}{\tilde{x}}_p\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">(</mo>1}-\mathrm{<mo\; stretchy="false">(</mo>1}-J_{ij}\mathrm{<mo\; stretchy="false">)</mo>}\odot \Phi \mathrm{<mo\; stretchy="false">(</mo>}\frac{{\stackrel{⌣}{x}}_{ij}}{h_n}\mathrm{<mo>|</mo>}{\stackrel{⌣}{x}}_p\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}\text{ (9)}$

Where Ω is the within individual variance matrix, ($\phi \mathrm{<mo\; stretchy="false">(</mo>}\cdot \mathrm{<mo>|</mo>}{\tilde{x}}_p\mathrm{<mo\; stretchy="false">)</mo>}$ and $\Phi \mathrm{<mo\; stretchy="false">(</mo>}\cdot \mathrm{<mo>|</mo>}{\stackrel{⌣}{x}}_p\mathrm{<mo\; stretchy="false">)</mo>}$are the adjusted quantities as those in (4). And $I_{ir}\odot \Omega \odot I_{ir}$I is the sub-matrix of Ω with rows and columns corresponding to the non-zero elements of I_im.

We estimate

We estimate Ω by$\widehat{\Omega }\mathrm{<mo>=</mo><mo\; stretchy="false">(</mo>}{\widehat{\omega }}_{rs}\mathrm{<mo\; stretchy="false">)</mo>}$ with

${\widehat{\omega }}_{rs}\mathrm{<mo>=</mo>}\frac{1}{n_{rs}-1}{\displaystyle \sum _{l\mathrm{<mo>=</mo>1}}^{n_{rs}}}\mathrm{<mo\; stretchy="false">(</mo>}{z}_{rl}-{\overline{z}}_l\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">(</mo>}{z}_{sl}-{\overline{z}}_s\mathrm{<mo\; stretchy="false">)</mo>,<mi>\; </mi><mi>\; </mi><mi>\; </mi><mi>\; </mi>}i\mathrm{,}j\mathrm{<mo>=</mo>1,...,}d\mathrm{,}$

where nrs is the total number of individuals with non-missing (r,s)- th components, zrls are the rearrangement of the r-th component of $y_{it}-{\displaystyle \sum _{j=1}^k{\pi }_k}{I}_{it}\odot \mu (j)$ for which the (r,s)-th components are non-missing.

Let $\widehat{\mu }\mathrm{<mo\; stretchy="false">(</mo>}r\mathrm{,}x\mathrm{<mo\; stretchy="false">)</mo>,}\widehat{\pi }\mathrm{<mo\; stretchy="false">)</mo>}$ be the minimize of (9), where the full $\widehat{\mu }\mathrm{<mo\; stretchy="false">(</mo>}r\mathrm{,}x\mathrm{<mo\; stretchy="false">)</mo>}$ depends on the genotype r and the point value x. It has the intercept term ${\widehat{\mu }}_0\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo\; stretchy="false">)</mo>}$ (recall (3)), and $\widehat{\mu }\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo\; stretchy="false">)</mo>}$ is approximated by setting $\widehat{\mu }\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{\widehat{\mu }}_0\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo\; stretchy="false">)</mo>}$ Direct computation of $\widehat{\mu }\mathrm{<mo\; stretchy="false">(</mo>}r\mathrm{,}x\mathrm{<mo\; stretchy="false">)</mo>,}\widehat{\pi }\mathrm{<mo\; stretchy="false">)</mo>}$ in (7) is not easy, instead we use an iterative procedure as in the following steps.

Select starting values π⁽⁰⁾ for π. With this π⁽⁰⁾, compute Ω⁽⁰⁾, and T⁽⁰⁾(r)s. Let η=(μ₀,α,β) be the full representation of the regression parameters, X_ir(r=f,m,1,…,n_i) be the corresponding design matrix for the r-th individual in the i-th family. In iterations 1-m do the following

(i) Fix π(i) , Ω(i) and T(i)(r) s minimize (7) with respect to η to get

${\eta }^{\mathrm{<mo\; stretchy="false">(</mo>}i\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo>=</mo><mo\; stretchy="false">(</mo>}{\displaystyle \sum _{i\mathrm{<mo>=</mo>1}}^n}{\displaystyle \sum _{r\mathrm{<mo>=</mo>1}}^k}{\displaystyle \sum _l}{\varpi }_{lr}^{\mathrm{<mo\; stretchy="false">(</mo>}i\mathrm{<mo\; stretchy="false">)</mo>2}}{X}_{il}{X^{\prime }}_{il}{J}_{il}\odot \phi \mathrm{<mo\; stretchy="false">(</mo>}\frac{{\tilde{x}}_{il}}{h_n}\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">(</mo>1}-\mathrm{<mo\; stretchy="false">(</mo>1}-J_{il}\mathrm{<mo\; stretchy="false">)</mo>}\odot \Phi \mathrm{<mo\; stretchy="false">(</mo>}\frac{{\stackrel{⌣}{x}}_{il}}{h_n}{\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}}^{-1}$

$\times {\displaystyle \sum _{i\mathrm{<mo>=</mo>1}}^n}{\displaystyle \sum _{r\mathrm{<mo>=</mo>1}}^k}{\displaystyle \sum _l}{\varpi }_{lr}^{\mathrm{<mo\; stretchy="false">(</mo>}i\mathrm{<mo\; stretchy="false">)</mo>2}}{y}_{il}{X^{\prime }}_{il}{J}_{il}\odot \phi \mathrm{<mo\; stretchy="false">(</mo>}\frac{{\tilde{x}}_{il}}{h_n}\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">(</mo>1}-\mathrm{<mo\; stretchy="false">(</mo>1}-J_{il}\mathrm{<mo\; stretchy="false">)</mo>}\odot \Phi \mathrm{<mo\; stretchy="false">(</mo>}\frac{{\stackrel{⌣}{x}}_{il}}{h_n}\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>,<mo\; stretchy="false">(</mo>10<mo\; stretchy="false">)</mo>}$

where ${\varpi }_{lr}^{\mathrm{<mo\; stretchy="false">(</mo>}i\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo>=</mo>}{\pi }_r^{\mathrm{<mo\; stretchy="false">(</mo>}i\mathrm{<mo\; stretchy="false">)</mo>}}$ for l=f,m, and ${\varpi }_{lr}^{\mathrm{<mo\; stretchy="false">(</mo>}i\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo>=</mo>}{T}^{\mathrm{<mo\; stretchy="false">(</mo>}i\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">(</mo>}r\mathrm{<mo\; stretchy="false">)</mo>}$ for l=1,…,n_i. (ii) Fix _(i), minimize (7) with respect to π, with the constraint ${\displaystyle {\sum }_{j\mathrm{<mo>=</mo>1}}^k}{\pi }_j\mathrm{<mo>=</mo>1}$ to get

${\pi }_r^{\mathrm{<mo\; stretchy="false">(</mo>}i+\mathrm{1<mo\; stretchy="false">)</mo>}}\mathrm{<mo>=</mo>}\frac{{\displaystyle \sum _{i\mathrm{<mo>=</mo>1}}^n}{\displaystyle \sum _{l\mathrm{<mo>=</mo>}f\mathrm{,}m}}\tilde{y}{\mathrm{<mo>\text{'}</mo>}}_{ilr}\mathrm{<mo\; stretchy="false">(</mo>}{I}_{il}\odot \Omega \odot I_{il}\mathrm{<mo\; stretchy="false">)</mo>}{\tilde{y}}_{ilr}{J}_{il}\odot \phi \mathrm{<mo\; stretchy="false">(</mo>}\frac{{\tilde{x}}_{il}}{h_n}\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">(</mo>1}-\mathrm{<mo\; stretchy="false">(</mo>1}-J_{il}\mathrm{<mo\; stretchy="false">)</mo>}\odot \Phi \mathrm{<mo\; stretchy="false">(</mo>}\frac{{\stackrel{⌣}{x}}_{il}}{h_n}\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}}{{\displaystyle \sum _{i\mathrm{<mo>=</mo>1}}^n}{\displaystyle \sum _{l\mathrm{<mo>=</mo>}f\mathrm{,}m}}{\displaystyle \sum _{r\mathrm{<mo>=</mo>1}}^k}\tilde{y}{\mathrm{<mo>\text{'}</mo>}}_{ilr}\mathrm{<mo\; stretchy="false">(</mo>}{I}_{il}\odot \Omega \odot I_{il}\mathrm{<mo\; stretchy="false">)</mo>}{\tilde{y}}_{ilr}{J}_{il}\odot \phi \mathrm{<mo\; stretchy="false">(</mo>}\frac{{\tilde{x}}_{il}}{h_n}\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">(</mo>1}-\mathrm{<mo\; stretchy="false">(</mo>1}-J_{il}\mathrm{<mo\; stretchy="false">)</mo>}\odot \Phi \mathrm{<mo\; stretchy="false">(</mo>}\frac{{\stackrel{⌣}{x}}_{il}}{h_n}\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}}\mathrm{,}\text{ (11)}$

(r=1,…,k) (11)

and update Ω⁽ⁱ⁺¹⁾, and T^(i+1(r) with ⁽ⁱ⁺¹⁾.For some pre-specified ∈>0, when the relative errors

$\frac{\mathrm{<mo>|</mo><mo\; stretchy="false">(</mo>}{\mu }^{\mathrm{<mo\; stretchy="false">(</mo>}m\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">(</mo>}r\mathrm{,}x\mathrm{<mo\; stretchy="false">)</mo>,}{\pi }^{\mathrm{<mo\; stretchy="false">(</mo>}m\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">)</mo>}-\mathrm{<mo\; stretchy="false">(</mo>}{\mu }^{\mathrm{<mo\; stretchy="false">(</mo>}m-\mathrm{1<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">(</mo>}r\mathrm{,}x\mathrm{<mo\; stretchy="false">)</mo>,}{\pi }^{\mathrm{<mo\; stretchy="false">(</mo>}m-\mathrm{1<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">)</mo><mo>|</mo>}}{\mathrm{<mo>|</mo><mo\; stretchy="false">(</mo>}{\mu }^{\mathrm{<mo\; stretchy="false">(</mo>}m-\mathrm{1<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">(</mo>}r\mathrm{,}x\mathrm{<mo\; stretchy="false">)</mo>,}{\pi }^{\mathrm{<mo\; stretchy="false">(</mo>}m-\mathrm{1<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">)</mo><mo>|</mo>}}\le \epsilon$

we stop the process at the last step m, and take

$\mathrm{<mo\; stretchy="false">(</mo>}\widehat{\mu }\mathrm{<mo\; stretchy="false">(</mo>}r\mathrm{,}x\mathrm{<mo\; stretchy="false">)</mo>,}\widehat{\pi }\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo><mo\; stretchy="false">(</mo>}{\mu }^{\mathrm{<mo\; stretchy="false">(</mo>}m\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">(</mo>}r\mathrm{,}x\mathrm{<mo\; stretchy="false">)</mo>,}{\pi }^{\mathrm{<mo\; stretchy="false">(</mo>}m\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">)</mo>}$

For arbitrary kernel and reasonably chosen band width h_n, various asymptotic results are established in case of standard nonmixture data. We conjecture that similar results will hold under some regularity conditions.

Lastly, the bandwidth determines the smoothness of the estimate. Interesting research that addresses the crucial problem of bandwidth selection can be found in [35]. There are considerable literatures for automatic methods that attempt to minimize a lack-of-fit criterion such as an integrated squared error. But most of the methods provide an optimal hn determined by some unknown quantities. For simplicity, let k=|J_ij| be the dimension of the observed covariate of the j-th (j=f,m,1,…,n_i) individual in the i-th family, for the corresponding kernel, we choose h_n=Cn^-1/(k+1), for some constant C>0, and C can be selected through numerical trial.

Variance components model

As an alternative to the mixture models considered above, the Variance Components (VC) model [36,37] has received much attention recently due to very efficient in computation as well as relatively robustness to model misspecification [38-46].

Let y_i be the trait vector of the i-th individual in the family, in case without censoring and missing records, the commonly used VC model describing the trait value is

y_i=μ+g_i+G_i+ηx_i+e_i

Where μ is the overall mean, gi is the unobserved random vector of major gene effects at the trait locus with alleles A and B, G_i is the unobserved polygenic effects vector, the ηj’s are effects associated with the covariates x_ij’s, and ei is the residual random error vector. The usual assumption is that g_i,G_i and e_i are uncorrelated and E(g_i)=E(G_i)=E(e_i)=0. When missing records are present, the model is modified as

$y_i\mathrm{<mo>=</mo>}{I}_i\odot \mathrm{<mo\; stretchy="false">(</mo>}\mu +g_i+G_i+\eta \mathrm{<mo\; stretchy="false">)</mo>}{J}_i\odot x_i+I_i\odot e_i\left(12\right)$

In this model, the parameters of interests are specified in the family variance matrix, thus computation can be carried out efficiently without the multiple mixing. Let y_k,π_k and Ω_k be the observation, its mean and variance matrix of the k-th family. We can define I_k, d_k and J_k accordingly. The commonly used model for quantitative traits is the multivariate normal distribution, thus the total likelihood is

$L\mathrm{<mo\; stretchy="false">(</mo>}z\mathrm{<mo>|</mo>}\theta \mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{\displaystyle \sum _{k\mathrm{<mo>=</mo>1}}^K}{\delta }_k{I}_k\odot \phi \mathrm{<mo\; stretchy="false">(</mo>}{y}_k-{\mu }_k\mathrm{<mo>|</mo>}{\Omega }_k\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">(</mo>1}-{\delta }_k\mathrm{<mo\; stretchy="false">)</mo>}{I}_k\odot \Phi y_k-{\mu }_k\mathrm{<mo>|</mo>}{\Omega }_k\mathrm{<mo\; stretchy="false">)</mo>.}$

Here f is the distribution function of the normal distribution with mean 0 and variance Ω.

The key lies in the specification of the variance matrices Ω_ks, which we illustrate in the following settings.

In the simplest case of Hardy-Weinberg equilibrium among locus alleles without linkage to marker, and without censoring and missing records, the covariance matrix between individuals i and j of a given family can be found, for example, in [38]. Modified to our case, it is

$Cov\mathrm{<mo\; stretchy="false">(</mo>}{Y}_{ki}\mathrm{,}{Y}_{kj}\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}\{\begin{array}{ll}{I}_i\odot \mathrm{<mo\; stretchy="false">(</mo>}{\sigma }_a^2+{\sigma }_d^2+{\sigma }_G^2+{\sigma }_e^2\mathrm{<mo\; stretchy="false">)</mo>}\odot I_j\hfill & if\mathrm{<mi>\; </mi><mi>\; </mi>}i\mathrm{<mo>=</mo>}j\hfill \\ 2{\Phi }_{ij}{I}_i\odot \mathrm{<mo\; stretchy="false">(</mo>}{\sigma }_a^2+{\Delta }_{7ij}{\sigma }_d^2+2{\Phi }_{ij}{\sigma }_G^2\mathrm{<mo\; stretchy="false">)</mo>}\odot I_j\mathrm{,}\hfill & if\mathrm{<mi>\; </mi><mi>\; </mi>}i\ne j\hfill \end{array}\left(13\right)$

where ${\sigma }_a^2$ is the additive genetic variance matrix due to the locus, ${\sigma }_d^2$ is the dominant genetic variance matrix, ${\Phi }_{ij}\mathrm{<mo>=</mo>}{\Delta }_{7ij}\mathrm{<mo>/</mo>2}+{\Delta }_{8ij}\mathrm{<mo>/</mo>4}$is the kinship coefficient between individuals i and j [47], and Δ_7ij,Δ_8ij,Δ_9ij, etc. are the condensed kinship coefficient of Jacquard [48], between individuals i and j.

In the more general Hardy-Weinberg disequilibrium case, let f be the within population inbreeding coefficient f at the trait locus [49- 51]. Introduced CV model in this case, which modified in our case is

$Cov\mathrm{<mo\; stretchy="false">(</mo>}{Y}_i\mathrm{,}{Y}_j\mathrm{<mo>|</mo>}f\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}\{\begin{array}{ll}{I}_i\odot \mathrm{<mo\; stretchy="false">(</mo><mo\; stretchy="false">(</mo>1}+\frac{f}{2}\mathrm{<mo\; stretchy="false">)</mo>}{\sigma }_a^2+\mathrm{<mo\; stretchy="false">(</mo>1}-f\mathrm{<mo\; stretchy="false">)</mo>}{\sigma }_d^2+f{\sigma }_0^2+{\sigma }_G^2+{\sigma }_e^2\mathrm{<mo\; stretchy="false">)</mo>}\odot I_j\mathrm{,}\hfill & if\mathrm{<mi>\; </mi><mi>\; </mi>}i\mathrm{<mo>=</mo>}j\hfill \\ I_i\odot \mathrm{<mo\; stretchy="false">(</mo>}{\Delta }_{7ij}{\gamma }_7\mathrm{<mo\; stretchy="false">(</mo>}f\mathrm{<mo\; stretchy="false">)</mo>}+{\Delta }_{8ij}{\gamma }_8\mathrm{<mo\; stretchy="false">(</mo>}f\mathrm{<mo\; stretchy="false">)</mo>}+2{\Phi }_{ij}{\sigma }_G^2\mathrm{<mo\; stretchy="false">)</mo>}\odot I_j\mathrm{,}\hfill & if\mathrm{<mi>\; </mi><mi>\; </mi>}i\ne j\hfill \end{array}\text{ (14)}$

where γ_l(f)s are matrices determined by 2 ${\sigma }_a^2$ , ${\sigma }_d^2$ and f etc., see there for details.

In the case of linkage to marker with both Hardy-Weinberg and linkage equilibrium, the covariance in our case can be specified based on that of, for example [40], in the same way as above. In the case of linkage to marker with either one or both Hardy-Weinberg and linkage disequilibrium, the covariance in our case can be specified based on that of [51], in the same way.

Competing risks

Now suppose that the response y_i is the failure time and only the failure for one of the d diseases is observed for each individual. For the i-th family, the data have the form (y_i,d_i,j_i,x_i), where y_i=(y_if,y_im,y_i1,… ,y_ib_i), similarly for d_i,j_i, and x_i, where j_i is the observed disease type indicator. For example if the observed disease for the father is type 2, then j_if=2. Given the data (y_i,d_i,j_i,x_i)s, we like to investigate the objective of interests for each of the d disease. This problem is that of the competing risks. Note here the response for each individual is one-dimensional, and hence the corresponding quantities have simple notations. We are interested in the genetic regression analysis for the competing risks. We use a variant of the proportional hazards model. The mean of the j-th type, r-th member of the i-th family is specified as

${\mu }_j\mathrm{<mo\; stretchy="false">(</mo>}{g}_{ir}\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{I}_{ir}\odot \mathrm{<mo\; stretchy="false">(</mo>}{\mu }_{0j}+{\alpha }_j\chi \mathrm{<mo\; stretchy="false">(</mo>}{g}_{ir}\mathrm{<mo\; stretchy="false">)</mo>}+{\beta }_i\mathrm{<mo\; stretchy="false">)</mo>}{J}_{ir}\odot x_{ir}\mathrm{,<mi>\; </mi><mi>\; </mi><mi>\; </mi><mi>\; </mi><mo\; stretchy="false">(</mo>}r\mathrm{<mo>=</mo>}f\mathrm{,}m\mathrm{,1,...,}{b}_i\mathrm{<mo\; stretchy="false">)</mo>.}$

For a reasonably chosen function h(⋅), we specify

$h\mathrm{<mo\; stretchy="false">(</mo>}{\mu }_i\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo><mo\; stretchy="false">(</mo>}{\displaystyle \sum _{l\mathrm{<mo>=</mo>1}}^k}{\pi }_{l}h\mathrm{<mo\; stretchy="false">(</mo>}{\mu }_{j_f}\mathrm{<mo\; stretchy="false">(</mo>}{g}_l\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo><mo\; stretchy="false">(</mo>}{\displaystyle \sum _{l\mathrm{<mo>=</mo>1}}^k}{\pi }_{l}h\mathrm{<mo\; stretchy="false">(</mo>}{\mu }_{j_m}\mathrm{<mo\; stretchy="false">(</mo>}{g}_l\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}{\displaystyle \prod _{l\mathrm{<mo>=</mo>1}}^{b_i}}\mathrm{<mo\; stretchy="false">(</mo>}{\displaystyle \sum _{s\mathrm{<mo>=</mo>1}}^k}T\mathrm{<mo\; stretchy="false">(</mo>}s\mathrm{<mo\; stretchy="false">)</mo>}h\mathrm{<mo\; stretchy="false">(</mo>}{\mu }_{j_l}\mathrm{<mo\; stretchy="false">(</mo>}{g}_s\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>.}$

More convenient below is to use the notations

$h_r\mathrm{<mo\; stretchy="false">(</mo>}{\mu }_i\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{\displaystyle \sum _{l\mathrm{<mo>=</mo>1}}^k}{\pi }_{l}h\mathrm{<mo\; stretchy="false">(</mo>}{\mu }_{j_r}\mathrm{<mo\; stretchy="false">(</mo>}{g}_l\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>,<mi>\; </mi><mi>\; </mi><mi>\; </mi><mi>\; </mi><mo\; stretchy="false">(</mo>}r\mathrm{<mo>=</mo>}f\mathrm{,}m\mathrm{<mo\; stretchy="false">)</mo>}$

and

$h_r\mathrm{<mo\; stretchy="false">(</mo>}{\mu }_i\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{\displaystyle \sum _{s\mathrm{<mo>=</mo>1}}^k}T\mathrm{<mo\; stretchy="false">(</mo>}s\mathrm{<mo\; stretchy="false">)</mo>}h\mathrm{<mo\; stretchy="false">(</mo>}{\mu }_{j_r}\mathrm{<mo\; stretchy="false">(</mo>}{g}_s\mathrm{<mo\; stretchy="false">)</mo><mo>|</mo>}{y}_p\mathrm{<mo\; stretchy="false">)</mo>,<mi>\; </mi><mi>\; </mi><mi>\; </mi><mi>\; </mi><mo\; stretchy="false">(</mo>}r\mathrm{<mo>=</mo>1,...,}{b}_i\mathrm{<mo\; stretchy="false">)</mo>.}$

Let y_j1<…<y_jkj be the k_j failures of type j(j=1,…,d), R(y_ji),be the risk set at y_ji, the partial likelihood is

$L\mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{<mo>|</mo>}\theta \mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{\displaystyle \prod _{j\mathrm{<mo>=</mo>1}}^d}{\displaystyle \prod _{i\mathrm{<mo>=</mo>1}}^{k_j}}\frac{h_{r_i}\mathrm{<mo\; stretchy="false">(</mo>}{\mu }_j\mathrm{<mo\; stretchy="false">)</mo>}}{{\displaystyle \sum _{l\in R\mathrm{<mo\; stretchy="false">(</mo>}{y}_{ji}\mathrm{<mo\; stretchy="false">)</mo>}}}{h}_{r_l}\mathrm{<mo\; stretchy="false">(</mo>}{\mu }_j\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{.<mo\; stretchy="false">(</mo>15<mo\; stretchy="false">)</mo>}$

Asymptotic heuristic

For IID data, various asymptotic results can be obtained. The results from the score function, the likelihood ratio statistic, and the MLE are equivalent. These results can be used to establish confidence intervals or hypothesis testing, etc. for θ. Here we are more interested in using the MLE. For general dependence model, usually the treatment is non-standard. But for our model, since the log-likelihood is in the form of several additive pieces, standard method can be used to derive the asymptotic distribution of the MLE. Let zi=(yi,di)(i=1,…,n). For the IID data, it is well known that under mild regularity conditions, the MLE ${\widehat{\theta }}_n$nis strongly consistent and asymptotically distributed normal with mean at the true parameter value θ0, and variance matrix given by the inverse of the Fisher information. Here the observations are unbalanced, the asymptotic variance is the Fisher information times a weight matrix. To derive it, we need some notations, and mainly concentrate on model (4).

Let Nip be the total number of parents with the i-th measurement non-missing (i=1,…,d), Niss be those for the siblings, N be the total number of individuals in the study, γN_ir=N_ir/N(r=p,s). Assume lim_{N→8γNir=γir>0} exists (r=p,s;i=1,…,d). Let Y_p and Y_j be general random vectors associated with a parent and sib respectively, and Δ_p and Δ_j be the corresponding random vectors.

For model (4), let

$H\mathrm{<mo\; stretchy="false">(</mo>}\theta \mathrm{<mo>|</mo>}{Y}_p\mathrm{,}{Y}_j\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{\Delta }_p\odot log\mathrm{<mo\; stretchy="false">(</mo>}{\displaystyle \sum _{l\mathrm{<mo>=</mo>1}}^k}{\pi }_{l}f\mathrm{<mo\; stretchy="false">(</mo>}{Y}_p\mathrm{<mo>|</mo>}\theta \mathrm{,}l\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}+\mathrm{<mo\; stretchy="false">(</mo>1}-{\Delta }_p\mathrm{<mo\; stretchy="false">)</mo>}\odot log\mathrm{<mo\; stretchy="false">(</mo>}{\displaystyle \sum _{l\mathrm{<mo>=</mo>1}}^k}{\pi }_{l}S{Y}_p\mathrm{<mo>|</mo>}\theta \mathrm{,}l\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}$

$+{\Delta }_j\odot log\mathrm{<mo\; stretchy="false">(</mo>}{\displaystyle \sum _{l\mathrm{<mo>=</mo>1}}^k}T\mathrm{<mo\; stretchy="false">(</mo>}l\mathrm{<mo\; stretchy="false">)</mo>}f\mathrm{<mo\; stretchy="false">(</mo>}{Y}_j\mathrm{<mo>|</mo>}G\mathrm{,}\theta \mathrm{,}l\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}+\mathrm{<mo\; stretchy="false">(</mo>1}-{\Delta }_j\mathrm{<mo\; stretchy="false">)</mo>}\odot log\mathrm{<mo\; stretchy="false">(</mo>}{\displaystyle \sum _{l\mathrm{<mo>=</mo>1}}^k}T\mathrm{<mo\; stretchy="false">(</mo>}l\mathrm{<mo\; stretchy="false">)</mo>}S\mathrm{<mo\; stretchy="false">(</mo>}{Y}_j\mathrm{<mo>|</mo>}G\mathrm{,}\theta \mathrm{,}l\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}\text{ (16)}$

here we use the notation Δ_p v(⋅) to represent the marginal version of v(⋅) corresponding to the non-zero components of Δ_p. The Fisher information matrix is

$I\mathrm{<mo\; stretchy="false">(</mo>}\theta \mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}-E\mathrm{<mo\; stretchy="false">(</mo>}\frac{{\partial }^{2}H\mathrm{<mo\; stretchy="false">(</mo>}\theta \mathrm{<mo>|</mo>}{Y}_p\mathrm{,}{Y}_j\mathrm{<mo\; stretchy="false">)</mo>}}{\partial {\theta }^{\prime }\partial \theta }\mathrm{<mo\; stretchy="false">)</mo>.<mo\; stretchy="false">(</mo>17<mo\; stretchy="false">)</mo>}$

The above expectation is more involved than it looks, since that involves summations of all possible combinations of non-zero elements of it with respect to Δ_r, and also the unknown distribution of it. Instead, an empirical version of it has a known form

$I_N\mathrm{<mo\; stretchy="false">(</mo>}\theta \mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}-\frac{1}{N}\frac{{\partial }^{2}logL\mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{<mo>|</mo>}\theta \mathrm{<mo\; stretchy="false">)</mo>}}{\partial {\theta }^{\prime }\partial \theta }{\mathrm{<mo>|</mo>}}_{\theta \mathrm{<mo>=</mo>}{\widehat{\theta }}_n}\mathrm{,<mo\; stretchy="false">(</mo>18<mo\; stretchy="false">)</mo>}$

where L(y\θ) is given by (4). At the true data generating parameter θ₀, I_N is strongly consistent for I(⋅). To obtain the weight matrix, we need to specify the parameter order in θ. We arrange the first k-1 entry to be π1,…,π_k-1, next we arrange all the regression parameters for the first response variable,..., all the regression parameters for the last response variable, then all the independent parameters in the variance matrix S and covariance matrix &Omega; in the similar order. It is clear that, in the weight matrix W, for (i,j) corresponding to the first k-1 components in θ, the weight should be γ_p=γ_1p+…+γ_dp; for (i,j) corresponding to the r-th and the l-th regression parameters, the weight is $\sqrt{(\Upsilon rp+\Upsilon rs)(\Upsilon lp+\Upsilon ls)}$; for (i,j) corresponding to the (a,b)-th and the (u,v)-th variance or covariance, the weight is [(γ_ap+γ_as) (γ_bp+γ_bs) (γ_up+γ_us) (γ_vp+γ_vs)]1/4. Let θ₀ be the true unknown data generating parameter, d→ stands for convergence in distribution. Then, we have

$\sqrt{N}\mathrm{<mo\; stretchy="false">(</mo>}{\widehat{\theta }}_N-{\theta }_0\mathrm{<mo\; stretchy="false">)</mo>}\stackrel{d}{\to }N\mathrm{<mo\; stretchy="false">(</mo>0,}{I}^{-1}\mathrm{<mo\; stretchy="false">(</mo>}{\theta }_0\mathrm{<mo\; stretchy="false">)</mo>}\otimes W\mathrm{<mo\; stretchy="false">)</mo>}\text{ (19)}$

where A⊗B stands for the Kronecker product of matrices A and B. Since I(⋅) involves unknown quantities, equivalently

$\sqrt{N}\mathrm{<mo\; stretchy="false">(</mo>}{\widehat{\theta }}_N-{\theta }_0\mathrm{<mo\; stretchy="false">)</mo>}:N\mathrm{<mo\; stretchy="false">(</mo>0,}{I}_N^{-1}\mathrm{<mo\; stretchy="false">(</mo>}{\widehat{\theta }}_n\mathrm{<mo\; stretchy="false">)</mo>}\otimes W_N\mathrm{<mo\; stretchy="false">)</mo>}\text{ (20)}$

where W_N is W with the γs replaced by the γNs.

The above ideal applies to the other models in this paper, but the results will be more involved, and we only discuss them briefly.

For the proportional hazards model, even for the IID data case, the conditional likelihood looks much different from the full likelihood. Interestingly, the MLE from this model (under the assumption of correct model specification and some regularity conditions) has the same asymptotic distribution as that from the full likelihood [21,52]. For the proportional hazards model, it is noted [4] that the survival function can be written as

S(y\x,β)=-S0(y)^h(β'x),

and f(y\x,β)=-dS(y\x,β)/dy. S_o to get the full log-likelihood (14), we need the estimate of S0(y) for the d-dimensional case [29,30]. Proposed the multi-dimensional generalization of the Kaplan-Meyer nonparametric estimator of S(y), similar technique can be used here for the construction of S₀(y). Then (18) continues to hold in this case. Due to technical involvement, we will not pursue the details here.

For the least squares estimator, since the weight involves the kernel and h_n, the treatment is different from those above, and in the case of full observation generally the asymptotic result is of the form

$\sqrt{n{h}_n^d}\mathrm{<mo\; stretchy="false">(</mo>}{\widehat{\theta }}_n-{\theta }_0-h_n^{2}C-o_p\mathrm{<mo\; stretchy="false">(</mo>}{h}_n^2\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}\stackrel{d}{\to }N\mathrm{<mo\; stretchy="false">(</mo>0,}\Omega \mathrm{<mo\; stretchy="false">)</mo>,}$

for some constant C and matrix Ω determined by the kernel and the true (unknown) data and censoring distributions [53,54]. In our case of partial observation, the above result holds with Ω replaced by Ω⊗W.

For the competing risks model, the structure is similar to that of the proportional hazards model. Here the response is one dimensional so that the survival function can be estimated by the Kaplan-Meyer estimator and the weight matrix W is the identity.

Discussion

We have considered several statistical methods, parametric, semi parametric and nonparametric models, for the genetic regression analysis of familial multiple endpoints data, with possible missing records. Here we only considered the case of nuclear families and the parameters are independent of time. The cases of arbitrary pedigrees and/or the time dependent parameter can be treated similarly. The variance components method can also be applied to the proportional hazards model and in the analysis of competing risks. There are some marginal models for the multiple endpoints data, which work well in practice. But we think the joint model is more appropriate when the within responses structure is important in the analysis. Another commonly used method to deal with the missing data is the EM algorithm [55], which can be implemented into the models considered here. But for the multiple endpoints data, the proportion of missing part is usually large; the EM algorithm may not be efficient. When the missing pattern is non-ignorable, more complicated approaches need to be considered to reduce potential biases. Hypothesis testing for parameters of interests can be conducted using the likelihood ratio statistics based on the parametric models. We only derived the basic forms of these models; more features can be implemented to them in particular applications.

References

1. Kaplan EL, Meier P. Nonparametric estimation from incomplete observations. Journal of the American Statistical Association. 1958; 53: 457-481.
2. Breslow N. Covariance analysis of censored survival data. Biometrics. 1974; 30: 89-99.
3. Fleming TR, Lin DY. Survival analysis in clinical trials: past developments and future directions. Biometrics. 2000; 56: 971-983.
4. Kalbfleisch JD, Prentice RL. The Statistical Analysis of Failure Time Data, Wiley. 1980.
5. Andersen PK, Gill RD. Cox’s regression model for counting processes: A large sample study. Annals of Statistics. 1982; 10: 1100-1120.
6. Arjas E. Survival models and martingale dynamics (with discussion). Scandinavian Journal of Statistics. 1989; 16: 177-225.
7. Flemming TR, Harrington DP. Counting Processes and Survival Analysis. Wiley, New York. 1991.
8. Andersen PK, Borgan O, Gill RD, Keiding N. Statistical Models Based on Counting Processes. Springer-Verlag. 1991
9. Printice R, Williams BJ, Peterson AV. On the regression analysis of multiple failure time data. Biometrika. 1981; 68: 373-379.
10. Wei LJ, Lin DY, Weissfeld L. Regression analysis of multivariate incomplete failure time data by modeling marginal distributions. Journal of the American Statistical Association. 1989; 84: 1065-1073.
11. Freedman LS, Graubard BI, Schatzkin A. Statistical validation of intermediate endpoints for chronic diseases. Stat Med. 1992; 11: 167-178.
12. Fleming TR, DeMets DL. Surrogate end points in clinical trials: are we being misled? Ann Intern Med. 1996; 125: 605-613.
13. Zhang J, Quan H, Ng J, Stepanavage ME. Some statistical methods for multiple endpoints in clinical trials. Controlled Clinical Trials. 1997; 18: 204-221.
14. Wei LJ, Glidden DV. An overview of statistical methods for multiple failure time data in clinical trials. Stat Med. 1997; 16: 833-839.
15. Remington DL. Effects of genetic and environmental factors on trait network predictions from quantitative trait locus data. Genetics. 2009; 181: 1087-1099.
16. Knüppel S, Rohde K, Meidtner K, Drogan D, Holzhütter HG, Boeing H, et al. Evaluation of 41 candidate gene variants for obesity in the EPIC-Potsdam cohort by multi-locus stepwise regression. PLoS One. 2013; 8: e68941.
17. Elston RC, Stewart J. A general model for the genetic analysis of pedigree data. Hum Hered. 1971; 21: 523-542.
18. Bonney GE. Compound regressive models for family data. Hum Hered. 1992; 42: 28-41.
19. Yuan A, Bonney GE. Two new recursive likelihood calculation methods for genetic analysis. Hum Hered. 2002; 54: 82-98.
20. Cox DR. Regression models and life-tables. Journal of the Royal Statistical Society. 1972; 34: 187-202.
21. Cox DR. Partial likelihood. Biometrika. 1975; 62: 269-276.
22. Clayton DG. A Monte Carlo method for Bayesian inference in frailty models. Biometrics. 1991; 47: 467-485.
23. Hougaard P. A class of multivariate failure time distributions. Biometrika. 1986; 73: 671-678.
24. Holt JD, Prentice RL. Survival analysis in twin studies and matched pairs experiments. Biometrika. 1974; 61: 17-30.
25. Oakes D. A model for association in bivariate survival data. Journal of the Royal Statistical Society. 1982; 44: 414-422.
26. Wild CJ. Failure time models with matched data. Biometrika. 1983; 70: 633-641.
27. Stute W. Consistent estimation under random censorship when covariable is present. Journal of Multivariate Analysis. 1993; 45: 89-103.
28. Stute W. Distributional convergence under random censorship when covariables are present.Scandinavian Journal of Statistics. 1996; 23: 462-471.
29. Dabrowska DM. Kaplan Meier estimate on the plane. Annals of Statistics. 1988; 16: 1475-1489.
30. Dabrowska DM. Kaplan Meier estimate on the plane: weak convergence, LIL and the bootstrap. Journal of Multivariate Analysis. 1988; 29: 308-325.
31. Gill R. Multivariate survival analysis. Theory of Probability and Its Applications. 1992a; 37: 18-31.
32. Gill R. Multivariate survival analysis: Part 2. Theory of Probability and Its Applications. 1992b; 37: 284-301.
33. Fan J. Locally linear regression smoothers and their minimax efficiencies. Annals of Statistics. 1993; 21: 196-216.
34. Ruppert D, Wand MP. Multivariate locally weighted least squares regression. Annals of Statistics. 1994; 22: 1346-1370.
35. Fan J, Gijbels I. Variable bandwidth and local linear regression smoothers. Annals of Statistics. 1992; 20: 2008-2036.
36. Fisher RA. The correlation between relatives on the supposition of Mendel inheritance. Trans Roy Soc. Edinb. 1918; 52: 399-433.
37. Harris DL. Genotypic Covariances Between Inbred Relatives. Genetics. 1964; 50: 1319-1348.
38. Lange K, Boehnke M. Extensions to pedigree analysis IV. Covariance components models for multivariate traits. Am J Med Genet. 1983; 14: 513-524.
39. Goldgar DE. Multipoint analysis of human quantitative genetic variation. Am J Hum Genet. 1990; 47: 957-967.
40. Amos CI. Robust variance-components approach for assessing genetic linkage in pedigrees. Am J Hum Genet. 1994; 54: 535-543.
41. Almasy L, Blangero J. Multipoint quantitative-trait linkage analysis in general pedigrees. Am J Hum Genet. 1998; 62: 1198-1211.
42. Fulker DW, Cherny SS, Sham PC, Hewitt JK. Combined linkage and association sib-pair analysis for quantitative traits. Am J Hum Genet. 1999; 64: 259-267.
43. Amos CI, Gu X, Chen J, Davis BR. Least squares estimation of variance components for linkage. Genet Epidemiol. 2000; 19: S1-7.
44. Blangero J, Williams JT, Almasy L. Robust LOD scores for variance component-based linkage analysis. Genet Epidemiol. 2000; 19: S8-14.
45. Sham PC, Purcell S. Equivalence between Haseman-Elston and variance-components linkage analyses for sib pairs. Am J Hum Genet. 2001; 68: 1527-1532.
46. de Andrade M, Guéguen R, Visvikis S, Sass C, Siest G, Amos CI. Extension of variance components approach to incorporate temporal trends and longitudinal pedigree data analysis. Genet Epidemiol. 2002; 22: 221-232.
47. Lange K. Statistics for Biology and Health. Springer. 1997.
48. Jacquard A. The genetic structure of populations. 1974.
49. Wright S. The genetical structure of populations. Ann Eugen. 1951; 15: 323-354.
50. Cockerham CC. Variance of gene frequencies. Evolution. 1969; 23: 72-84.
51. Yuan A, Chen G, Yang Q, Rotimi C, Bonney G. Variance components model with disequilibria. Eur J Hum Genet. 2006; 14: 941-952.
52. Tsiatis A. A large sample study of the estimate for the integrated hazard function in Cox’s regression model for survival data. Technical Report No. 562. 1978; 172-180.
53. Fan J, Hu I, Truong Y. Robust nonparametric function estimation. Scandinavian Journal of Statistics. 1994; 21: 433-446.
54. Cai Z. Weighted local linear approach to censored nonparametric regression. Akritas MG, Politis DM. In: Recent Advances and Trends in Nonparametric Statistics. 2003; 217-231.
55. Dempster AP, Laird NM, Rubin DB. Maximum likelihood from incomplete data via the EM algorithm. Journal of the Royal Statistical Society. 1977; 39: 1-38.