######################################################
## Examples for binary responses
###########################################
## Wheeze among Steubenville (see [3]):
## Latent Beta-distributed propensity
data(wheeze)
fit1 <- drm(wheeze~I(age>9)+smoking+cluster(id),data=wheeze,dep="B", print=0)
## Obesity among Muscatine children (see [2]):
## Analysis for completers: M2 for girls
data(obese)
fit2 <- drm(obese~age+cluster(id)+Time(year), subset=sex=="female",
dep="M2",data=obese)
## Not run:
# ## Muscatine children continued (see [3]):
# ## LM for boys and girls separately
# fit3 <- drm(obese~age+cluster(id)+Time(age), subset=sex=="male",
# dep="LM",data=obese)
#
# fit4 <- drm(obese~age+cluster(id)+Time(age), subset=sex=="female",
# dep="LM",data=obese)
# ## End(Not run)
############################################
## Examples for ordinal responses
############################################
## Movie critic example (see [6]):
## Latent Dirichlet propensities with baseline category link.
data(movie)
options(contrasts=c("contr.treatment","contr.treatment"))
fit5 <- drm(y~critic+cluster(movie), data=movie, dep="D", link="bcl")
## Longitudinal dataset on teenage marijuana use (see [6]):
## Superposition of structures N, L and M for the girls.
data(marijuana)
fit6 <- drm(y~age+cluster(id)+Time(age), data=marijuana,
subset=sex=="female", dep=list("NLM", ~kappa1==1,
~kappa2==0, ~tau12==1, ~tau21==1, ~tau11==tau22))
## Parameter restrictions with functions using M-structure for the boys.
## Plot the second order dependence ratios:
plot(depratio(y~cluster(id)+Time(age), data=marijuana,
subset=sex=="male"))
## fit the model in [6]:
fit7 <- drm(y~age+cluster(id)+Time(age), data=marijuana,
subset=sex=="male", dep=list("M",
tau12~function(a=1,b=0) a+b*c(0:3),
tau21~function(a=1,b=0) a+b*c(0:3)))
## Not run:
# ##############################################
# ## Covariates for the association (see [7]):
# ##############################################
# data(madras)
#
# ## plot empirical 2nd order dependence ratios with bootstrap CI's
# tau.madras <- depratio(symptom~cluster(id)+Time(month), data=madras,
# boot.ci = TRUE, n.boot = 1000)
# plot(tau.madras, log="y", ylim=c(1,40), plot.ci=TRUE)
#
# ## create matrix for profiles:
# W.madras <- profiles.drm(n.categories=2, n.repetitions=12, "M")
#
# ## create four-level covariate, combining age and sex:
# madras$age.sex <- factor(paste(madras$age,madras$sex,sep="."))
#
# ## fit the model in [7], Section 4:
# fit8 <- drm(symptom~age+sex+month+month:age+month:sex+cluster(id)+Time(month),
# data=madras, Ncond=FALSE, save.profiles=FALSE, pmatrix="W.madras",
# dep=list("NM",nu~nu:age.sex,
# tau~function(a0=0,a1=0) 1+a0*exp(a1*c(0:10))), print=2)
#
# ###################################################
# ## Dropout model on top of regression & association
# ###################################################
# ## Continue with the madras data.
# ## fit a model without the dropout model:
# fit9 <- drm(symptom~age+sex+month+month:age+month:sex+cluster(id)+Time(month),
# data=madras, save.profiles=FALSE, pmatrix="W.madras", print=0,
# dep=list("NM", tau~function(a0=0,a1=0) 1+a0*exp(a1*c(0:10))))
#
# ## A dropout model assuming MCAR for the thought disorders:
#
# mcar <- drm(symptom~age+sex+month+month:age+month:sex+cluster(id)+Time(month),
# data=madras, save.profiles=FALSE, pmatrix="W.madras",
# dep=list("NM", tau~function(a0=0,a1=0) 1+a0*exp(a1*c(0:10)),
# ~symptom.cur==0,~symptom.prev==0),
# dropout=TRUE, start=c(coef(fit9), -4))
#
# ## A dropout model assuming MAR; including sex as a covariate:
#
# mar <- drm(symptom~age+sex+month+month:age+month:sex+cluster(id)+Time(month),
# data=madras, save.profiles=FALSE, pmatrix="W.madras",
# dep=list("NM", tau~function(a0=0,a1=0) 1+a0*exp(a1*c(0:10)),
# ~symptom.cur==0), dropout=TRUE, drop.x=sex,
# start=c(coef(mcar),0,0))
#
# ## A dropout model assuming MNAR and sex as a covariate:
#
# mnar <- drm(symptom~age+sex+month+month:age+month:sex+cluster(id)+Time(month),
# data=madras, save.profiles=FALSE, pmatrix="W.madras",
# dep=list("NM", tau~function(a0=0,a1=0) 1+a0*exp(a1*c(0:10))),
# dropout=TRUE, drop.x=sex, start=c(coef(mcar),0,0,0))
#
# ## print out coefficients and std.errors:
# coef(summary(mnar))
# ## End(Not run)
## std.error of `symptom.cur' all over the place; too few dropouts
## for a comprehensive evaluation of the dropout mechanism
Run the code above in your browser using DataLab