sim_pop <- function(n, beta1, beta2, af, h2)
{
g <- rbinom(n, 2, af)
prs <- g * beta1
vg <- rnorm(n, 0, h2)
v <- rnorm(n, 0, beta2 * g)
ve <- rnorm(n, 0, sqrt(1 - var(vg) - var(v) - var(prs)))
y <- prs + v + vg + ve
return(tibble(
g, y
))
}
sim_mz <- function(n, beta1, beta2, af, h2)
{
g <- rbinom(n, 2, af)
prs <- g * beta1
vg <- rnorm(n, 0, h2)
v1 <- rnorm(n, 0, beta2 * g)
ve1 <- rnorm(n, 0, sqrt(1 - var(vg) - var(v1) - var(prs)))
y1 <- prs + v1 + vg + ve1
v2 <- rnorm(n, 0, beta2 * g)
ve2 <- rnorm(n, 0, sqrt(1 - var(vg) - var(v2) - var(prs)))
y2 <- prs + v2 + vg + ve2
return(tibble(
g, y1, y2
))
}
test_drm <- function(g, y)
{
y.i <- tapply(y, g, median, na.rm=T)
z.ij <- abs(y - y.i[g+1])
fast_assoc(z.ij, g) %>%
as_tibble() %>%
mutate(method="drm")
}
test_mz <- function(g, y1, y2)
{
yd1 <- abs(y1-y2)
r1 <- fast_assoc(yd1, g) %>%
as_tibble() %>%
mutate(method="mzdiff")
r1
}
test_drm_lme4 <- function(g, y1, y2)
{
dat <- tibble(g=c(g,g), y=c(y1,y2), twin=rep(1:length(y1), each=2))
y.i <- tapply(dat$y, dat$g, median, na.rm=T)
dat$z.ij <- abs(dat$y - y.i[dat$g+1])
out <- lmer(z.ij ~ g + (1 | twin), dat) %>%
summary() %>%
coef() %>%
{
tibble(ahat=.[1,1], bhat=.[2,1], se=.[2,2], fval=.[2,4]^2, pval=.[2,5], n=nrow(dat), method="drm_lmer")
}
return(out)
}