# NOT RUN {
mu.mod<-function(a,l,delta,height){
mu <- as.numeric(l%*%c(0.2,0.2,0.3,-0.1))+triangle(a-2.5,delta,height)+a*(-0.1*l[,1]+0.1*l[,3])
return(mu)
}
triangle <- function(a,delta,height){
y <- exp(-a^2/((delta/2)^2))*height
return(y)
}
set.seed(2000)
n <- 500
d <- 4
sigma <- 0.5
delta <- 1
height <- 0
arange<-c(0.01,4.99)
l <- matrix(rnorm(n*d),ncol=d)
colnames(l) <- paste("l",1:4,sep="")
logit.lambda <- as.numeric(l%*%c(0.1,0.1,-0.1,0.2))
lambda <- exp(logit.lambda)/(1+exp(logit.lambda))
a <- rbeta(n, shape1 = lambda, shape2 =1-lambda)*5
mu <- mu.mod(a,l,delta,height)
residual.list <- rnorm(n,mean=0,sd=sigma)
y <- mu+residual.list
## We use the oracal propensity score and outcome regression for illustration
pilist <- dbeta(a/5, shape1=lambda, shape2 = 1-lambda)/5
varpilist <- colMeans(matrix(dbeta(rep(a,each=n)/5,
shape1=rep(lambda,n),
shape2 = 1-rep(lambda,n))/5, nrow=n))
mulist <- mu
malist <-colMeans(matrix(mu.mod(rep(a,each=n),l[rep(1:n,n),],delta,height),nrow=n))
out <- drdrtest.base(y,a,pilist,varpilist,mulist,malist,arange)
# }
Run the code above in your browser using DataLab