solvewtsc: SOLVING THE WEIGHTED SCORES EQUATIONS WITH INPUTS OF THE WEIGHT MATRICES AND THE DATA

Description

Solving the weighted scores equations with inputs of the weight matrices and the data.

Usage

solvewtsc(start,WtScMat,xdat,ydat,id,tvec,margmodel,link)
solvewtsc.ord(start,WtScMat,xdat,ydat,id,tvec,link)

Arguments

start

A starting value of the vector of regression and not regression parameters. The CL1 estimates of regression and not regression parameters is a good starting value.

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

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.

Value

root: The weighted scores estimates.
f.root: The value of the wtsc function evaluated at the root.
iter: The number of iterations used.
estim.precis: The estimated precision for root.

Details

Obtain robust estimates ${\hat a}$ of the univariate parameters solving the weighted scores equation, $ g_1= g_1(a)=\sum_{i=1}^n X_i^T\, W_{i,\rm working}^{-1}\, s_i( a)=0, $ where $ W_{i,\rm working}^{-1}=\Delta_i\Omega_{i,\rm working}^{-1}= \Delta_i({\tilde a})\Omega_i({\tilde a},{\tilde R})^{-1}$ is based on the covariance matrix of $s_i(a)$ computed from the fitted discretized MVN model with estimated parameters ${\tilde a}, {\tilde R}$. A reliable non-linear system solver is used.

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

References

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

################################################################################
#                           NB2 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="nb2"
################################################################################
#                      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,link)
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)
################################################################################
#                         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="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)
# ## End(Not run)

Run the code above in your browser using DataLab