washb_lincom

Usage

washb_lincom(lc = NULL, fit, vcv, measure = "RR", flag = NULL)

Arguments

lc

: Index vector of coefficients from glm modellinear combination of coefficients

fit

: Model object returned from coeftest (class=coeftest) command within the washb_glm function. Accessed with $fit from washb_glm output.

vcv

: variance-covariance matrix of coefficients. Accessed with $vcv from washb_glm output.

measure

measure of effect. RR = risk ratio, RD = risk difference

flag

Internal argument used to flag and suppress printing if the washb_lincom function is called within another function.

Examples

library(washb)
data(washb_anthro)
Warning message:
data set ‘washb_anthro’ not found
data(washb_bd_enrol)

 # drop svydate and month because they are superceded in the child level diarrhea data
washb_bd_enrol$svydate <- NULL
washb_bd_enrol$month <- NULL

# merge the baseline dataset to the follow-up dataset
ad <- merge(washb_bd_enrol,washb_bd_diar,by=c("dataid","clusterid","block","tr"),all.x=F,all.y=T)

# subset to the relevant measurement
# Year 1 or Year 2
ad <- subset(ad,svy==1|svy==2)

#subset the diarrhea to children <36 mos at enrollment
### (exlude new births that are not target children)
ad <- subset(ad,sibnewbirth==0)
ad <- subset(ad,gt36mos==0)

# Exclude children with missing data
ad <- subset(ad,!is.na(ad$diar7d))

#Re-order the tr factor for convenience
ad$tr <- factor(ad$tr,levels=c("Control","Water","Sanitation","Handwashing","WSH","Nutrition","Nutrition + WSH"))

#Ensure that month is coded as a factor
ad$month <- factor(ad$month)

#Sort the data for perfect replication when using V-fold cross-validation
ad <- ad[order(ad$block,ad$clusterid,ad$dataid,ad$childid),]

#Estimate subgroup analysis glm with washb_glm
glm.C.N.byChildType <- washb_glm(Y=ad$diar7d,tr=ad$tr,pair=ad$block, W=W_tchild, V="tchild", id=ad$clusterid, contrast=c("Control","Nutrition"), family=binomial(link='log'), print=FALSE)


#The lincom function is called within washb_glm, so for standard subgroup analysis, there is no need to use the lincom()
#function separately. But it can be used on its own to create a custom combination of regression coefficients


#Examine the treatment effect across subgroups with `objectname'$lincom
glm.C.N.byChildType$lincom
  Tr vs. C by Subgroup    est se.est est.lb est.ub       Z      P
1              Sibling 0.6422 0.2986 0.3577 1.1531 -1.4830 0.1381
2         Target child 1.2734 0.1475 0.9537 1.7003  1.6387 0.1013


#Below is an example of using the `lincom()` function to externally replicate the target child/sibling subgroup analysis of the diarrheal disease outcome. First, recall the glm model contrasting diarrheal disease prevalence ratio between target child and sibling fitted to diarrheal disease:
#Use the lincom function to verify the subgroup analysis calculated within the `washb_glm()` function. `washb_lincom()` takes as argements the fit and variance-covariance matrix returned by the glm function, and an index vector `lc`.
#`lc`= vector indexing coefficients to be linearly combined. If a single coefficient is chosen, it will return the same coefficient and confidence interval as the glm coefficient. Example:
#Create lc vector of 0's equal in length to the number of coefficients from the glm model.
lc=rep(0,nrow(glm.C.N.byChildType$fit))
#Examine model coefficients (minus the pair-matched block estimates) to determine the position of coefficients to combine.
glm.C.N.byChildType$fit[1:6,]
                      PR       2.5%       97.5%    Estimate Std. Error     z value       Pr(>|z|)
(Intercept)   0.01335657 0.00229925  0.07758973 -4.31574660  0.8976665 -4.80773914 0.000001526468
trNutrition   0.64223112 0.35771147  1.15305449 -0.44280703  0.2985824 -1.48303118 0.138066125872
VTarget child 1.39433404 0.96442019  2.01589249  0.33241691  0.1880842  1.76738315 0.077164083197
pair2         1.07697676 0.11608458  9.99167117  0.07415782  1.1365276  0.06524946 0.947975377785
pair3         1.15694133 0.11228806 11.92035232  0.14577973  1.1900344  0.12250043 0.902502702220
pair4         4.55031058 0.68843416 30.07597167  1.51519549  0.9635363  1.57253603 0.115826290776

#Replace the second and fourth position in the vector with 1 (the positions of the treatment coefficient and interaction term in the model)
lc[c(2,4)]<-1

#Run the lincom function and compare output to the $lincom output from the GLM model.
washb_lincom(lc=lc,fit=glm.C.N.byChildType$fit,vcv=glm.C.N.byChildType$vcv, measure="RR")

Linear combination of coefficients:
[1] "trNutrition + pair2"

         RR  se.RR  RR.lb RR.ub       Z      P
[1,] 0.6917 1.1671 0.0702 6.813 -0.3159 0.7521


# The `lincom()` function can also be used to calculate risk difference
#Calculate risk difference using the risk difference output from the washb_glm function:
washb_lincom(lc=lc,fit=glm.C.N.byChildType$RDfit,vcv=glm.C.N.byChildType$vcvRD, measure="RD")
Error in t(lc) %*% x[, 1]: requires numeric/complex matrix/vector arguments
#Compare to RD by subgroup from the washb_glm function:
glm.C.N.byChildType$lincomRD
NULL