sim2 <- function(nc, n, betas = runif(nc, -2, 2)) {
g <- rbinom(n, 2, 0.4)
m <- sapply(1:nc, function(i)
{
g * betas[i] + rnorm(n)
})
# for each individual sample cell type proportions
cellprop <- sapply(1:n, function(x) {a <- runif(nc); a/sum(a)}) %>% t()
# weighted sum
M <- (scale(m) * cellprop) %>% rowSums
res <- lapply(1:nc, function(i)
{
o <- summary(lm(M ~ g * cellprop[,i]))
tibble(cell=i, b=o$coef[4,1], se=o$coef[4,2], pval=o$coef[4,4], method="interaction")
}) %>% bind_rows()
res2 <- lapply(1:nc, function(i)
{
o <- summary(lm(m[,i] ~ g))
tibble(cell=i, b=o$coef[2,1], se=o$coef[2,2], pval=o$coef[2,4], method="marginal")
}) %>% bind_rows()
betas <- tibble(cell=1:nc, b=betas, method="true")
cellprop <- as_tibble(cellprop)
m <- as_tibble(m)
names(m) <- paste0("m", 1:nc)
names(cellprop) <- paste0("cp", 1:nc)
cellprop$M <- M
cellprop$g <- g
cellprop <- bind_cols(cellprop, m)
return(cellprop)
}
ana <- function(d) {
nc <- length(grep("cp", names(d)))
# per cell estimates
o1 <- lapply(1:nc, \(x) {
f <- paste0("M ~ g * cp", x)
o <- summary(lm(f, data=d))$coef
n <- rownames(o)
o <- as_tibble(o)
names(o) <- c("b", "se", "tval", "pval")
o$term <- n
o$cell <- x
o$mod <- "interaction"
o
}) %>% bind_rows()
# cell int estimates
o2 <- lapply(1:nc, \(x) {
f <- paste0("m", x, " ~ g")
o <- summary(lm(f, data=d))$coef
n <- rownames(o)
o <- as_tibble(o)
names(o) <- c("b", "se", "tval", "pval")
o$term <- n
o$cell <- x
o$mod <- "marginal"
o
}) %>% bind_rows()
bind_rows(o1, o2)
}
d <- sim2(2, 1000)
