godambe: INVERSE GODAMBE MATRIX

Description

Inverse Godambe matrix.

Usage

godambe(param,WtScMat,xdat,ydat,id,tvec,margmodel,link)
godambe.ord(param,WtScMat,xdat,ydat,id,tvec,link)

Arguments

param

The weighted scores estimates of regression and not regression parameters.

WtScMat

A list containing the following components. omega: The array with the $\Omega_i,\,i=1,\ldots,n$ matrices; delta: The array with the $\Delta_i,\,i=1,\ldots,n$ matrices; X: The array with the $X_i,\,i=1,\ldots,n$ matrices.

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

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 observation

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 Triv

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

The inverse Godambe matrix.

Details

If the $W_{i,\rm working}$ are assumed fixed for the second stage of solving the weighted scores equations $$g_1= g_1(a)=\sum_{i=1}^n X_i^T\, W_{i,\rm working}^{-1}\, s_i( a)=0,$$ the asymptotic covariance matrix of the solution $\widehat a_1$ is $$V_1=(-D_{g_1})^{-1}M_{g_1}(-D^T_{g_1})^{-1}$$ with $$-D_{g_1} =\sum_{i=1}^n X_i^T W_{i,\rm working}^{-1}\Delta_i X_i,$$ $$M_{ g_1} = \sum_{i=1}^n X_i^T W_{i,\rm working}^{-1} \Omega_{i,\rm true}( W_{i,\rm working}^{-1})^T X_i,$$ where $\Omega_{i,\rm true}$ is the true covariance matrix of $s_i(a)$. The inverse of $V_1$ is known as Godambe information matrix (Godambe, 1991).

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

References

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.

Nikoloulopoulos, A.K. (2015a) Correlation structure and variable selection in generalized estimating equations via composite likelihood information criteria. Arxiv e-prints.

Nikoloulopoulos, A.K. (2015b) Weighted scores estimating equations for longitudinal ordinal data. Arxiv e-prints.

Examples

Run this code

################################################################################
#                           Poisson regression 
################################################################################
################################################################################
#                      read and set up the 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="poisson"
################################################################################
#                      select the  correlation structure
################################################################################
corstr="exch"
################################################################################
#                      perform CL1 estimation
################################################################################
i.est<-iee(xdat,ydat,margmodel)
cat("iest: IEE estimates
")
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("est.rho: CL1 estimates
")
print(est.rho$e)
################################################################################
#                      obtain the fixed weight matrices
################################################################################
WtScMat<-weightMat(b=i.est$reg,gam=i.est$gam,rh=est.rho$e,
xdat,ydat,id,tvec,margmodel,corstr)
################################################################################
#                      obtain the weighted scores estimates
################################################################################
# solve the nonlinear system of equations
ws<-solvewtsc(start=c(i.est$reg,i.est$gam),WtScMat,xdat,ydat,id,
tvec,margmodel,link)
cat("ws=parameter estimates
")
print(ws$r)
################################################################################
#                      obtain the inverse Godambe matrix
################################################################################
acov<-godambe(ws$r,WtScMat,xdat,ydat,id,tvec,margmodel)
cat("acov: inverse Godambe matrix with W based on first-stage wt
matrices
")
print(acov)
################################################################################
#                         Ordinal regression 
################################################################################
################################################################################
#                      read and set up data set
################################################################################
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="probit"
################################################################################
#                      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)
################################################################################
#                      obtain the fixed weight matrices
################################################################################
WtScMat<-weightMat.ord(b=i.est$reg,gam=i.est$gam,rh=est.rho$e,xdat,ydat,id,
tvec,corstr,link)
################################################################################
#                      obtain the weighted scores estimates
################################################################################
# solve the nonlinear system of equations
ws<-solvewtsc.ord(start=c(i.est$reg,i.est$gam),WtScMat,xdat,ydat,id,
tvec,link)
cat("ws=parameter estimates\n")
print(ws$r)
################################################################################
#                      obtain the inverse Godambe matrix
################################################################################
acov<-godambe.ord(ws$r,WtScMat,xdat,ydat,id,tvec,link)
cat("\nacov: inverse Godambe matrix with W based on first-stage wt
matrices\n")
print(acov)

Run the code above in your browser using DataLab