# NOT RUN {
# Full model: Y|X, B
# Reduced model: Y|X
# X,B follows normal distribution with mean zero, variance one and correlation 0.3
# Y|X, B follows N(-1-0.5X+0.5B, 1)
set.seed(2333)
n = 800
data.n = data.frame(matrix(ncol = 3, nrow = n))
colnames(data.n) = c('Y', 'X', 'B')
data.n[,c('X', 'B')] = MASS::mvrnorm(n, rep(0,2), diag(0.7,2)+0.3)
data.n$Y = rnorm(n, -1 - 0.5*data.n$X + 0.5*data.n$B, 1)
# Generate the beta estimates from the external reduced model:
# generate a data of size m from the full model first, then fit the reduced regression
# to obtain the beta estiamtes and the corresponsing estimated variance
m = 30000
data.m = data.frame(matrix(ncol = 3, nrow = m))
names(data.m) = c('Y', 'X', 'B')
data.m[,c('X', 'B')] = MASS::mvrnorm(m, rep(0,2), diag(0.7,2)+0.3)
data.m$Y = rnorm(m, -1 - 0.5*data.m$X + 0.5*data.m$B, 1)
#fit Y|X to obtain the external beta estimates, save the beta estiamtes and the
# corresponsing estimated variance
fit.E = lm(Y ~ X, data = data.m)
beta.E = coef(fit.E)
names(beta.E) = c('int', 'X')
V.E = vcov(fit.E)
#get full model estimate from direct regression using the internal data only
fit.gamma.I = lm(Y ~ X + B, data = data.n)
gamma.I = coef(fit.gamma.I)
#Get CML estimates
gamma.CML = fxnCC_LinReg(p=2,
q=3,
YInt=data.n$Y,
XInt=data.n[,'X'],
BInt=data.n[,'B'],
betaHatExt=beta.E,
gammaHatInt=gamma.I,
n=nrow(data.n),
tol=1e-8,
maxIter=400,
factor=1)[["gammaHat"]]
# }
Run the code above in your browser using DataLab