tests/testthat/DeprecatedTestCode.R
In schildjs/ods4lda: Outcome Dependent Sampling (ODS) Designs for Linear Mixed Models Data Analysis
## Run the code for one realization from simulation scenarion e) in Schildcrout et al (2015) with the exception that we
## estimate the parameter associated with C_i here. In the paper it was not estimated.
##
## Expect Fit.Biv to take substantially longer to fit than the others
setwd("~/rsch/OutDepSamp/Software/TestODS4LDA")
rm(list=ls())
library(ODS4LDA)
odsInt <- read.csv("odsInt.csv")
odsSlp <- read.csv("odsSlp.csv")
odsBiv <- read.csv("odsBiv.csv")
date()
Fit.Int <- acml.linear(y=odsInt$Y,
x=as.matrix(cbind(1, odsInt[,c("time","snp","snptime","confounder")])),
z=as.matrix(cbind(1, odsInt$time)),
id=odsInt$id,
InitVals=c(5, 1, -2.5, 0.75, 0, 1.6094379, 0.2231436, -0.5108256, 1.6094379),
ProfileCol=NA,
cutpoints=c(-2.569621, 9.992718),
SampProb=c(1, 0.1228, 1),
w.function="intercept")
date()
Fit.Int2 <- acml.lmem(y=odsInt$Y,
x=as.matrix(cbind(1, odsInt[,c("time","snp","snptime","confounder")])),
z=as.matrix(cbind(1, odsInt$time)),
id=odsInt$id,
InitVals=c(5, 1, -2.5, 0.75, 0, 1.6094379, 0.2231436, -0.5108256, 1.6094379),
ProfileCol=NA,
cutpoints=c(-2.569621, 9.992718),
SampProb=c(1, 0.1228, 1),
w.function="intercept")
date()
Fit.Slp <- acml.linear(y=odsSlp$Y,
x=as.matrix(cbind(1, odsSlp[,c("time","snp","snptime","confounder")])),
z=as.matrix(cbind(1, odsSlp$time)),
id=odsSlp$id,
InitVals=c(5, 1, -2.5, 0.75, 0, 1.6094379, 0.2231436, -0.5108256, 1.6094379),
ProfileCol=NA,
cutpoints=c(-0.7488912, 3.4557775),
SampProb=c(1, 0.1228, 1),
w.function="slope")
date()
Fit.Slp2 <- acml.lmem(y=odsSlp$Y,
x=as.matrix(cbind(1, odsSlp[,c("time","snp","snptime","confounder")])),
z=as.matrix(cbind(1, odsSlp$time)),
id=odsSlp$id,
InitVals=c(5, 1, -2.5, 0.75, 0, 1.6094379, 0.2231436, -0.5108256, 1.6094379),
ProfileCol=NA,
cutpoints=c(-0.7488912, 3.4557775),
SampProb=c(1, 0.1228, 1),
w.function="slope")
date()
Fit.Biv <- acml.linear(y=odsBiv$Y,
x=as.matrix(cbind(1, odsBiv[,c("time","snp","snptime","confounder")])),
z=as.matrix(cbind(1, odsBiv$time)),
id=odsBiv$id,
InitVals=c(5, 1, -2.5, 0.75, 0, 1.6094379, 0.2231436, -0.5108256, 1.6094379),
ProfileCol=NA,
cutpoints=c(-4.413225, 11.935188, -1.390172, 4.084768),
SampProb=c(0.122807, 1),
w.function="bivar")
date()
Fit.Biv2 <- acml.lmem(y=odsBiv$Y,
x=as.matrix(cbind(1, odsBiv[,c("time","snp","snptime","confounder")])),
z=as.matrix(cbind(1, odsBiv$time)),
id=odsBiv$id,
InitVals=c(5, 1, -2.5, 0.75, 0, 1.6094379, 0.2231436, -0.5108256, 1.6094379),
ProfileCol=NA,
cutpoints=c(-4.413225, 11.935188, -1.390172, 4.084768),
SampProb=c(0.122807, 1),
w.function="bivar")
date()
Fit.Int
Fit.Slp
Fit.Biv
library(MASS)
GenPopnData <- function(N = 1000, n = 11, beta = c(1, 0.25, 0, 0.25),
sig.b0 = 0.25, sig.b1 = 0.25, rho = 0, sig.e = 0.5,
pr.conf=0.25, pr.grp.0 = 0.1, pr.grp.1 = 0.2, n.low=11, n.high=11){
id <- rep(1:N, each=n)
conf.i <- rbinom(N, 1, pr.conf)
pr.grp.i <- ifelse(conf.i==1, pr.grp.1, pr.grp.0)
grp.i <- rbinom(N,1,pr.grp.i)
conf <- rep(conf.i, each=n)
grp <- rep(grp.i, each=n)
#grp <- rep(rbinom(N, 1, prob.grp), each=n)
time <- rep(seq(-1,1, length.out=n), N)
cov.mat <- matrix(c(sig.b0^2, rho*sig.b0*sig.b1, rho*sig.b0*sig.b1, sig.b1^2),2,2)
bi <- mvrnorm(N, mu=c(0,0), Sigma=cov.mat)
b <- cbind(rep(bi[,1], each=n),rep(bi[,2], each=n))
error <- rnorm(N*n, 0, sig.e)
X <- as.matrix(cbind(1, time, grp, time*grp, conf))
Z <- X[,c(1:2)]
Y <- X %*% beta + Z[,1]*b[,1] + Z[,2]*b[,2] + error
## Induce MCAR dropout
n.obs <- rep(sample(rep(c(n.low:n.high),2), N, replace=TRUE), each=n)
obs.num <- rep(c(1:n), N)
X <- X[obs.num<= n.obs,]
Z <- Z[obs.num<= n.obs,]
Y <- Y[obs.num<= n.obs]
id <- id[obs.num<= n.obs]
list(id=id, X=X, Y=Y, Z=Z, ni=c(unlist(tapply(Y,id,length))),
N=N, n=n, beta=beta,sig.b0=sig.b0, sig.b1=sig.b1, rho=rho, sig.e=sig.e, pr.grp.1=pr.grp.1, pr.grp.0=pr.grp.0)
}
LinRegFn <- function(data){ X <- cbind(1, data[,3])
Y <- data$Y
solve(t(X)%*%X) %*% t(X) %*% Y}
est.quants <- function(N, n, beta, sig.b0, sig.b1, rho, sig.e, pr.grp.0, pr.grp.1, pr.conf, quant, n.low=11, n.high=11){
d <- GenPopnData(N=N, n=n, beta=beta, sig.b0=sig.b0, sig.b1=sig.b1, rho=rho, sig.e=sig.e,
pr.conf=pr.conf, pr.grp.0 = pr.grp.0, pr.grp.1 = pr.grp.1)
data.tmp <- data.frame(id=d$id, Y=d$Y, X.time=d$X[,2])
out <- matrix(unlist(lapply(split(data.tmp, data.tmp$id), LinRegFn)), byrow=TRUE, ncol=2)
out <- cbind(c(tapply(data.tmp$Y, data.tmp$id, mean)), out)
r <- rbind( quant,
quantile(out[,1], probs=quant),
quantile(out[,2], probs=quant),
quantile(out[,3], probs=quant))
r}
## Do a search to find the quantiles that correspond to the central rectangle that contains 60 and 80 percent of
## the subject specific intercepts and slopes. This is not necessary if slope and intercepts are independent
## but with unequal followup they were positiviely correlated. Searches for the smallest 'rectangle' defined
## by quantiles that contains 60 and 80 percent of the data
est.bivar.lims <- function(N, n, beta, sig.b0, sig.b1, rho, sig.e, pr.grp.0, pr.grp.1, pr.conf, quants, n.low=11, n.high=11){
d <- GenPopnData(N=N, n=n, beta=beta, sig.b0=sig.b0, sig.b1=sig.b1, rho=rho, sig.e=sig.e,
pr.conf=pr.conf, pr.grp.0 = pr.grp.0, pr.grp.1 = pr.grp.1)
#d <- GenPopnData(N, n, beta, sig.b0, sig.b1, rho, sig.e, prob.grp, n.low=n.low, n.high=n.high)
data.tmp <- data.frame(id=d$id, Y=d$Y, X.time=d$X[,2])
out <- matrix(unlist(lapply(split(data.tmp, data.tmp$id), LinRegFn)), byrow=TRUE, ncol=2)
out <- cbind(c(tapply(data.tmp$Y, data.tmp$id, mean)), out)
print(cor(out[,2], out[,3]))
q1 <- .99
Del <- 1
while (Del>0.001){ q1 <- q1-.00025
Del <- abs(mean(out[,2] > quantile(out[,2], probs=1-q1) &
out[,2] < quantile(out[,2], probs=q1) &
out[,3] > quantile(out[,3], probs=1-q1) &
out[,3] < quantile(out[,3], probs=q1)) - quants[1])
#print(c(q1,Del))
q1}
q2 <- .99
Del <- 1
while (Del>0.001){ q2 <- q2-.00025
Del <- abs(mean(out[,2] > quantile(out[,2], probs=1-q2) &
out[,2] < quantile(out[,2], probs=q2) &
out[,3] > quantile(out[,3], probs=1-q2) &
out[,3] < quantile(out[,3], probs=q2)) - quants[2])
q2}
rbind( c( quantile(out[,2], probs=1-q1), quantile(out[,2], probs=q1), quantile(out[,3], probs=1-q1), quantile(out[,3], probs=q1)),
c( quantile(out[,2], probs=1-q2), quantile(out[,2], probs=q2), quantile(out[,3], probs=1-q2), quantile(out[,3], probs=q2)))
}
ODS.Sampling.Bivar <- function(dat, ## a list generated from the GenPopnData() function
PopnQuantsBivar, ## a matrix from est.bivar.lims function
PopnPropInRectangle, ## Proportion of subjects in the central rectangle
TargetNSampledPerStratum, ## Theoretical (and possibly observed) number sampled per stratum
SamplingStrategy){ ## Options are "IndepODS" and "DepODS"
NCohort <- length(unique(dat$id))
NStratumThry <- round(NCohort*c(PopnPropInRectangle, (1-PopnPropInRectangle)))
SampProbThry <- TargetNSampledPerStratum / NStratumThry
Lims <- (PopnPropInRectangle==.6)*PopnQuantsBivar[1,] + (PopnPropInRectangle==.8)*PopnQuantsBivar[2,]
uid <- unique(dat$id)
ni <- c(unlist(tapply(dat$id, dat$id, length)))
SampVar <- NULL
for(i in uid){ yi <- dat$Y[dat$id==i]
xi <- dat$X[dat$id==i,1:2]
SampVar <- rbind(SampVar, c(mean(yi), solve(t(xi)%*%xi) %*% t(xi) %*% yi))
}
print(sum(SampVar[,2]>Lims[1] & SampVar[,2]<Lims[2]))
print(sum(SampVar[,3]>Lims[3] & SampVar[,3]<Lims[4]))
SampStratum <- ifelse(SampVar[,2]>Lims[1] & SampVar[,2]<Lims[2] &
SampVar[,3]>Lims[3] & SampVar[,3]<Lims[4], 1,2)
print(table(SampStratum))
NperStratum <- unlist(tapply(uid, SampStratum, length))
SampProbiThry <- ifelse(SampStratum==1, SampProbThry[1],
ifelse(SampStratum==2, SampProbThry[2], NA))
Sampled <- rbinom(length(SampProbiThry), 1, SampProbiThry)
SampProbObs <- c(tapply(Sampled, SampStratum, mean))
SampProbiObs <- ifelse(SampStratum==1, SampProbObs[1],
ifelse(SampStratum==2, SampProbObs[2], NA))
## Independent Sampling
if (SamplingStrategy=="IndepODS") InODS <- uid[ Sampled==1]
TheSample <- dat$id %in% InODS
X.ods <- dat$X[TheSample,]
Y.ods <- dat$Y[TheSample]
Z.ods <- dat$Z[TheSample,]
id.ods <- dat$id[TheSample]
if (SamplingStrategy=="IndepODS"){
SampProbi.ods <- rep(SampProbiThry, ni)[TheSample]
SampProb.ods <- SampProbThry
}
#if (SamplingStrategy=="DepODS"){
# SampProbi.ods <- rep(SampProbiObs, ni)[TheSample]
# SampProb.ods <- SampProbObs
#}
SampStrat.ods <- SampStratum[InODS]
Qi <- SampVar[InODS,2:3]
dup.id <- duplicated(dat$id)
dat.univariate <- dat$X[!dup.id,]
dat.univ.ods <- dat.univariate[InODS,]
list(#X=X.ods, Y=Y.ods, Z=Z.ods, id=id.ods,
X=dat$X, Y=dat$Y, Z=dat$Z, id=dat$id,
SampProb=SampProb.ods, SampProbi=SampProbi.ods,
N=dat$N, n=dat$n, beta=dat$beta, sig.b0=dat$sig.b0,
sig.b1=dat$sig.b1, rho=dat$rho, sig.e=dat$sig.e,
prob.grp=dat$prob.grp,w.function="bivar",
#cutpoint=c(C1,C2),
SampStratum=SampStrat.ods, Qi=Qi, dat.univ.ods=dat.univ.ods,
cutpoint=Lims, InSample=TheSample)
}
ODS.Sampling <- function(dat, ## a list generated from the GenPopnData() function
PopnQuants, ## a matrix from est.quants function
w.function, ## Response summary to sample on ("mean","intercept", or "slope")
quants, ## Population quantiles to define the theoretical sampling strata
TargetNSampledPerStratum, ## Theoretical (and possibly observed) number sampled per stratum
SamplingStrategy, ## Options are "IndepODS" and "DepODS"
Univariate=FALSE){
NCohort <- length(unique(dat$id))
NStratumThry <- round(NCohort*c(quants[1], quants[2]-quants[1], 1-quants[2]))
print(TargetNSampledPerStratum)
print(NStratumThry)
SampProbThry <- TargetNSampledPerStratum / NStratumThry
C1 <- ifelse(w.function=="mean", PopnQuants[2,match(quants[1], PopnQuants[1,])],
ifelse(w.function=="intercept", PopnQuants[3,match(quants[1], PopnQuants[1,])],
ifelse(w.function=="slope", PopnQuants[4,match(quants[1], PopnQuants[1,])])))
C2 <- ifelse(w.function=="mean", PopnQuants[2,match(quants[2], PopnQuants[1,])],
ifelse(w.function=="intercept", PopnQuants[3,match(quants[2], PopnQuants[1,])],
ifelse(w.function=="slope", PopnQuants[4,match(quants[2], PopnQuants[1,])])))
# C1 <- ifelse(quants[1]==.1 & w.function=="mean", PopnQuants[1,1],
# ifelse(quants[1]==.1 & w.function=="intercept", PopnQuants[2,1],
# ifelse(quants[1]==.1 & w.function=="slope", PopnQuants[3,1],
# ifelse(quants[1]==.125 & w.function=="mean", PopnQuants[1,2],
# ifelse(quants[1]==.125 & w.function=="intercept", PopnQuants[2,2],
# ifelse(quants[1]==.125 & w.function=="slope", PopnQuants[3,2],
# ifelse(quants[1]==.2 & w.function=="mean", PopnQuants[1,3],
# ifelse(quants[1]==.2 & w.function=="intercept", PopnQuants[2,3],
# ifelse(quants[1]==.2 & w.function=="slope", PopnQuants[3,3])))))))))
# C2 <- ifelse(quants[2]==.8 & w.function=="mean", PopnQuants[1,4],
# ifelse(quants[2]==.8 & w.function=="intercept", PopnQuants[2,4],
# ifelse(quants[2]==.8 & w.function=="slope", PopnQuants[3,4],
# ifelse(quants[2]==.875 & w.function=="mean", PopnQuants[1,5],
# ifelse(quants[2]==.875 & w.function=="intercept", PopnQuants[2,5],
# ifelse(quants[2]==.875 & w.function=="slope", PopnQuants[3,5],
# ifelse(quants[2]==.9 & w.function=="mean", PopnQuants[1,6],
# ifelse(quants[2]==.9 & w.function=="intercept", PopnQuants[2,6],
# ifelse(quants[2]==.9 & w.function=="slope", PopnQuants[3,6])))))))))
uid <- unique(dat$id)
ni <- c(unlist(tapply(dat$id, dat$id, length)))
SampVar <- NULL
for(i in uid){ yi <- dat$Y[dat$id==i]
xi <- dat$X[dat$id==i,1:2]
SampVar <- rbind(SampVar, c(mean(yi), solve(t(xi)%*%xi) %*% t(xi) %*% yi))
}
SampVar <- (w.function=="mean")*SampVar[,1] +
(w.function=="intercept")*SampVar[,2] +
(w.function=="slope")*SampVar[,3]
SampStratum <- ifelse(SampVar<C1, 1,
ifelse(SampVar<C2, 2, 3))
NperStratum <- unlist(tapply(uid, SampStratum, length))
SampProbiThry <- ifelse(SampVar<C1, SampProbThry[1],
ifelse(SampVar<C2, SampProbThry[2], SampProbThry[3]))
Sampled <- rbinom(length(SampProbiThry), 1, SampProbiThry)
SampProbObs <- c(tapply(Sampled, SampStratum, mean))
SampProbiObs <- ifelse(SampVar<C1, SampProbObs[1],
ifelse(SampVar<C2, SampProbObs[2], SampProbObs[3]))
#print(rbind(SampProbThry, SampProbObs))
#print(cbind(SampProbiThry,SampProbiObs))
## Independent Sampling
if (SamplingStrategy=="IndepODS") InODS <- uid[ Sampled==1]
TheSample <- dat$id %in% InODS
X.ods <- dat$X[TheSample,]
Y.ods <- dat$Y[TheSample]
Z.ods <- dat$Z[TheSample,]
id.ods <- dat$id[TheSample]
if (SamplingStrategy=="IndepODS"){
SampProbi.ods <- rep(SampProbiThry, ni)[TheSample]
SampProb.ods <- SampProbThry
}
SampStrat.ods <- SampStratum[InODS]
Qi <- SampVar[InODS]
dup.id <- duplicated(dat$id)
dat.univariate <- dat$X[!dup.id,]
dat.univ.ods <- dat.univariate[InODS,]
SampProb <- SampProbThry
SampProbi <- rep(SampProbiThry, each=ni)
Qi <- SampVar
list(X=dat$X, Y=dat$Y, Z=dat$Z, id=dat$id,
SampProb=SampProb, SampProbi=SampProbi,
N=dat$N, n=dat$n, beta=dat$beta, sig.b0=dat$sig.b0,
sig.b1=dat$sig.b1, rho=dat$rho, sig.e=dat$sig.e,
prob.grp=dat$prob.grp, w.function=w.function,
cutpoint=c(C1,C2), SampStratum=SampStratum, Qi=Qi, InSample=TheSample)
}
#########################################################
#########################################################
#########################################################
## note that we really do not need quants, PopnQuants, and w.function but to make the fitting function work
Random.Sampling <- function(d, quants=c(.1, .9), PopnQuants, w.function="mean", n=225){
s <- sample(unique(d$id), n)
TheSample <- d$id %in% s
X.rand <- d$X
Y.rand <- d$Y
Z.rand <- d$Z
id.rand <- d$id
C1 <- PopnQuants[2,match(quants[1], PopnQuants[1,])]
C2 <- PopnQuants[2,match(quants[2], PopnQuants[1,])]
list(X=X.rand, Y=Y.rand, Z=Z.rand, id=id.rand, InSample=TheSample,
N=d$N, n=d$n, beta=d$beta, sig.b0=d$sig.b0,
sig.b1=d$sig.b1, rho=d$rho, sig.e=d$sig.e,
prob.grp=d$prob.grp,
## output below is only needed for the fitting function to run. They are not used.
SampProb=c(1,1,1), cutpoint=c(C1,C2), SampProbi=rep(1, length(Y.rand)), w.function=w.function)
}
Fit <- acml.lmem(y=odsSlp$Y,
x=as.matrix(cbind(1, odsSlp[,c("time","snp","snptime","confounder")])),
z=as.matrix(cbind(1, odsSlp$time)),
id=odsSlp$id,
InitVals=c(5, 1, -2.5, 0.75, 0, 1.6094379, 0.2231436, -0.5108256, 1.6094379),
ProfileCol=NA,
cutpoints=c(-0.7488912, 3.4557775),
SampProb=c(1, 0.1228, 1),
w.function="slope")
date()
params <- Fit$Ests
vcovs <- Fit$covar
Xe <- x1.insample[,3]
Xo <- x1.insample[,5]
Xe.Xo.Seq1.mod <- glm(Xe ~ Xo, family=binomial)
print(paste(date(), "Start MI", sep="-----"))
Ests.mi <- Var.mi <- NULL
Fit <- ACML.LME(y=y.insample, x=x.insample, z=z.insample, id=id.insample, InitVals=inits, ProfileCol=ProfileCol.s,
cutpoints=CutPoint.s, SampProb=SampProb.s, w.function=w.function.s)
## Get estimates and variance-covariance matrix
if (is.na(ProfileCol.s)){ n.par <- length(Fit$Ests)
params <- Fit$Ests
vcovs <- Fit$covar
for (m1 in 2:n.par){ for (m2 in 1:(m1-1)){ #print(c(m1,m2))
vcovs[m1,m2] <- vcovs[m2,m1] }}
}else{ n.par <- length(Fit$Ests)
params <- c(Fit$Ests[1:(ProfileCol.s-1)], 0, Fit$Ests[ProfileCol.s:n.par])
vcovs <- cbind(Fit$covar[,1:(ProfileCol.s-1) ], 0, Fit$covar[,ProfileCol.s:n.par])
vcovs <- rbind(vcovs[1:(ProfileCol.s-1), ], 0, vcovs[ProfileCol.s:n.par,])
for (m1 in 2:(n.par+length(ProfileCol.s))){ for (m2 in 1:(m1-1)){ #print(c(m1,m2))
vcovs[m1,m2] <- vcovs[m2,m1] }}
}
Xe <- x1.insample[,3]
Xo <- x1.insample[,5]
Xe.Xo.Seq1.mod <- glm(Xe ~ Xo, family=binomial)
print(paste(date(), "Start MI", sep="-----"))
Ests.mi <- Var.mi <- NULL
for (j in 1:n.imp){ #print(j)
## Draw a sample from the parameter estimate distribution for the model Xe | Xo, S=1
Xe.Xo.Seq1.MI.params <- rmvnorm(1, Xe.Xo.Seq1.mod$coef, summary(Xe.Xo.Seq1.mod)$cov.unscaled)
prXeEq1.Xo <- expit(cbind(1,x1.outsample[,5]) %*% t(Xe.Xo.Seq1.MI.params))
## Draw a sample from the parameter estimate distribution for the model Y | Xe, Xo, S=1
Y.Xe.Xo.Seq1.MI.params <- rmvnorm(1, params, vcovs)
## Condtl log likelihood if Xe=0: log(pr(y | Xo, Xe=0, S=1))
x.outsample.0 <- x.outsample
x.outsample.0[,3] <- 0
x.outsample.0[,4] <- x.outsample.0[,2] * x.outsample.0[,3]
prY.Xo.S1.Xe0 <- LogLikeAndScore(params= Y.Xe.Xo.Seq1.MI.params, y=y.outsample, x=x.outsample.0, z=z.outsample, id=id.outsample,
w.function=w.function.s, cutpoints=CutPoint.s, SampProb = SampProb.s, SampProbi=rep(1, length(y.outsample)),
##this w.function, cutpoints, SampProb, and SampProbi dont really factor into the calculation, but need tospecify them for the function to work
ProfileCol=ProfileCol.s, Keep.liC=TRUE)
## Condtl log likelihood if Xe=1: log(pr(y | Xo, Xe=1, S=1))
x.outsample.1 <- x.outsample
x.outsample.1[,3] <- 1
x.outsample.1[,4] <- x.outsample.1[,2] * x.outsample.1[,3]
prY.Xo.S1.Xe1 <- LogLikeAndScore(params= Y.Xe.Xo.Seq1.MI.params, y=y.outsample, x=x.outsample.1, z=z.outsample, id=id.outsample,
w.function=w.function.s, cutpoints=CutPoint.s, SampProb = SampProb.s, SampProbi=rep(1, length(y.outsample)),
##this w.function, cutpoints, SampProb, and SampProbi dont really factor into the calculation, but need tospecify them for the function to work
ProfileCol=ProfileCol.s, Keep.liC=TRUE)
odds <- exp(prY.Xo.S1.Xe1)*prXeEq1.Xo/(exp(prY.Xo.S1.Xe0)*(1-prXeEq1.Xo))
probs <- Odds2Prob(odds)
Impute.Xe.1 <- rbinom(length(probs), 1, probs)
y.mi <- Sampled.InAndOut$Y
x.mi <- Sampled.InAndOut$X
z.mi <- Sampled.InAndOut$Z
id.mi <- Sampled.InAndOut$id
for (j in id1.outsample){ x.mi[id.mi == j,3] <- Impute.Xe.1[id1.outsample==j] }
x.mi[,4] <- x.mi[,2]*x.mi[,3]
lme.mi <- lmer(y.mi ~ x.mi[,2]+x.mi[,3]+x.mi[,4]+x.mi[,5]+(z.mi[,2] | id.mi), REML=FALSE)
Ests.mi <- rbind(Ests.mi, lme.mi@fixef )
Var.mi <- rbind(Var.mi,diag(vcov(lme.mi)) )
# ACML.mi <- ACML.LME(y=y.mi, x=x.mi, z=z.mi, id=id.mi, InitVals=ACML$Ests)
# Ests.mi <- rbind(Ests.mi, ACML.mi$Ests)
# Var.mi <- rbind(Var.mi, diag(ACML.mi$covar))
}
schildjs/ods4lda documentation built on March 16, 2020, 8:16 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
## Run the code for one realization from simulation scenarion e) in Schildcrout et al (2015) with the exception that we
## estimate the parameter associated with C_i here. In the paper it was not estimated.
##
## Expect Fit.Biv to take substantially longer to fit than the others
setwd("~/rsch/OutDepSamp/Software/TestODS4LDA")
rm(list=ls())
library(ODS4LDA)
odsInt <- read.csv("odsInt.csv")
odsSlp <- read.csv("odsSlp.csv")
odsBiv <- read.csv("odsBiv.csv")
date()
Fit.Int <- acml.linear(y=odsInt$Y,
x=as.matrix(cbind(1, odsInt[,c("time","snp","snptime","confounder")])),
z=as.matrix(cbind(1, odsInt$time)),
id=odsInt$id,
InitVals=c(5, 1, -2.5, 0.75, 0, 1.6094379, 0.2231436, -0.5108256, 1.6094379),
ProfileCol=NA,
cutpoints=c(-2.569621, 9.992718),
SampProb=c(1, 0.1228, 1),
w.function="intercept")
date()
Fit.Int2 <- acml.lmem(y=odsInt$Y,
x=as.matrix(cbind(1, odsInt[,c("time","snp","snptime","confounder")])),
z=as.matrix(cbind(1, odsInt$time)),
id=odsInt$id,
InitVals=c(5, 1, -2.5, 0.75, 0, 1.6094379, 0.2231436, -0.5108256, 1.6094379),
ProfileCol=NA,
cutpoints=c(-2.569621, 9.992718),
SampProb=c(1, 0.1228, 1),
w.function="intercept")
date()
Fit.Slp <- acml.linear(y=odsSlp$Y,
x=as.matrix(cbind(1, odsSlp[,c("time","snp","snptime","confounder")])),
z=as.matrix(cbind(1, odsSlp$time)),
id=odsSlp$id,
InitVals=c(5, 1, -2.5, 0.75, 0, 1.6094379, 0.2231436, -0.5108256, 1.6094379),
ProfileCol=NA,
cutpoints=c(-0.7488912, 3.4557775),
SampProb=c(1, 0.1228, 1),
w.function="slope")
date()
Fit.Slp2 <- acml.lmem(y=odsSlp$Y,
x=as.matrix(cbind(1, odsSlp[,c("time","snp","snptime","confounder")])),
z=as.matrix(cbind(1, odsSlp$time)),
id=odsSlp$id,
InitVals=c(5, 1, -2.5, 0.75, 0, 1.6094379, 0.2231436, -0.5108256, 1.6094379),
ProfileCol=NA,
cutpoints=c(-0.7488912, 3.4557775),
SampProb=c(1, 0.1228, 1),
w.function="slope")
date()
Fit.Biv <- acml.linear(y=odsBiv$Y,
x=as.matrix(cbind(1, odsBiv[,c("time","snp","snptime","confounder")])),
z=as.matrix(cbind(1, odsBiv$time)),
id=odsBiv$id,
InitVals=c(5, 1, -2.5, 0.75, 0, 1.6094379, 0.2231436, -0.5108256, 1.6094379),
ProfileCol=NA,
cutpoints=c(-4.413225, 11.935188, -1.390172, 4.084768),
SampProb=c(0.122807, 1),
w.function="bivar")
date()
Fit.Biv2 <- acml.lmem(y=odsBiv$Y,
x=as.matrix(cbind(1, odsBiv[,c("time","snp","snptime","confounder")])),
z=as.matrix(cbind(1, odsBiv$time)),
id=odsBiv$id,
InitVals=c(5, 1, -2.5, 0.75, 0, 1.6094379, 0.2231436, -0.5108256, 1.6094379),
ProfileCol=NA,
cutpoints=c(-4.413225, 11.935188, -1.390172, 4.084768),
SampProb=c(0.122807, 1),
w.function="bivar")
date()
Fit.Int
Fit.Slp
Fit.Biv
library(MASS)
GenPopnData <- function(N = 1000, n = 11, beta = c(1, 0.25, 0, 0.25),
sig.b0 = 0.25, sig.b1 = 0.25, rho = 0, sig.e = 0.5,
pr.conf=0.25, pr.grp.0 = 0.1, pr.grp.1 = 0.2, n.low=11, n.high=11){
id <- rep(1:N, each=n)
conf.i <- rbinom(N, 1, pr.conf)
pr.grp.i <- ifelse(conf.i==1, pr.grp.1, pr.grp.0)
grp.i <- rbinom(N,1,pr.grp.i)
conf <- rep(conf.i, each=n)
grp <- rep(grp.i, each=n)
#grp <- rep(rbinom(N, 1, prob.grp), each=n)
time <- rep(seq(-1,1, length.out=n), N)
cov.mat <- matrix(c(sig.b0^2, rho*sig.b0*sig.b1, rho*sig.b0*sig.b1, sig.b1^2),2,2)
bi <- mvrnorm(N, mu=c(0,0), Sigma=cov.mat)
b <- cbind(rep(bi[,1], each=n),rep(bi[,2], each=n))
error <- rnorm(N*n, 0, sig.e)
X <- as.matrix(cbind(1, time, grp, time*grp, conf))
Z <- X[,c(1:2)]
Y <- X %*% beta + Z[,1]*b[,1] + Z[,2]*b[,2] + error
## Induce MCAR dropout
n.obs <- rep(sample(rep(c(n.low:n.high),2), N, replace=TRUE), each=n)
obs.num <- rep(c(1:n), N)
X <- X[obs.num<= n.obs,]
Z <- Z[obs.num<= n.obs,]
Y <- Y[obs.num<= n.obs]
id <- id[obs.num<= n.obs]
list(id=id, X=X, Y=Y, Z=Z, ni=c(unlist(tapply(Y,id,length))),
N=N, n=n, beta=beta,sig.b0=sig.b0, sig.b1=sig.b1, rho=rho, sig.e=sig.e, pr.grp.1=pr.grp.1, pr.grp.0=pr.grp.0)
}
LinRegFn <- function(data){ X <- cbind(1, data[,3])
Y <- data$Y
solve(t(X)%*%X) %*% t(X) %*% Y}
est.quants <- function(N, n, beta, sig.b0, sig.b1, rho, sig.e, pr.grp.0, pr.grp.1, pr.conf, quant, n.low=11, n.high=11){
d <- GenPopnData(N=N, n=n, beta=beta, sig.b0=sig.b0, sig.b1=sig.b1, rho=rho, sig.e=sig.e,
pr.conf=pr.conf, pr.grp.0 = pr.grp.0, pr.grp.1 = pr.grp.1)
data.tmp <- data.frame(id=d$id, Y=d$Y, X.time=d$X[,2])
out <- matrix(unlist(lapply(split(data.tmp, data.tmp$id), LinRegFn)), byrow=TRUE, ncol=2)
out <- cbind(c(tapply(data.tmp$Y, data.tmp$id, mean)), out)
r <- rbind( quant,
quantile(out[,1], probs=quant),
quantile(out[,2], probs=quant),
quantile(out[,3], probs=quant))
r}
## Do a search to find the quantiles that correspond to the central rectangle that contains 60 and 80 percent of
## the subject specific intercepts and slopes. This is not necessary if slope and intercepts are independent
## but with unequal followup they were positiviely correlated. Searches for the smallest 'rectangle' defined
## by quantiles that contains 60 and 80 percent of the data
est.bivar.lims <- function(N, n, beta, sig.b0, sig.b1, rho, sig.e, pr.grp.0, pr.grp.1, pr.conf, quants, n.low=11, n.high=11){
d <- GenPopnData(N=N, n=n, beta=beta, sig.b0=sig.b0, sig.b1=sig.b1, rho=rho, sig.e=sig.e,
pr.conf=pr.conf, pr.grp.0 = pr.grp.0, pr.grp.1 = pr.grp.1)
#d <- GenPopnData(N, n, beta, sig.b0, sig.b1, rho, sig.e, prob.grp, n.low=n.low, n.high=n.high)
data.tmp <- data.frame(id=d$id, Y=d$Y, X.time=d$X[,2])
out <- matrix(unlist(lapply(split(data.tmp, data.tmp$id), LinRegFn)), byrow=TRUE, ncol=2)
out <- cbind(c(tapply(data.tmp$Y, data.tmp$id, mean)), out)
print(cor(out[,2], out[,3]))
q1 <- .99
Del <- 1
while (Del>0.001){ q1 <- q1-.00025
Del <- abs(mean(out[,2] > quantile(out[,2], probs=1-q1) &
out[,2] < quantile(out[,2], probs=q1) &
out[,3] > quantile(out[,3], probs=1-q1) &
out[,3] < quantile(out[,3], probs=q1)) - quants[1])
#print(c(q1,Del))
q1}
q2 <- .99
Del <- 1
while (Del>0.001){ q2 <- q2-.00025
Del <- abs(mean(out[,2] > quantile(out[,2], probs=1-q2) &
out[,2] < quantile(out[,2], probs=q2) &
out[,3] > quantile(out[,3], probs=1-q2) &
out[,3] < quantile(out[,3], probs=q2)) - quants[2])
q2}
rbind( c( quantile(out[,2], probs=1-q1), quantile(out[,2], probs=q1), quantile(out[,3], probs=1-q1), quantile(out[,3], probs=q1)),
c( quantile(out[,2], probs=1-q2), quantile(out[,2], probs=q2), quantile(out[,3], probs=1-q2), quantile(out[,3], probs=q2)))
}
ODS.Sampling.Bivar <- function(dat, ## a list generated from the GenPopnData() function
PopnQuantsBivar, ## a matrix from est.bivar.lims function
PopnPropInRectangle, ## Proportion of subjects in the central rectangle
TargetNSampledPerStratum, ## Theoretical (and possibly observed) number sampled per stratum
SamplingStrategy){ ## Options are "IndepODS" and "DepODS"
NCohort <- length(unique(dat$id))
NStratumThry <- round(NCohort*c(PopnPropInRectangle, (1-PopnPropInRectangle)))
SampProbThry <- TargetNSampledPerStratum / NStratumThry
Lims <- (PopnPropInRectangle==.6)*PopnQuantsBivar[1,] + (PopnPropInRectangle==.8)*PopnQuantsBivar[2,]
uid <- unique(dat$id)
ni <- c(unlist(tapply(dat$id, dat$id, length)))
SampVar <- NULL
for(i in uid){ yi <- dat$Y[dat$id==i]
xi <- dat$X[dat$id==i,1:2]
SampVar <- rbind(SampVar, c(mean(yi), solve(t(xi)%*%xi) %*% t(xi) %*% yi))
}
print(sum(SampVar[,2]>Lims[1] & SampVar[,2]<Lims[2]))
print(sum(SampVar[,3]>Lims[3] & SampVar[,3]<Lims[4]))
SampStratum <- ifelse(SampVar[,2]>Lims[1] & SampVar[,2]<Lims[2] &
SampVar[,3]>Lims[3] & SampVar[,3]<Lims[4], 1,2)
print(table(SampStratum))
NperStratum <- unlist(tapply(uid, SampStratum, length))
SampProbiThry <- ifelse(SampStratum==1, SampProbThry[1],
ifelse(SampStratum==2, SampProbThry[2], NA))
Sampled <- rbinom(length(SampProbiThry), 1, SampProbiThry)
SampProbObs <- c(tapply(Sampled, SampStratum, mean))
SampProbiObs <- ifelse(SampStratum==1, SampProbObs[1],
ifelse(SampStratum==2, SampProbObs[2], NA))
## Independent Sampling
if (SamplingStrategy=="IndepODS") InODS <- uid[ Sampled==1]
TheSample <- dat$id %in% InODS
X.ods <- dat$X[TheSample,]
Y.ods <- dat$Y[TheSample]
Z.ods <- dat$Z[TheSample,]
id.ods <- dat$id[TheSample]
if (SamplingStrategy=="IndepODS"){
SampProbi.ods <- rep(SampProbiThry, ni)[TheSample]
SampProb.ods <- SampProbThry
}
#if (SamplingStrategy=="DepODS"){
# SampProbi.ods <- rep(SampProbiObs, ni)[TheSample]
# SampProb.ods <- SampProbObs
#}
SampStrat.ods <- SampStratum[InODS]
Qi <- SampVar[InODS,2:3]
dup.id <- duplicated(dat$id)
dat.univariate <- dat$X[!dup.id,]
dat.univ.ods <- dat.univariate[InODS,]
list(#X=X.ods, Y=Y.ods, Z=Z.ods, id=id.ods,
X=dat$X, Y=dat$Y, Z=dat$Z, id=dat$id,
SampProb=SampProb.ods, SampProbi=SampProbi.ods,
N=dat$N, n=dat$n, beta=dat$beta, sig.b0=dat$sig.b0,
sig.b1=dat$sig.b1, rho=dat$rho, sig.e=dat$sig.e,
prob.grp=dat$prob.grp,w.function="bivar",
#cutpoint=c(C1,C2),
SampStratum=SampStrat.ods, Qi=Qi, dat.univ.ods=dat.univ.ods,
cutpoint=Lims, InSample=TheSample)
}
ODS.Sampling <- function(dat, ## a list generated from the GenPopnData() function
PopnQuants, ## a matrix from est.quants function
w.function, ## Response summary to sample on ("mean","intercept", or "slope")
quants, ## Population quantiles to define the theoretical sampling strata
TargetNSampledPerStratum, ## Theoretical (and possibly observed) number sampled per stratum
SamplingStrategy, ## Options are "IndepODS" and "DepODS"
Univariate=FALSE){
NCohort <- length(unique(dat$id))
NStratumThry <- round(NCohort*c(quants[1], quants[2]-quants[1], 1-quants[2]))
print(TargetNSampledPerStratum)
print(NStratumThry)
SampProbThry <- TargetNSampledPerStratum / NStratumThry
C1 <- ifelse(w.function=="mean", PopnQuants[2,match(quants[1], PopnQuants[1,])],
ifelse(w.function=="intercept", PopnQuants[3,match(quants[1], PopnQuants[1,])],
ifelse(w.function=="slope", PopnQuants[4,match(quants[1], PopnQuants[1,])])))
C2 <- ifelse(w.function=="mean", PopnQuants[2,match(quants[2], PopnQuants[1,])],
ifelse(w.function=="intercept", PopnQuants[3,match(quants[2], PopnQuants[1,])],
ifelse(w.function=="slope", PopnQuants[4,match(quants[2], PopnQuants[1,])])))
# C1 <- ifelse(quants[1]==.1 & w.function=="mean", PopnQuants[1,1],
# ifelse(quants[1]==.1 & w.function=="intercept", PopnQuants[2,1],
# ifelse(quants[1]==.1 & w.function=="slope", PopnQuants[3,1],
# ifelse(quants[1]==.125 & w.function=="mean", PopnQuants[1,2],
# ifelse(quants[1]==.125 & w.function=="intercept", PopnQuants[2,2],
# ifelse(quants[1]==.125 & w.function=="slope", PopnQuants[3,2],
# ifelse(quants[1]==.2 & w.function=="mean", PopnQuants[1,3],
# ifelse(quants[1]==.2 & w.function=="intercept", PopnQuants[2,3],
# ifelse(quants[1]==.2 & w.function=="slope", PopnQuants[3,3])))))))))
# C2 <- ifelse(quants[2]==.8 & w.function=="mean", PopnQuants[1,4],
# ifelse(quants[2]==.8 & w.function=="intercept", PopnQuants[2,4],
# ifelse(quants[2]==.8 & w.function=="slope", PopnQuants[3,4],
# ifelse(quants[2]==.875 & w.function=="mean", PopnQuants[1,5],
# ifelse(quants[2]==.875 & w.function=="intercept", PopnQuants[2,5],
# ifelse(quants[2]==.875 & w.function=="slope", PopnQuants[3,5],
# ifelse(quants[2]==.9 & w.function=="mean", PopnQuants[1,6],
# ifelse(quants[2]==.9 & w.function=="intercept", PopnQuants[2,6],
# ifelse(quants[2]==.9 & w.function=="slope", PopnQuants[3,6])))))))))
uid <- unique(dat$id)
ni <- c(unlist(tapply(dat$id, dat$id, length)))
SampVar <- NULL
for(i in uid){ yi <- dat$Y[dat$id==i]
xi <- dat$X[dat$id==i,1:2]
SampVar <- rbind(SampVar, c(mean(yi), solve(t(xi)%*%xi) %*% t(xi) %*% yi))
}
SampVar <- (w.function=="mean")*SampVar[,1] +
(w.function=="intercept")*SampVar[,2] +
(w.function=="slope")*SampVar[,3]
SampStratum <- ifelse(SampVar<C1, 1,
ifelse(SampVar<C2, 2, 3))
NperStratum <- unlist(tapply(uid, SampStratum, length))
SampProbiThry <- ifelse(SampVar<C1, SampProbThry[1],
ifelse(SampVar<C2, SampProbThry[2], SampProbThry[3]))
Sampled <- rbinom(length(SampProbiThry), 1, SampProbiThry)
SampProbObs <- c(tapply(Sampled, SampStratum, mean))
SampProbiObs <- ifelse(SampVar<C1, SampProbObs[1],
ifelse(SampVar<C2, SampProbObs[2], SampProbObs[3]))
#print(rbind(SampProbThry, SampProbObs))
#print(cbind(SampProbiThry,SampProbiObs))
## Independent Sampling
if (SamplingStrategy=="IndepODS") InODS <- uid[ Sampled==1]
TheSample <- dat$id %in% InODS
X.ods <- dat$X[TheSample,]
Y.ods <- dat$Y[TheSample]
Z.ods <- dat$Z[TheSample,]
id.ods <- dat$id[TheSample]
if (SamplingStrategy=="IndepODS"){
SampProbi.ods <- rep(SampProbiThry, ni)[TheSample]
SampProb.ods <- SampProbThry
}
SampStrat.ods <- SampStratum[InODS]
Qi <- SampVar[InODS]
dup.id <- duplicated(dat$id)
dat.univariate <- dat$X[!dup.id,]
dat.univ.ods <- dat.univariate[InODS,]
SampProb <- SampProbThry
SampProbi <- rep(SampProbiThry, each=ni)
Qi <- SampVar
list(X=dat$X, Y=dat$Y, Z=dat$Z, id=dat$id,
SampProb=SampProb, SampProbi=SampProbi,
N=dat$N, n=dat$n, beta=dat$beta, sig.b0=dat$sig.b0,
sig.b1=dat$sig.b1, rho=dat$rho, sig.e=dat$sig.e,
prob.grp=dat$prob.grp, w.function=w.function,
cutpoint=c(C1,C2), SampStratum=SampStratum, Qi=Qi, InSample=TheSample)
}
#########################################################
#########################################################
#########################################################
## note that we really do not need quants, PopnQuants, and w.function but to make the fitting function work
Random.Sampling <- function(d, quants=c(.1, .9), PopnQuants, w.function="mean", n=225){
s <- sample(unique(d$id), n)
TheSample <- d$id %in% s
X.rand <- d$X
Y.rand <- d$Y
Z.rand <- d$Z
id.rand <- d$id
C1 <- PopnQuants[2,match(quants[1], PopnQuants[1,])]
C2 <- PopnQuants[2,match(quants[2], PopnQuants[1,])]
list(X=X.rand, Y=Y.rand, Z=Z.rand, id=id.rand, InSample=TheSample,
N=d$N, n=d$n, beta=d$beta, sig.b0=d$sig.b0,
sig.b1=d$sig.b1, rho=d$rho, sig.e=d$sig.e,
prob.grp=d$prob.grp,
## output below is only needed for the fitting function to run. They are not used.
SampProb=c(1,1,1), cutpoint=c(C1,C2), SampProbi=rep(1, length(Y.rand)), w.function=w.function)
}
Fit <- acml.lmem(y=odsSlp$Y,
x=as.matrix(cbind(1, odsSlp[,c("time","snp","snptime","confounder")])),
z=as.matrix(cbind(1, odsSlp$time)),
id=odsSlp$id,
InitVals=c(5, 1, -2.5, 0.75, 0, 1.6094379, 0.2231436, -0.5108256, 1.6094379),
ProfileCol=NA,
cutpoints=c(-0.7488912, 3.4557775),
SampProb=c(1, 0.1228, 1),
w.function="slope")
date()
params <- Fit$Ests
vcovs <- Fit$covar
Xe <- x1.insample[,3]
Xo <- x1.insample[,5]
Xe.Xo.Seq1.mod <- glm(Xe ~ Xo, family=binomial)
print(paste(date(), "Start MI", sep="-----"))
Ests.mi <- Var.mi <- NULL
Fit <- ACML.LME(y=y.insample, x=x.insample, z=z.insample, id=id.insample, InitVals=inits, ProfileCol=ProfileCol.s,
cutpoints=CutPoint.s, SampProb=SampProb.s, w.function=w.function.s)
## Get estimates and variance-covariance matrix
if (is.na(ProfileCol.s)){ n.par <- length(Fit$Ests)
params <- Fit$Ests
vcovs <- Fit$covar
for (m1 in 2:n.par){ for (m2 in 1:(m1-1)){ #print(c(m1,m2))
vcovs[m1,m2] <- vcovs[m2,m1] }}
}else{ n.par <- length(Fit$Ests)
params <- c(Fit$Ests[1:(ProfileCol.s-1)], 0, Fit$Ests[ProfileCol.s:n.par])
vcovs <- cbind(Fit$covar[,1:(ProfileCol.s-1) ], 0, Fit$covar[,ProfileCol.s:n.par])
vcovs <- rbind(vcovs[1:(ProfileCol.s-1), ], 0, vcovs[ProfileCol.s:n.par,])
for (m1 in 2:(n.par+length(ProfileCol.s))){ for (m2 in 1:(m1-1)){ #print(c(m1,m2))
vcovs[m1,m2] <- vcovs[m2,m1] }}
}
Xe <- x1.insample[,3]
Xo <- x1.insample[,5]
Xe.Xo.Seq1.mod <- glm(Xe ~ Xo, family=binomial)
print(paste(date(), "Start MI", sep="-----"))
Ests.mi <- Var.mi <- NULL
for (j in 1:n.imp){ #print(j)
## Draw a sample from the parameter estimate distribution for the model Xe | Xo, S=1
Xe.Xo.Seq1.MI.params <- rmvnorm(1, Xe.Xo.Seq1.mod$coef, summary(Xe.Xo.Seq1.mod)$cov.unscaled)
prXeEq1.Xo <- expit(cbind(1,x1.outsample[,5]) %*% t(Xe.Xo.Seq1.MI.params))
## Draw a sample from the parameter estimate distribution for the model Y | Xe, Xo, S=1
Y.Xe.Xo.Seq1.MI.params <- rmvnorm(1, params, vcovs)
## Condtl log likelihood if Xe=0: log(pr(y | Xo, Xe=0, S=1))
x.outsample.0 <- x.outsample
x.outsample.0[,3] <- 0
x.outsample.0[,4] <- x.outsample.0[,2] * x.outsample.0[,3]
prY.Xo.S1.Xe0 <- LogLikeAndScore(params= Y.Xe.Xo.Seq1.MI.params, y=y.outsample, x=x.outsample.0, z=z.outsample, id=id.outsample,
w.function=w.function.s, cutpoints=CutPoint.s, SampProb = SampProb.s, SampProbi=rep(1, length(y.outsample)),
##this w.function, cutpoints, SampProb, and SampProbi dont really factor into the calculation, but need tospecify them for the function to work
ProfileCol=ProfileCol.s, Keep.liC=TRUE)
## Condtl log likelihood if Xe=1: log(pr(y | Xo, Xe=1, S=1))
x.outsample.1 <- x.outsample
x.outsample.1[,3] <- 1
x.outsample.1[,4] <- x.outsample.1[,2] * x.outsample.1[,3]
prY.Xo.S1.Xe1 <- LogLikeAndScore(params= Y.Xe.Xo.Seq1.MI.params, y=y.outsample, x=x.outsample.1, z=z.outsample, id=id.outsample,
w.function=w.function.s, cutpoints=CutPoint.s, SampProb = SampProb.s, SampProbi=rep(1, length(y.outsample)),
##this w.function, cutpoints, SampProb, and SampProbi dont really factor into the calculation, but need tospecify them for the function to work
ProfileCol=ProfileCol.s, Keep.liC=TRUE)
odds <- exp(prY.Xo.S1.Xe1)*prXeEq1.Xo/(exp(prY.Xo.S1.Xe0)*(1-prXeEq1.Xo))
probs <- Odds2Prob(odds)
Impute.Xe.1 <- rbinom(length(probs), 1, probs)
y.mi <- Sampled.InAndOut$Y
x.mi <- Sampled.InAndOut$X
z.mi <- Sampled.InAndOut$Z
id.mi <- Sampled.InAndOut$id
for (j in id1.outsample){ x.mi[id.mi == j,3] <- Impute.Xe.1[id1.outsample==j] }
x.mi[,4] <- x.mi[,2]*x.mi[,3]
lme.mi <- lmer(y.mi ~ x.mi[,2]+x.mi[,3]+x.mi[,4]+x.mi[,5]+(z.mi[,2] | id.mi), REML=FALSE)
Ests.mi <- rbind(Ests.mi, lme.mi@fixef )
Var.mi <- rbind(Var.mi,diag(vcov(lme.mi)) )
# ACML.mi <- ACML.LME(y=y.mi, x=x.mi, z=z.mi, id=id.mi, InitVals=ACML$Ests)
# Ests.mi <- rbind(Ests.mi, ACML.mi$Ests)
# Var.mi <- rbind(Var.mi, diag(ACML.mi$covar))
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.