# This is R code used in Shepherd, Shaw, and Dodd (2012) to simulate a setting where there are data errors in some records, an audit is performed to find the error rate and magnitude, and then this information is used to correct estimates of the relationship between X, Z, and Y.

rm(list=ls())

library("mvtnorm")

audit.correct.cov<-function(ystar,w,z,v,x,y) {

N<-length(ystar)
nv<-sum(v)

if(!is.matrix(z)) z = matrix(z,ncol=1)
if((!is.matrix(w)) w = matrix(w,ncol=1)
if((!is.matrix(x)) x = matrix(x,ncol=1)

dim.z<-ncol(z)
dim.w<-ncol(w)

t<-matrix(NA,nrow=N,ncol=dim.w)
for (j in 1:ncol(w)) {
t[,j]<-ifelse(v==1,w[,j]-x[,j],NA)
x[,j]<-ifelse(v==1,x[,j],NA)
}
tstar<-ifelse(v==1,ystar-y,NA)
y<-ifelse(v==1,y,NA)

mu.w<-apply(w,2,mean)
sigma.ww<-cov(w,w)
mu.z<-apply(z,2,mean)
sigma.zz<-cov(z,z)
sigma.wz<-cov(w,z)
sigma.tt<-cov(t,t,use="complete.obs")
mu.ystar<-mean(ystar)
sigma.wystar<-cov(w,ystar)
sigma.zystar<-cov(z,ystar)
sigma.ttstar<-cov(t,tstar,use="complete.obs")
sigma.ztstar<-cov(z,tstar,use="complete.obs")
sigma.zt<-cov(z,t,use="complete.obs")
sigma.xtstar<-cov(x,tstar,use="complete.obs")
sigma.yt<-cov(y,t,use="complete.obs")


C<-matrix(NA,nrow=dim.w+dim.z,ncol=dim.w+dim.z)
C[1:dim.w,1:dim.w]<-sigma.ww-sigma.tt
C[1:dim.w,(dim.w+1):(dim.w+dim.z)]<-sigma.wz-sigma.zt
C[(dim.w+1):(dim.w+dim.z),1:dim.w]<-t(sigma.wz-sigma.zt)
C[(dim.w+1):(dim.w+dim.z),(dim.w+1):(dim.w+dim.z)]<-sigma.zz

D<-c(sigma.wystar-sigma.xtstar-sigma.yt-sigma.ttstar,sigma.zystar-sigma.ztstar)

C.inv<-solve(C)

est<-C.inv%*%D
C.inv%*%D


mu.tstar<-mean(tstar,na.rm=TRUE)
mu.t<-apply(t,2,mean,na.rm=TRUE)

psi.mu.tstar<-ifelse(v==1,tstar-mu.tstar,0)
psi.mu.t<-psi.mu.x<-matrix(NA,N,dim.w)
for (j in 1:dim.w) {
psi.mu.t[,j]<-ifelse(v==1,t[,j]-mu.t[j],0)
}

psi.mu.ystar<-ystar-mu.ystar

psi.mu.w<-matrix(NA,N,dim.w)
junk<-matrix(NA,nrow=N,ncol=(dim.w^2))
junk1<-matrix(NA,nrow=N,ncol=dim.w^2)
psi.sigma.wystar<-matrix(NA,nrow=N,ncol=dim.w)
for (j in 1:dim.w) {
psi.mu.w[,j]<-w[,j]-mu.w[j]
psi.sigma.wystar[,j]<-psi.mu.ystar*psi.mu.w[,j]-sigma.wystar[j,]
}
for (j in 1:dim.w) {
for (k in j:dim.w) {
junk[,k+dim.w*(j-1)]<-psi.mu.w[,j]*psi.mu.w[,k]-sigma.ww[j,k] # Would sum to 0 if multiply by (N-1)/N. Similar elsewhere
junk1[,k+dim.w*(j-1)]<-ifelse(v==1,psi.mu.t[,j]*psi.mu.t[,k]-sigma.tt[j,k],0)
}
}
psi.sigma.ww<-junk
psi.sigma.tt<-junk1
if (dim.w>=2) {
psi.sigma.ww<-junk[,-which(is.na(junk[1,]))]
psi.sigma.tt<-junk1[,-which(is.na(junk[1,]))]
}

psi.mu.z<-matrix(NA,N,dim.z)
junk<-matrix(NA,nrow=N,ncol=(dim.z^2))
psi.sigma.zystar<-matrix(NA,nrow=N,ncol=dim.z)
for (j in 1:dim.z) {
psi.mu.z[,j]<-z[,j]-mu.z[j]
psi.sigma.zystar[,j]<-psi.mu.ystar*psi.mu.z[,j]-sigma.zystar[j,]
}
for (j in 1:dim.z) {
for (k in j:dim.z) {
junk[,k+dim.z*(j-1)]<-psi.mu.z[,j]*psi.mu.z[,k]-sigma.zz[j,k]
}
}
psi.sigma.zz<-junk
if (dim.z>=2) {
psi.sigma.zz<-junk[,-which(is.na(junk[1,]))]
}

psi.sigma.wz<-matrix(NA,nrow=N,ncol=dim.w*dim.z)
for (j in 1:dim.w) {
for (k in 1:dim.z) {
psi.sigma.wz[,k+dim.z*(j-1)]<-psi.mu.w[,j]*psi.mu.z[,k]-sigma.wz[j,k]
}
}
psi.sigma.ttstar<-matrix(0,nrow=N,ncol=dim.w) #### Typically, we only identify a single dimensional tstar=ystar-y.
for (j in 1:dim.w) {
psi.sigma.ttstar[,j]<-ifelse(v==1,psi.mu.t[,j]*psi.mu.tstar-sigma.ttstar[j],0)
}


mu.x<-apply(x,2,mean,na.rm=TRUE)
mu.y<-mean(y,na.rm=TRUE)

psi.mu.y<-ifelse(v==1,y-mu.y,0)

psi.sigma.ztstar<-matrix(NA,nrow=N,ncol=dim.z)
for (j in 1:dim.z) {
psi.sigma.ztstar[,j]<-(psi.mu.tstar*psi.mu.z[,j]-sigma.ztstar[j,])*v
}

psi.sigma.zt<-matrix(NA,nrow=N,ncol=dim.w*dim.z)
for (j in 1:dim.w) {
for (k in 1:dim.z) {
psi.sigma.zt[,k+dim.z*(j-1)]<-(psi.mu.t[,j]*psi.mu.z[,k]-sigma.zt[j,k])*v
}
}

psi.mu.t<-psi.mu.x<-matrix(NA,N,dim.w)
psi.sigma.xtstar<-psi.sigma.yt<-matrix(NA,nrow=N,ncol=dim.w)
for (j in 1:dim.w) {
psi.mu.t[,j]<-ifelse(v==1,t[,j]-mu.t[j],0)
psi.mu.x[,j]<-ifelse(v==1,x[,j]-mu.x[j],0)
psi.sigma.xtstar[,j]<-(psi.mu.tstar*psi.mu.x[,j]-sigma.xtstar[j,])*v
psi.sigma.yt[,j]<-(psi.mu.y*psi.mu.t[,j]-sigma.yt[j,])*v
}

psi<-cbind(psi.mu.w,
psi.sigma.ww,
psi.mu.z,
psi.sigma.zz,
psi.sigma.wz,
psi.sigma.tt,
psi.mu.ystar,
psi.sigma.wystar,
psi.sigma.zystar,
psi.sigma.ttstar,
psi.mu.tstar,
psi.mu.t,
psi.mu.x,
psi.mu.y,
psi.sigma.ztstar,
psi.sigma.zt,
psi.sigma.xtstar,
psi.sigma.yt)

B<-t(psi)%*%psi/N


dim.A<-dim.w+sum(1:dim.w)+dim.z+sum(1:dim.z)+dim.w*dim.z+sum(1:dim.w)+1+dim.w+dim.z+dim.w+
1+dim.w+dim.w+1+dim.z+dim.w*dim.z+dim.w+dim.w ##### new stuff

A<-matrix(0,dim.A,dim.A)
for (j in 1:(dim.w+sum(1:dim.w)+dim.z+sum(1:dim.z)+dim.w*dim.z) ) {
A[j,j]<-1
}
for (j in (dim.w+sum(1:dim.w)+dim.z+sum(1:dim.z)+dim.w*dim.z+1):
(dim.w+sum(1:dim.w)+dim.z+sum(1:dim.z)+dim.w*dim.z+sum(1:dim.w))) {
A[j,j]<-nv/N
}
for (j in (dim.w+sum(1:dim.w)+dim.z+sum(1:dim.z)+dim.w*dim.z+sum(1:dim.w)+1):
(dim.A-5*dim.w-2-dim.z*dim.w-dim.z)) {
A[j,j]<-1
}
for (j in (dim.A-5*dim.w-2-dim.z*dim.w-dim.z+1):dim.A) {
A[j,j]<-nv/N
}

var.theta<-solve(A)%*%B%*%solve(A)/N

var.theta.diag<-diag(var.theta)


dbeta.dmu.w<- matrix(0,nrow=dim.w+dim.z,ncol=dim.w)

dbeta.dsigma.ww<-dbeta.dsigma.tt<-matrix(NA,nrow=dim.w+dim.z,ncol=dim(psi.sigma.ww)[2])
for (j in 1:dim(psi.sigma.ww)[2]) {
junk<-rep(0,dim(psi.sigma.ww)[2])
junk[j]<-1
Dmat<-Dmat1<-matrix(0,dim.w+dim.z,dim.w+dim.z)
Dmat[1,1:dim.w]<-junk[1:dim.w]
if (dim.w>=2){
Dmat[2,2:dim.w]<-junk[(dim.w+1):(2*dim.w-1)]
if (dim.w>=3) {
Dmat[3,3:dim.w]<-junk[(2*dim.w):(3*dim.w-3)]
if (dim.w>=4) {
Dmat[4,4:dim.w]<-junk[(3*dim.w-2):(4*dim.w-6)]
if (dim.w>=5) {
Dmat[5,5:dim.w]<-junk[(4*dim.w-5):(5*dim.w-10)]
if (dim.w==6) {
Dmat[6,6:dim.w]<-junk[(5*dim.w-9):(6*dim.w-15)]
}}}}}
if (sum(diag(Dmat))==0) {
Dmat<-Dmat+t(Dmat)
}
dbeta.dsigma.ww[,j]<- -C.inv%*%Dmat%*%C.inv%*%D
dbeta.dsigma.tt[,j]<- - dbeta.dsigma.ww[,j]
}

dbeta.dmu.z<- matrix(0,nrow=dim.w+dim.z,ncol=dim.z)

dbeta.dsigma.zz<-matrix(NA,nrow=dim.w+dim.z,ncol=dim(psi.sigma.zz)[2])
for (j in 1:dim(psi.sigma.zz)[2]) {
junk<-rep(0,dim(psi.sigma.zz)[2])
junk[j]<-1
Dmat<-matrix(0,dim.w+dim.z,dim.w+dim.z)
Dmat[dim.w+1,(dim.w+1):(dim.w+dim.z)]<-junk[1:dim.z]
if (dim.z>=2){
Dmat[dim.w+2,(dim.w+2):(dim.w+dim.z)]<-junk[(dim.z+1):(2*dim.z-1)]
if (dim.z>=3) {
Dmat[dim.w+3,(dim.w+3):(dim.w+dim.z)]<-junk[(2*dim.z):(3*dim.z-3)]
if (dim.z>=4) {
Dmat[dim.w+4,(dim.w+4):(dim.w+dim.z)]<-junk[(3*dim.z-2):(4*dim.z-6)]
if (dim.z>=5) {
Dmat[dim.w+5,(dim.w+5):(dim.w+dim.z)]<-junk[(4*dim.z-5):(5*dim.z-10)]
if (dim.z==6) {
Dmat[dim.w+6,(dim.w+6):(dim.w+dim.z)]<-junk[(5*dim.z-9):(6*dim.z-15)]
}}}}}
if (sum(diag(Dmat))==0) {
Dmat<-Dmat+t(Dmat)
}
dbeta.dsigma.zz[,j]<- -C.inv%*%Dmat%*%C.inv%*%D
}

dbeta.dsigma.wz<-matrix(NA,nrow=dim.w+dim.z,ncol=dim(psi.sigma.wz)[2])
for (j in 1:dim(psi.sigma.wz)[2]) {
junk<-rep(0,dim(psi.sigma.wz)[2])
junk[j]<-1
Dmat12<-matrix(junk,nrow=dim.w,ncol=dim.z,byrow=TRUE)
Dmat21<-t(Dmat12)
Dmat<-matrix(0,dim.w+dim.z,dim.w+dim.z)
Dmat[1:dim.w,(dim.w+1):(dim.w+dim.z)]<-Dmat12
Dmat[(dim.w+1):(dim.w+dim.z),1:dim.w]<-Dmat21
dbeta.dsigma.wz[,j]<- -C.inv%*%Dmat%*%C.inv%*%D
}

dbeta.dmu.ystar<- matrix(0,nrow=dim.w+dim.z,ncol=1)

stuff<-matrix(0,nrow=dim.w+dim.z,ncol=dim.w+dim.z)
diag(stuff)<-1
dbeta.dsigma.wzystar<- C.inv%*%stuff

stuff<-matrix(0,nrow=dim.w+dim.z,ncol=dim.w)
diag(stuff)<- -1
dbeta.dsigma.ttstar<- C.inv%*%stuff


dbeta.dmu.tstar<- matrix(0,nrow=dim.w+dim.z,ncol=1)
dbeta.dmu.t<- matrix(0,nrow=dim.w+dim.z,ncol=dim.w)
dbeta.dmu.x<- matrix(0,nrow=dim.w+dim.z,ncol=dim.w)
dbeta.dmu.y<- matrix(0,nrow=dim.w+dim.z,ncol=1)

stuff<-matrix(0,nrow=dim.w+dim.z,ncol=dim.w)
diag(stuff)<- -1
dbeta.dsigma.xtstar <- dbeta.dsigma.yt <- dbeta.dsigma.ztstar <- C.inv%*%stuff

dbeta.dsigma.zt<-matrix(NA,nrow=dim.w+dim.z,ncol=dim(psi.sigma.zt)[2])
for (j in 1:dim(psi.sigma.zt)[2]) {
junk<-rep(0,dim(psi.sigma.zt)[2])
junk[j]<- -1
Dmat12<-matrix(junk,nrow=dim.w,ncol=dim.z,byrow=TRUE)
Dmat21<-t(Dmat12)
Dmat<-matrix(0,dim.w+dim.z,dim.w+dim.z)
Dmat[1:dim.w,(dim.w+1):(dim.w+dim.z)]<-Dmat12
Dmat[(dim.w+1):(dim.w+dim.z),1:dim.w]<-Dmat21
dbeta.dsigma.zt[,j]<- -C.inv%*%Dmat%*%C.inv%*%D
}


dg.theta<-cbind(dbeta.dmu.w, dbeta.dsigma.ww, dbeta.dmu.z, dbeta.dsigma.zz, dbeta.dsigma.wz, dbeta.dsigma.tt,
dbeta.dmu.ystar, dbeta.dsigma.wzystar, dbeta.dsigma.ttstar,
dbeta.dmu.tstar, dbeta.dmu.t, dbeta.dmu.x, dbeta.dmu.y,
dbeta.dsigma.ztstar, dbeta.dsigma.zt, dbeta.dsigma.xtstar, dbeta.dsigma.yt)

var.est<-dg.theta%*%var.theta%*%t(dg.theta)

return(list("CorrectEst"=est,"VarCorrectEst"=var.est))

}





date()
N<-100 ##### varies 100, 1000
num.audit<-50 ##### varies 25, 50, 100, 300, 500

sigma<-.5
sigma.x<-50
mean.x<-200
beta0<-6
beta1<- -.01
beta2<- 1
sigma.uy<-.5
sigma.ustar<-.5
mean.ustar<-0

py<- 0.2
sigma.u<-50
cor.u.ustar<- .5
p<- .2

##### Variable parameters

mean.xz<- c(0,0,-50)
mean.uyz <- c(0,1,1)
n.scenarios<-3

Nimpute<-20
Nrep<-1000

perc.bias.fixedx<-mean.sd.fixedx<-cov.prob.fixedx<-mse.fixedx<-emp.sd.fixedx<-matrix(NA,n.scenarios,length(num.audit))
perc.bias.fixed2x<-mean.sd.fixed2x<-cov.prob.fixed2x<-mse.fixed2x<-emp.sd.fixed2x<-matrix(NA,n.scenarios,length(num.audit))
perc.bias.fixed3x<-mean.sd.fixed3x<-cov.prob.fixed3x<-mse.fixed3x<-emp.sd.fixed3x<-matrix(NA,n.scenarios,length(num.audit))
perc.bias.attx<-mean.sd.attx<-cov.prob.attx<-mse.attx<-emp.sd.attx<-NULL
perc.bias.correctx<-mean.sd.correctx<-cov.prob.correctx<-mse.correctx<-emp.sd.correctx<-NULL

perc.bias.fixedz<-mean.sd.fixedz<-cov.prob.fixedz<-mse.fixedz<-emp.sd.fixedz<-matrix(NA,n.scenarios,length(num.audit))
perc.bias.fixed2z<-mean.sd.fixed2z<-cov.prob.fixed2z<-mse.fixed2z<-emp.sd.fixed2z<-matrix(NA,n.scenarios,length(num.audit))
perc.bias.fixed3z<-mean.sd.fixed3z<-cov.prob.fixed3z<-mse.fixed3z<-emp.sd.fixed3z<-matrix(NA,n.scenarios,length(num.audit))
perc.bias.attz<-mean.sd.attz<-cov.prob.attz<-mse.attz<-emp.sd.attz<-NULL
perc.bias.correctz<-mean.sd.correctz<-cov.prob.correctz<-mse.correctz<-emp.sd.correctz<-NULL

for (ii in 1:n.scenarios) {

correctx<-attenuatedx<-NULL
var.gammax<-var.correctx<-NULL
correctz<-attenuatedz<-NULL
var.gammaz<-var.correctz<-NULL
check<-matrix(NA,Nrep,length(num.audit))
var.fixedx <- fixedx<-lower.fixedx <-upper.fixedx <-matrix(NA,Nrep,length(num.audit))
var.fixed2x <- fixed2x<-lower.fixed2x <-upper.fixed2x <-matrix(NA,Nrep,length(num.audit))
var.fixed3x <- fixed3x<-lower.fixed3x <-upper.fixed3x <-matrix(NA,Nrep,length(num.audit))
var.fixedz <- fixedz<-lower.fixedz <-upper.fixedz <-matrix(NA,Nrep,length(num.audit))
var.fixed2z <- fixed2z<-lower.fixed2z <-upper.fixed2z <-matrix(NA,Nrep,length(num.audit))
var.fixed3z <- fixed3z<-lower.fixed3z <-upper.fixed3z <-matrix(NA,Nrep,length(num.audit))

lower.attenx <- upper.attenx<- lower.correctx<- upper.correctx<-NULL
lower.attenz <- upper.attenz<- lower.correctz<- upper.correctz<-NULL

for (i in 1:Nrep) {

z<-rbinom(N,1,.5)
x<-rnorm(N,mean.x+mean.xz[ii]*z,sigma.x)
s<-rbinom(N,1,p)
stuff<-rmvnorm(N, mean=c(0,mean.ustar),
sigma=matrix(c(sigma.u^2,cor.u.ustar*sigma.u*sigma.ustar,cor.u.ustar*sigma.u*sigma.ustar,sigma.ustar^2),
nrow=2))
u<-stuff[,1]
ustar<-stuff[,2]
w<-ifelse(s==1,x+u,x)
e<-rnorm(N,0,sigma)
y<-beta0+beta1*x+beta2*z+e
sy<-rbinom(N,1,py)
uy<-rnorm(N,mean.uyz[ii]*z,sigma.uy)
ystar<-y+sy*uy+s*ustar



mod.correct<-lm(y~x+z)
mod.atten<-lm(ystar~w+z)


gammax<-mod.atten$coeff[2]
attenuatedx[i]<-gammax
var.gammax[i]<-summary(mod.atten)$coeff[2,2]^2
correctx[i]<-mod.correct$coeff[2]
var.correctx[i]<-summary(mod.correct)$coeff[2,2]^2

gammaz<-mod.atten$coeff[3]
attenuatedz[i]<-gammaz
var.gammaz[i]<-summary(mod.atten)$coeff[3,2]^2
correctz[i]<-mod.correct$coeff[3]
var.correctz[i]<-summary(mod.correct)$coeff[3,2]^2

lower.attenx[i]<-attenuatedx[i]-1.96*sqrt(var.gammax[i])
upper.attenx[i]<-attenuatedx[i]+1.96*sqrt(var.gammax[i])
lower.correctx[i]<-correctx[i]-1.96*sqrt(var.correctx[i])
upper.correctx[i]<-correctx[i]+1.96*sqrt(var.correctx[i])

lower.attenz[i]<-attenuatedz[i]-1.96*sqrt(var.gammaz[i])
upper.attenz[i]<-attenuatedz[i]+1.96*sqrt(var.gammaz[i])
lower.correctz[i]<-correctz[i]-1.96*sqrt(var.correctz[i])
upper.correctz[i]<-correctz[i]+1.96*sqrt(var.correctz[i])

for (j in 1:length(num.audit)) {

nv<-num.audit[j]
v<-c(rep(1,nv),rep(0,N-nv))

ests<-audit.correct.cov(ystar,w,z,v,x,y)

fixedx[i,j]<- ests$CorrectEst[1,1]
fixedz[i,j]<- ests$CorrectEst[2,1]
var.fixedx[i,j]<- ests$VarCorrectEst[1,1]
var.fixedz[i,j]<- ests$VarCorrectEst[2,2]


ypart<-ifelse(v==1,y,ystar)
xpart<-ifelse(v==1,x,w)

mod.naive2<-lm(ypart~xpart+z)

fixed2x[i,j]<- mod.naive2$coeff[2]
fixed2z[i,j]<- mod.naive2$coeff[3]
var.fixed2x[i,j]<-summary(mod.naive2)$coeff[2,2]^2
var.fixed2z[i,j]<-summary(mod.naive2)$coeff[3,2]^2

check[i,j]<-sum(s*v)

if (check[i,j]>0) {

#### multiple imputation stuff

modx<-lm(x~w+ystar+z,subset=v==1)
sigmas<-sqrt(summary(modx)$sigma^2*rchisq(Nimpute, summary(modx)$df[2], ncp=0)/summary(modx)$df[2])
desmatrix <- rmvnorm(Nimpute, mean=modx$coeff, sigma=summary(modx)$sigma^2*summary(modx)$cov)

mody<-lm(y~w+ystar+x+z,subset=v==1)
sigmasy<-sqrt(summary(mody)$sigma^2*rchisq(Nimpute, summary(mody)$df[2], ncp=0)/summary(mody)$df[2])
desmatrixy <- rmvnorm(Nimpute, mean=mody$coeff, sigma=summary(mody)$sigma^2*summary(mody)$cov)

x.coeff.mi1<-x.coeff.var.mi1<-NULL
x.coeff.mi2<-x.coeff.var.mi2<-NULL
for (k in 1:Nimpute) {
pred.x<-as.vector(desmatrix[k,]%*%t(cbind(1,w,ystar,z))+rnorm(N,0,sigmas[k]))
x.impute<-ifelse(v==0,pred.x,x)
pred.y<-as.vector(desmatrixy[k,]%*%t(cbind(1,w,ystar,pred.x,z))+rnorm(N,0,sigmasy[k]))
y.impute<-ifelse(v==0,pred.y,y)
mod.imputed<-lm(y.impute~x.impute+z)
x.coeff.mi1[k]<-mod.imputed$coeff[2]
x.coeff.var.mi1[k]<-summary(mod.imputed)$coeff[2,2]^2
x.coeff.mi2[k]<-mod.imputed$coeff[3]
x.coeff.var.mi2[k]<-summary(mod.imputed)$coeff[3,2]^2
}
fixed3x[i,j]<-mean(x.coeff.mi1)
var.fixed3x[i,j]<-mean(x.coeff.var.mi1)+(Nimpute+1)*var(x.coeff.mi1)/Nimpute
fixed3z[i,j]<-mean(x.coeff.mi2)
var.fixed3z[i,j]<-mean(x.coeff.var.mi2)+(Nimpute+1)*var(x.coeff.mi2)/Nimpute

} else {
fixed3x[i,j]<-gammax
fixed3z[i,j]<-gammaz
var.fixed3x[i,j]<-var.gammax[i]
var.fixed3z[i,j]<-var.gammaz[i]
}

lower.fixedx[i,j]<-fixedx[i,j]-1.96*sqrt(var.fixedx[i,j])
upper.fixedx[i,j]<-fixedx[i,j]+1.96*sqrt(var.fixedx[i,j])
lower.fixedz[i,j]<-fixedz[i,j]-1.96*sqrt(var.fixedz[i,j])
upper.fixedz[i,j]<-fixedz[i,j]+1.96*sqrt(var.fixedz[i,j])

lower.fixed2x[i,j]<-fixed2x[i,j]-1.96*sqrt(var.fixed2x[i,j])
upper.fixed2x[i,j]<-fixed2x[i,j]+1.96*sqrt(var.fixed2x[i,j])
lower.fixed2z[i,j]<-fixed2z[i,j]-1.96*sqrt(var.fixed2z[i,j])
upper.fixed2z[i,j]<-fixed2z[i,j]+1.96*sqrt(var.fixed2z[i,j])

lower.fixed3x[i,j]<-fixed3x[i,j]-1.96*sqrt(var.fixed3x[i,j])
upper.fixed3x[i,j]<-fixed3x[i,j]+1.96*sqrt(var.fixed3x[i,j])
lower.fixed3z[i,j]<-fixed3z[i,j]-1.96*sqrt(var.fixed3z[i,j])
upper.fixed3z[i,j]<-fixed3z[i,j]+1.96*sqrt(var.fixed3z[i,j])


}


}


perc.bias.attx[ii]<-(beta1-mean(attenuatedx))*100
mean.sd.attx[ii]<-mean((upper.attenx-lower.attenx)/(2*1.96))*100
cov.prob.attx[ii]<-sum(lower.attenx<beta1&upper.attenx>beta1)/Nrep
mse.attx[ii]<-mean((attenuatedx-beta1)^2)

perc.bias.fixedx[ii,]<-(beta1-apply(fixedx,2,mean))*100
mean.sd.fixedx[ii,]<-apply(sqrt(var.fixedx),2,mean)*100
cov.prob.fixedx[ii,]<-apply(lower.fixedx<beta1&upper.fixedx>beta1,2,mean)
mse.fixedx[ii,]<-apply((fixedx-beta1)^2,2,mean)
emp.sd.fixedx[ii,]<-sqrt(apply(fixedx,2,var))

perc.bias.fixed2x[ii,]<-(beta1-apply(fixed2x,2,mean))*100
mean.sd.fixed2x[ii,]<-apply(sqrt(var.fixed2x),2,mean)*100
cov.prob.fixed2x[ii,]<-apply(lower.fixed2x<beta1&upper.fixed2x>beta1,2,mean)
mse.fixed2x[ii,]<-apply((fixed2x-beta1)^2,2,mean)
emp.sd.fixed2x[ii,]<-sqrt(apply(fixed2x,2,var))

perc.bias.fixed3x[ii,]<-(beta1-apply(fixed3x,2,mean))*100
mean.sd.fixed3x[ii,]<-apply(sqrt(var.fixed3x),2,mean)*100
cov.prob.fixed3x[ii,]<-apply(lower.fixed3x<beta1&upper.fixed3x>beta1,2,mean)
mse.fixed3x[ii,]<-apply((fixed3x-beta1)^2,2,mean)
emp.sd.fixed3x[ii,]<-sqrt(apply(fixed3x,2,var))

perc.bias.correctx[ii]<-(beta1-mean(correctx))*100
mean.sd.correctx[ii]<-mean((upper.correctx-lower.correctx)/(2*1.96))*100
cov.prob.correctx[ii]<-sum(lower.correctx<beta1&upper.correctx>beta1)/Nrep
mse.correctx[ii]<-mean((correctx-beta1)^2)


perc.bias.attz[ii]<-(beta2-mean(attenuatedz))*100
mean.sd.attz[ii]<-mean((upper.attenz-lower.attenz)/(2*1.96))*100
cov.prob.attz[ii]<-sum(lower.attenz<beta2&upper.attenz>beta2)/Nrep
mse.attz[ii]<-mean((attenuatedz-beta2)^2)

perc.bias.fixedz[ii,]<-(beta2-apply(fixedz,2,mean))*100
mean.sd.fixedz[ii,]<-apply(sqrt(var.fixedz),2,mean)*100
cov.prob.fixedz[ii,]<-apply(lower.fixedz<beta2&upper.fixedz>beta2,2,mean)
mse.fixedz[ii,]<-apply((fixedz-beta2)^2,2,mean)
emp.sd.fixedz[ii,]<-sqrt(apply(fixedz,2,var))

perc.bias.fixed2z[ii,]<-(beta2-apply(fixed2z,2,mean))*100
mean.sd.fixed2z[ii,]<-apply(sqrt(var.fixed2z),2,mean)*100
cov.prob.fixed2z[ii,]<-apply(lower.fixed2z<beta2&upper.fixed2z>beta2,2,mean)
mse.fixed2z[ii,]<-apply((fixed2z-beta2)^2,2,mean)
emp.sd.fixed2z[ii,]<-sqrt(apply(fixed2z,2,var))

perc.bias.fixed3z[ii,]<-(beta2-apply(fixed3z,2,mean))*100
mean.sd.fixed3z[ii,]<-apply(sqrt(var.fixed3z),2,mean)*100
cov.prob.fixed3z[ii,]<-apply(lower.fixed3z<beta2&upper.fixed3z>beta2,2,mean)
mse.fixed3z[ii,]<-apply((fixed3z-beta2)^2,2,mean)
emp.sd.fixed3z[ii,]<-sqrt(apply(fixed3z,2,var))

perc.bias.correctz[ii]<-(beta2-mean(correctz))*100
mean.sd.correctz[ii]<-mean((upper.correctz-lower.correctz)/(2*1.96))*100
cov.prob.correctz[ii]<-sum(lower.correctz<beta2&upper.correctz>beta2)/Nrep
mse.correctz[ii]<-mean((correctz-beta2)^2)

}