library(MASS)
j<- 500
p_z <- 6 ## dimension of Z at each time point
n_t <- 4 ## number of time points
alphas <- list()
gammas <- list()
z_para <- c(-1/p_z, -1/p_z, -1/p_z, -1/p_z, -0.5/p_z,-0.5/p_z, -0.5/p_z,
-0.5/p_z)
Z <- list()
beta = c(0.2, -0.3, -0.01, 0.02, 0.03, 0.04, rep(rep(0.02,p_z), n_t))
beta_T = -0.2
sd_z_x = 0.4
X = mvrnorm(j, mu=c(1,5,6,7,8), Sigma=diag(1,5))
TRT = rbinom(j, size = 1, prob = 0.5)
Y_constant <- beta[1]+(X%*%beta[2:6])
Y0 <- 0
Y1 <- 0
A <- A1 <- A0 <- matrix(NA, nrow = j, ncol = n_t)
gamma <- c(1,-.1,-0.05,0.05,0.05,.05)
A0[,1] <- rbinom(j, size = 1, prob = 1/(1+exp(-(gamma[1] +
(X %*% gamma[2:6])))))
A1[,1] <- rbinom(j, size = 1, prob = 1/(1+exp(-(gamma[1] +
(X %*% gamma[2:6])))))
A[,1] <- A1[,1]*TRT + A0[,1]*(1-TRT)
for(i in 2:n_t){
alphas[[i]] <- matrix(rep(c(2.3, -0.3, -0.01, 0.02, 0.03, 0.04, -0.4),
p_z),ncol=p_z)
gammas[[i]] <- c(1, -0.1, 0.2, 0.2, 0.2, 0.2, rep(z_para[i],p_z))
Z0 <- alphas[[i]][1,]+(X%*%alphas[[i]][2:6,]) + mvrnorm(j, mu = rep(0,p_z)
, Sigma = diag(sd_z_x,p_z))
Z1 <- alphas[[i]][1,]+(X%*%alphas[[i]][2:6,])+alphas[[i]][7,] +
mvrnorm(j, mu = rep(0,p_z), Sigma = diag(sd_z_x,p_z))
Z[[i]] <- Z1*TRT + Z0*(1-TRT)
Y0 <- (Y0 + Z0 %*% matrix(beta[ (7 + (i-1)*p_z):
(6+p_z*i)],ncol = 1) )[,1]
Y1 <- (Y1 + Z1 %*% matrix(beta[ (7 + (i-1)*p_z):
(6+p_z*i)],ncol = 1) )[,1]
A0[,i] <- rbinom(j, size = 1,
prob = 1/(1+exp(-(gammas[[i]][1]+
(X%*%gammas[[i]][2:6])+Z0%*%matrix(gammas[[i]][7:
(7+p_z-1)], ncol=1))[,1])))*A0[,i-1]
A1[,i] <- rbinom(j, size = 1,
prob = 1/(1+exp(-(gammas[[i]][1]+
(X%*%gammas[[i]][2:6])+Z1%*%matrix(gammas[[i]][7:
(7+p_z-1)], ncol=1))[,1])))*A1[,i-1]
A[,i] <- A1[,i]*TRT + A0[,i]*(1-TRT)
}
Y0 <- Y0 + rnorm(j, mean = 0, sd = 0.3) + Y_constant
Y1 <- Y1 + + beta_T + rnorm(j, mean = 0, sd = 0.3) + Y_constant
Y <- as.vector( Y1*TRT+Y0*(1-TRT))
for(i in 2:n_t){
Z[[i]][A[,(i-1)]==0,] <- NA
}
Z[[1]] <- matrix(NA, nrow=nrow(Z1),ncol=ncol(Z1))
Y[A[,n_t] == 0] <- NA
# estimate the treatment difference
fit <- est_S_Star_Plus_MethodA(X, A, Z, Y, TRT)
fit
# Calculate the true values
true1 <- mean(Y1[A1[,n_t]==1])
true1
true0 <- mean(Y0[A1[,n_t]==1])
true0
true_d = true1 - true0
true_d
Run the code above in your browser using DataLab