godambe {weightedScores}R Documentation

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 Ωi,i=1,,n\Omega_i,\,i=1,\ldots,n matrices; delta: The array with the Δi,i=1,,n\Delta_i,\,i=1,\ldots,n matrices; X: The array with the Xi,i=1,,nX_i,\,i=1,\ldots,n matrices.

xdat

(x1,x2,,xn)(\mathbf{x}_1 , \mathbf{x}_2 , \ldots , \mathbf{x}_n )^\top, where the matrix xi,i=1,,n\mathbf{x}_i,\,i=1,\ldots,n for a given unit will depend on the times of observation for that unit (jij_i) and will have number of rows jij_i, each row corresponding to one of the jij_i elements of yiy_i and pp columns where pp 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×p),(N\times p), where N=i=1njiN =\sum_{i=1}^n j_i is the total number of observations from all units.

ydat

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

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).

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.

Details

If the Wi,workingW_{i,\rm working} are assumed fixed for the second stage of solving the weighted scores equations

g1=g1(a)=i=1nXiTWi,working1si(a)=0, 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 a^1\widehat a_1 is

V1=(Dg1)1Mg1(Dg1T)1 V_1=(-D_{g_1})^{-1}M_{g_1}(-D^T_{g_1})^{-1}

with

Dg1=i=1nXiTWi,working1ΔiXi, -D_{g_1} =\sum_{i=1}^n X_i^T W_{i,\rm working}^{-1}\Delta_i X_i,

Mg1=i=1nXiTWi,working1Ωi,true(Wi,working1)TXi, 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 Ωi,true\Omega_{i,\rm true} is the true covariance matrix of si(a)s_i(a). The inverse of V1V_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.

Value

The inverse Godambe matrix.

Author(s)

Aristidis K. Nikoloulopoulos A.Nikoloulopoulos@uea.ac.uk
Harry Joe harry.joe@ubc.ca

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. 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.

See Also

wtsc, solvewtsc, weightMat, wtsc.wrapper

Examples

################################################################################
#                           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("\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)
################################################################################
#                      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\n")
print(ws$r)
################################################################################
#                      obtain the inverse Godambe matrix
################################################################################
acov<-godambe(ws$r,WtScMat,xdat,ydat,id,tvec,margmodel)
cat("\nacov: inverse Godambe matrix with W based on first-stage wt
matrices\n")
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)


[Package weightedScores version 0.9.5.3 Index]