simulation-multinom-vs-logistic-for-mary.R

library(nnet)

nsnp <- 3
simdata <- function(n) {
  X <- matrix(pmin(abs(rnorm(n*nsnp)),3),nrow=n)

  p <- rbinom(n, 1, X[,1]/3)
  q <- rbinom(n, 1, X[,nsnp]/3)

  Y <- ifelse(p>0, 1, ifelse(q>0, 2, 0))
  table(Y)
  
return(list(X=as.data.frame(X),Y=Y))
}
r.binom <- function(d) {
y <- pmin(d$Y, 1)
m1 <- glm(y ~ ., data=d$X, subset=d$Y %in% c(0,1), family="binomial")
m2 <- glm(y ~ ., data=d$X, subset=d$Y %in% c(0,2), family="binomial")
ret <- as.data.frame(rbind(summary(m1)$coefficients[-1,1:2],
                           summary(m2)$coefficients[-1,1:2]))
rownames(ret) <- NULL
colnames(ret) <- c("beta","se")
ret$snp<-1:nsnp
ret$var <- c(rep(1,nsnp),rep(2,nsnp))
ret$analysis <- "binom"
return(ret)
}
r.multinom <- function(d) {
m <- multinom(d$Y ~ ., data=d$X)
sm <- summary(m)
ret <- as.data.frame(rbind(cbind(sm$coefficients["1",-1],sm$standard.errors["1",-1]),
                           cbind(sm$coefficients["2",-1],sm$standard.errors["2",-1])))
rownames(ret) <- NULL
colnames(ret) <- c("beta","se")
ret$snp <- 1:nsnp
ret$var <- c(rep(1,nsnp),rep(2,nsnp))
ret$analysis <- "multinom"
return(ret)
}

runone <- function(n) {
  d <- simdata(n)
  bd <- r.binom(d)
  md <- r.multinom(d)
  ret <- rbind(bd,md)
  ret$n <- n
  return(ret)
}

results <- do.call("rbind",lapply(as.list(seq(100,2000,100)), runone))

library(ggplot2)
ggplot(results,aes(x=n,y=beta,ymin=beta-1.96*se,ymax=beta+1.96*se,col=as.factor(analysis))) +
  geom_pointrange() +
  geom_line() + facet_grid(snp ~ var)