clogit_qgcomp.R

# using qgcomp.cox.noboot in a nested case-control analysis as a way to fit a conditional logit qgcomp model
library(qgcomp)
library(survival)

# read in data from a nested case control study
data(infert, package = "datasets")
head(infert)


### 1) check that it works, part A
# fit a conditional logistic model using the standard clogit model and checking against the qgcomp.cox.noboot function

infert$faketime = 1 # fake failure time to trick coxph

# checking that clogit gives the same answer
# this is just a trick to get qgcomp to run a standard Cox model (q=NULL, single exposure)
uvl.check = qgcomp.cox.noboot(Surv(faketime, case) ~ spontaneous + induced + strata(stratum), 
                              expnms = "spontaneous", data = infert, q=NULL)


uvl = clogit(case ~ spontaneous + induced + strata(stratum), data=infert)

# these should be equal
uvl$coefficients[1]
coef(uvl.check)

### 1) check that it works, part B

# now simulate some fake exposures just for this example 
# (the original data don't have continuous exposures)
N = nrow(infert)
set.seed(12312)
X <- data.frame(qgcomp:::.rmvnorm(N, c(0,0,0), diag(c(1,1,1))))
names(X) <- c("x.1", "x.2", "x.3")
infertX = cbind(infert,X)

# now running mixtures analysis: first a multivariable model
mcl = clogit(case ~ x.1 + x.2 + x.3 + spontaneous + induced + strata(stratum), data=infertX)

# now checking that qgcomp will fit the correct underlying model with non-quantized variables (q=NULL)
# note this estimates effect of 1 unit increase of all exposures (not 1 quantile)
mcl.check = qgcomp.cox.noboot(Surv(faketime, case) ~ x.1 + x.2 + x.3 + 
                                spontaneous + induced + strata(stratum), 
                              expnms = c("x.1", "x.2", "x.3"), data = infertX, q=NULL)

# these should be equal
sum(mcl$coefficients[1:3])
coef(mcl.check)

### 1) How to use qgcomp in practice to fit a conditional logistic model to matched case-control studies

# now running actual qgcomp 
# note this estimates effect of 1 quantile increase of all exposures 
qgc = qgcomp.cox.noboot(Surv(faketime, case) ~ x.1 + x.2 + x.3 + 
                          spontaneous + induced + strata(stratum), 
                        expnms = c("x.1", "x.2", "x.3"), data = infertX, 
                        q=4)
print(qgc)