# NOT RUN {
# Load data frame with (g, Y, X, Xtilde) values for 496 pools, list of C
# values for members of each pool, and list of Xtilde values where 25
# single-specimen pools have replicates. Xtilde values are affected by
# processing error and measurement error. True log-OR = 0.5, sigsq_p = 0.25,
# sigsq_m = 0.1.
data(dat_p_gdfa)
dat <- dat_p_gdfa$dat
reps <- dat_p_gdfa$reps
c.list <- dat_p_gdfa$c.list
# Unobservable truth estimator - use precise X's
fit.unobservable <- p_gdfa(
g = dat$g,
y = dat$y,
xtilde = dat$x,
c = c.list,
errors = "neither"
)
fit.unobservable$estimates
# Naive estimator - use imprecise Xtilde's, but treat as precise
fit.naive <- p_gdfa(
g = dat$g,
y = dat$y,
xtilde = dat$xtilde,
c = c.list,
errors = "neither"
)
fit.naive$estimates
# Corrected estimator - use Xtilde's and account for errors (not using
# replicates here)
# }
# NOT RUN {
fit.noreps <- p_gdfa(
g = dat$g,
y = dat$y,
xtilde = dat$xtilde,
c = c.list,
errors = "both"
)
fit.noreps$estimates
# Corrected estimator - use Xtilde's including 25 replicates
fit.reps <- p_gdfa(
g = dat$g,
y = dat$y,
xtilde = reps,
c = c.list,
errors = "both"
)
fit.reps$estimates
# Same as previous, but allowing for non-constant odds ratio.
fit.nonconstant <- p_gdfa(
g = dat$g,
y = dat$y,
xtilde = reps,
c = c.list,
constant_or = FALSE,
errors = "both",
hcubature_list = list(tol = 1e-4)
)
fit.nonconstant$estimates
# Visualize estimated log-OR vs. X based on previous model fit
p <- plot_gdfa(
estimates = fit.nonconstant$estimates,
varcov = fit.nonconstant$theta.var,
xrange = range(dat$xtilde[dat$g == 1]),
cvals = mean(unlist(c))
)
p
# }
# NOT RUN {
# }
Run the code above in your browser using DataLab