# simulate SBP data
simmeddat <- function(mu=144,bage=0.5,bsex=4.,bg=2.,pB=0.3,rvar=21^2,N=1000) {
ageb <- c(25,74)
pmale <- .5
htthresh <- 160
trprob <- .5
mutreff <- (-15.)
trvar <- 4^2
age <- runif(N,min=ageb[1],max=ageb[2])
sex <- 1*(runif(N)<=pmale)
gt <- rbinom(N,size=2,prob=pB)
y.true <- rnorm(N,mu,sqrt(rvar)) + bage*age + bsex*sex + bg*gt
d.true <- (y.true>=htthresh)
medication <- 1*d.true
medication[d.true] <- 1*(runif(sum(d.true))<=trprob)
treatm <- rnorm(sum(medication),mutreff,sqrt(trvar))
treatm[treatm>0] <- 0
treff <- rep(0,N)
treff[medication==1] <- treatm
y.obs <- y.true + treff
out <- data.frame(age,sex,gt,y.true,d.true,medication,treff,y.obs)
out
}
x <- simmeddat(bg=2.0,N=3000)
x[1:15,]
# substitute value of treated people
imptra <- npsubtreated(x$y.obs,x$medication)
imptra[1:15]
# Almost always, correlation should be higher for the "imputed" trait
cor(x$y.true,x$y.obs)
cor(x$y.true,imptra)
# see what comes out of regression
# analysis of true value
summary(lm(y.true~sex+age+gt,data=x))
# ignore treatment (as a rule, betas are underestimated;
# on average, power goes down)
summary(lm(y.obs~sex+age+gt,data=x))
# treatment as covariate (as a rule, betas are underestimated;
# on average, power goes down)
summary(lm(y.obs~sex+age+gt+medication,data=x))
# analyse "imputed" trait (as a rule betas are better; power
# approaches that of analysis of "true" trait)
summary(lm(imptra~sex+age+gt,data=x))
Run the code above in your browser using DataLab