CLIC: CL1 INFORMATION CRITERIA
In weightedScores: Weighted Scores Method for Regression Models with Dependent Data

Description Usage Arguments Details Value Author(s) References See Also Examples

Composite likelihood (CL1) information criteria.

1 2	clic1dePar(nbcl,r,b,gam,xdat,id,tvec,corstr,WtScMat,link,mvncmp) clic(nbcl,r,b,gam,xdat,id,tvec,corstr,WtScMat,link,mvncmp)

`nbcl`	The negative value of the sum of bivariate marginal log-likelihoods at CL1 estimates.
`r`	The CL1 estimates of the latent correlations.
`b`	The CL1 estimates of the regression coefficients.
`gam`	The CL1 estimates of the cutpoints.
`xdat`	(\mathbf{x}_1 , \mathbf{x}_2 , … , \mathbf{x}_n )^\top, where the matrix \mathbf{x}_i,\,i=1,…,n for a given unit will depend on the times of observation for that unit (j_i) and will have number of rows j_i, each row corresponding to one of the j_i elements of y_i and p columns where p is the number of covariates. This xdat matrix is of dimension (N\times p), where N =∑_{i=1}^n j_i is the total number of observations from all units.
`id`	An index for individuals or clusters.
`tvec`	A vector with the time indicator of individuals or clusters.
`WtScMat`	A list containing the following components. omega: The array with the Ω_i,\,i=1,…,n matrices; delta: The array with the Δ_i,\,i=1,…,n matrices; X: The array with the X_i,\,i=1,…,n matrices.
`corstr`	Indicates the latent correlation structure of normal copula. Choices are “exch”, “ar”, and “unstr” for exchangeable, ar(1) and unstructured correlation structure, respectively.
`link`	The link function. Choices are “logit” for the logit link function, and “probit” for the probit link function.
`mvncmp`	The method of computation of the MVN rectangle probabilities. Choices are 1 for `mvnapp` (faster), and 2 for `pmvnorm` (more accurate).

First, consider the sum of univariate log-likelihoods

L_1= ∑_{i=1}^n∑_{j=1}^d\, \log f_1(y_{ij};ν_{ij},\boldsymbol{γ})=∑_{i=1}^n∑_{j=1}^d\, \ell_1(ν_{ij},\boldsymbol{γ},\, y_{ij}),

and then the sum of bivariate log-likelihoods

L_2=∑_{i=1}^{n}∑_{j<k} \log{f_2(y_{ij},y_{ik};ν_{ij},ν_{ik},\boldsymbol{γ},ρ_{jk})}=∑_{i=1}^{n}∑_{j<k}\ell_2(ν_{ij},ν_{ik},\boldsymbol{γ},ρ_{jk};y_{ij},y_{ik}),

where \ell_2(\cdot)=\log f_2(\cdot) and

f_2(y_{ij},y_{ik};ν_{ij},ν_{ik},\boldsymbol{γ},ρ_{jk})=\int_{Φ^{-1}[F_1(y_{ij}-1;ν_{ij},\boldsymbol{γ})]}^{Φ^{-1}[F_1(y_{ij};ν_{ij},\boldsymbol{γ})]} \int_{Φ^{-1}[F_1(y_{ik}-1;ν_{ik}),\boldsymbol{γ}]}^{Φ^{-1}[F_1(y_{ik};ν_{ik},\boldsymbol{γ})]} φ_2(z_j,z_d;ρ_{jk}) dz_j dz_k;

φ_2(\cdot;ρ) denotes the standard bivariate normal density with correlation ρ.

Let \mathbf{a}^\top=(\boldsymbol{β}^\top,\boldsymbol{γ}^\top) be the column vector of all r=p+q univariate parameters. Differentiating L_1 with respect to \mathbf{a} leads to the independent estimating equations or univariate composite score functions:

\mathbf{g}_1=\mathbf{g}_1(\mathbf{a})= \frac{\partial L_1}{\partial \mathbf{a}}= ∑_{i=1}^n\mathbf{X}_i^\top\mathbf{s}_i^{(1)}(\mathbf{a})=\bf 0,

Differentiating L_2 with respect to \mathbf{R}=\bigl(ρ_{jk},1≤q j<k≤q d\bigr) leads to the bivariate composite score functions (Zhao and Joe, 2005):

\mathbf{g}_2=\frac{\partial L_2}{\partial \mathbf{R}}= ∑_{i=1}^{n}\mathbf{s}_i^{(2)}(\mathbf{a},\mathbf{R})= ∑_{i=1}^n\Bigl(\mathbf{s}_{i,jk}^{(2)}(\mathbf{a},ρ_{jk}),1≤q j<k≤q d\Bigr)=\mathbf{0},

where \mathbf{s}_i^{(2)}(\mathbf{a},\mathbf{R})=\frac{\partial∑_{j<k}\ell_2(ν_{ij},ν_{ik},\boldsymbol{γ},ρ_{jk};y_{ij},y_{ik})}{\partial \mathbf{R}} and \mathbf{s}_{i,jk}^{(2)}(\boldsymbol{γ},ρ_{jk})=\frac{\partial \ell_2(ν_{ij},ν_{ik},\boldsymbol{γ},ρ_{jk};y_{ij},y_{ik})}{\partial ρ_{jk}}. The CL1 estimates \widetilde{\mathbf{a}} and \widetilde{\mathbf{R}} of the discretized MVN model are obtained by solving the above CL1 estimating functions.

The asymptotic covariance matrix for the estimator that solves them, also known as the inverse Godambe (Godambe, 1991) information matrix, is

\mathbf{V}=(-\mathbf{H}_\mathbf{g})^{-1}\mathbf{J}_\mathbf{g}(-\mathbf{H}_\mathbf{g}^\top)^{-1},

where \mathbf{g}=(\mathbf{g}_1,\mathbf{g}_2)^\top. First set \boldsymbol{θ}=(\mathbf{a},\mathbf{R})^\top, then

-\mathbf{H}_\mathbf{g}=E\Bigl(\frac{\partial \mathbf{g}}{\partial\boldsymbol{θ}}\Bigr)= ≤ft[\begin{array}{cc} E\Bigl(\frac{\partial \mathbf{g}_1}{\partial\mathbf{a}}\Bigr) & E\Bigl(\frac{\partial \mathbf{g}_1}{\partial\mathbf{R}}\Bigr)\\ E\Bigl(\frac{\partial \mathbf{g}_2}{\partial\mathbf{a}}\Bigr) & E\Bigl(\frac{\partial \mathbf{g}_2}{\partial\mathbf{R}}\Bigr) \end{array}\right]= ≤ft[\begin{array}{cc} -\mathbf{H}_{\mathbf{g}_1}&\mathbf{0}\\ -\mathbf{H}_{\mathbf{g}_{2,1}}&-\mathbf{H}_{\mathbf{g}_2} \end{array}\right],

where -\mathbf{H}_{\mathbf{g}_1}=∑_i^n\mathbf{X}_i^\top\boldsymbol{Δ}_i^{(1)} \mathbf{X}_i, -\mathbf{H}_{\mathbf{g}_{2,1}}=∑_i^n\boldsymbol{Δ}_i^{(2,1)}\mathbf{X}_i, and -\mathbf{H}_{\mathbf{g}_2}=∑_i^n\boldsymbol{Δ}_i^{(2,2)}.

The covariance matrix \mathbf{J}_\mathbf{g} of the composite score functions \mathbf{g} is given as below

\mathbf{J}_\mathbf{g}=\mbox{Cov}(\mathbf{g})= ≤ft[\begin{array}{cc} \mbox{Cov}(\mathbf{g}_1) & \mbox{Cov}(\mathbf{g}_1,\mathbf{g}_2)\\ \mbox{Cov}(\mathbf{g}_2,\mathbf{g}_1) & \mbox{Cov}(\mathbf{g}_2) \end{array}\right]= ≤ft[\begin{array}{cc}\mathbf{J}_\mathbf{g}^{(1)}& \mathbf{J}_\mathbf{g}^{(1,2)}\\ \mathbf{J}_\mathbf{g}^{(2,1)}& \mathbf{J}_\mathbf{g}^{(2)}\end{array}\right]= ∑_i≤ft[\begin{array}{cc} \mathbf{X}_i^\top\boldsymbol{Ω}_i^{(1)} \mathbf{X}_i & \mathbf{X}_i^\top\boldsymbol{Ω}_i^{(1,2)}\\ \boldsymbol{Ω}_i^{(2,1)}\mathbf{X}_i& \boldsymbol{Ω}_i^{(2)}\end{array}\right],

where

≤ft[\begin{array}{cc} \boldsymbol{Ω}_i^{(1)}& \boldsymbol{Ω}_i^{(1,2)}\\ \boldsymbol{Ω}_i^{(2,1)}& \boldsymbol{Ω}_i^{(2)} \end{array}\right]= ≤ft[\begin{array}{cc} \mbox{Cov}\Bigl(\mathbf{s}_i^{(1)}(\mathbf{a})\Bigr) & \mbox{Cov}\Bigl(\mathbf{s}_i^{(1)}(\mathbf{a}), \mathbf{s}_i^{(2)}(\mathbf{a},\mathbf{R})\Bigr)\\ \mbox{Cov}\Bigl(\mathbf{s}_i^{(2)}(\mathbf{a},\mathbf{R}), \mathbf{s}_i^{(1)}(\mathbf{a})\Bigr) & \mbox{Cov}\Bigl(\mathbf{s}_i^{(2)}(\mathbf{a},\mathbf{R})\Bigr) \end{array}\right].

To this end, the composite AIC (Varin and Vidoni, 2005) and BIC (Gao and Song, 2011) criteria have the forms:

\mbox{CL1AIC} = -2L_2 + 2\mbox{tr}\Bigl(\mathbf{J}_\mathbf{g}\mathbf{H}_\mathbf{g}^{-1}\Bigr),

\mbox{CL1BIC} =-2L_2+\log(n)\mbox{tr}\Bigl(\mathbf{J}_\mathbf{g}\mathbf{H}_\mathbf{g}^{-1}\Bigr).

A list containing the following components:

`AIC`	The CL1AIC.
`BIC`	The CL1BIC.

Aristidis K. Nikoloulopoulos A.Nikoloulopoulos@uea.ac.uk

Gao, X. and Song, P.X.K. (2011). Composite likelihood EM algorithm with applications to multivariate hidden Markov model. Statistica Sinica 21, 165–185.

Godambe, V. P. (1991) Estimating Functions. Oxford: Oxford University Press

Nikoloulopoulos, A.K., Joe, H. and Chaganty, N.R. (2011) Weighted scores method for regression models with dependent data. Biostatistics, 12, 653–665. doi: 10.1093/biostatistics/kxr005.

Nikoloulopoulos, A.K. (2016) Correlation structure and variable selection in generalized estimating equations via composite likelihood information criteria. Statistics in Medicine, 35, 2377–2390. doi: 10.1002/sim.6871.

Nikoloulopoulos, A.K. (2017) Weighted scores method for longitudinal ordinal data. Arxiv e-prints, <arXiv:1510.07376>. https://arxiv.org/abs/1510.07376.

Varin, C. and Vidoni, P. (2005). A note on composite likelihood inference and model selection. Biometrika 92, 519–528.

Zhao, Y. and Joe, H. (2005) Composite likelihood estimation in multivariate data analysis. The Canadian Journal of Statistics, 33, 335–356.

solvewtsc, wtsc.wrapper

  ################################################################################
  #                           Binary regression 
  ################################################################################
  ################################################################################
  #                      read and set up the data set
  ################################################################################
  data(childvisit)
  # covariates
  season1<-childvisit$q
  season1[season1>1]<-0
  xdat<-cbind(childvisit$sex,childvisit$age,childvisit$m,season1)
  # response
  ydat<-childvisit$hosp
  ydat[ydat>0]=1
  ydat=2-ydat
  #id
  id<-childvisit$id
  #time
  tvec<-childvisit$q
  ################################################################################
  #                      select the link
  ################################################################################
  link="logit"
  ################################################################################
  #                      select the  correlation structure
  ################################################################################
  corstr="exch"
  ################################################################################
  #                      perform CL1 estimation
  ################################################################################
  i.est<-iee.ord(xdat,ydat,link)
  cat("\niest: IEE estimates\n")
  print(c(i.est$reg,i.est$gam))
  est.rho<-cl1.ord(b=i.est$reg,gam=i.est$gam,xdat,ydat,id,tvec,corstr,link)
  cat("\nest.rho: CL1 estimates\n")
  print(est.rho$e)
  cat("\nest.rho: negative CL1 log-likelhood\n")
  print(est.rho$m)
  ################################################################################
  #                      obtain the fixed weight matrices
  ################################################################################
  WtScMat<-weightMat.ord(b=i.est$reg,gam=i.est$gam,rh=0.1961,xdat,ydat,id,
                         tvec,corstr,link)
  ################################################################################
  #                      obtain the CL1 information criteria
  ################################################################################
  out<-clic1dePar(nbcl=est.rho$m,r=est.rho$e,i.est$r,i.est$g,xdat,id,tvec,corstr,WtScMat,link,1)