weightedScores (version 0.9.5.3)

cl1: OPTIMIZATION ROUTINE FOR BIVARIATE COMPOSITE LIKELIHOOD FOR MVN COPULA

Description

Optimization routine for bivariate composite likelihood for MVN copula.

Usage

cl1(b,gam,xdat,ydat,id,tvec,margmodel,corstr,link)
cl1.ord(b,gam,xdat,ydat,id,tvec,corstr,link)

Arguments

b

The regression coefficients.

gam

The uinivariate parameters that are not regression coefficients. That is the parameter \(\gamma\) of negative binomial distribution or the \(q\)-dimensional vector of the univariate cutpoints of ordinal model. \(\gamma\) is NULL for Poisson and binary regression.

xdat

\((\mathbf{x}_1 , \mathbf{x}_2 , \ldots , \mathbf{x}_n )^\top\), where the matrix \(\mathbf{x}_i,\,i=1,\ldots,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 including the unit first column to account for the intercept (except for ordinal regression where there is no intercept). This xdat matrix is of dimension \((N\times p),\) where \(N =\sum_{i=1}^n j_i\) is the total number of observations from all units.

ydat

\((y_1 , y_2 , \ldots , y_n )^\top\), where the response data vectors \(y_i,\,i=1,\ldots,n\) are of possibly different lengths for different units. In particular, we now have that \(y_i\) is (\(j_i \times 1\)), where \(j_i\) is the number of observations on unit \(i\). The total number of observations from all units is \(N =\sum_{i=1}^n j_i\). The ydat are the collection of data vectors \(y_i, i = 1,\ldots,n\) one from each unit which summarize all the data together in a single, long vector of length \(N\).

id

An index for individuals or clusters.

tvec

A vector with the time indicator of individuals or clusters.

margmodel

Indicates the marginal model. Choices are “poisson” for Poisson, “bernoulli” for Bernoulli, and “nb1” , “nb2” for the NB1 and NB2 parametrization of negative binomial in Cameron and Trivedi (1998).

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 “log” for the log link function, “logit” for the logit link function, and “probit” for the probit link function.

Value

A list containing the following components:

minimum

The negative value of the sum of bivariate marginal log-likelihoods at CL1 estimates.

estimate

The CL1 estimates.

gradient

The gradient at the estimated minimum of CL1.

code

An integer indicating why the optimization process terminated, same as in nlm.

Details

The CL1 composite likelihood method in Zhao and Joe (2005). The univariate parameters are estimated from the sum of univariate marginal log-likelihoods and then the dependence parameters are estimated from the sum of bivariate marginal log-likelihoods with the univariate parameters fixed from the first step.

Note that bcl.ord is a variant of the code for ordinal (probit and logistic) regression.

References

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

See Also

bcl iee

Examples

# NOT RUN {
################################################################################
#                      NB1 regression for count data
################################################################################
################################################################################
#                      read and set up data set
################################################################################
data(childvisit)
# covariates
season1<-childvisit$q
season1[season1>1]<-0
xdat<-cbind(1,childvisit$sex,childvisit$age,childvisit$m,season1)
# response
ydat<-childvisit$hosp
#id
id<-childvisit$id
#time
tvec<-childvisit$q
################################################################################
#                      select the marginal model
################################################################################
margmodel="nb1"
################################################################################
#                      select the  correlation structure
################################################################################
corstr="exch"
################################################################################
#                      perform CL1 estimation
################################################################################
i.est<-iee(xdat,ydat,margmodel)
cat("\niest: IEE estimates\n")
print(c(i.est$reg,i.est$gam))
est.rho<-cl1(b=i.est$reg,gam=i.est$gam,xdat,ydat,id,tvec,margmodel,corstr)
cat("\nest.rho: CL1 estimates\n")
print(est.rho$e)
################################################################################
#                         Ordinal regression 
################################################################################
################################################################################
#                      read and set up data set
################################################################################
# }
# NOT RUN {
data(arthritis)
nn=nrow(arthritis)
bas2<-bas3<-bas4<-bas5<-rep(0,nn)
bas2[arthritis$b==2]<-1
bas3[arthritis$b==3]<-1
bas4[arthritis$b==4]<-1
bas5[arthritis$b==5]<-1
t2<-t3<-rep(0,nn)
t2[arthritis$ti==3]<-1
t3[arthritis$ti==5]<-1
xdat=cbind(t2,t3,arthritis$trt,bas2,bas3,bas4,bas5,arthritis$age) 
ydat=arthritis$y
id<-arthritis$id
#time
tvec<-arthritis$time
################################################################################
#                      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)
# }