## R code used for the simulations in the manuscript "Accounting for data errors discovered from an audit in multiple linear regression" by Bryan E. Shepherd and Chang Yu.

## audit.correct() takes information from an audit, and corrects naive estimates when there are errors in the
## predictor variable. This function employs the methods given in Section 2.1, and returns the estimates
## CorrectEst1, CorrectEst2, CorrectEst3, and their variances, which are based on (lambda1, lambda2, and lambda3),
## as outlined in Section 2.1.


audit.correct<-function(y,w,v,x,s) {

mod.atten<-lm(y~w)
gamma1<-mod.atten$coeff[2]
var.gamma1<-summary(mod.atten)$coeff[2,2]^2
mu.w<-mean(w)
var.w<-var(w)
p.est<-mean(s[v==1])
u.est<-w[s*v==1]-x[s*v==1]
var.u<-ifelse(sum(s*v)==0,0,
ifelse(sum(s*v)==1,u.est^2,var(u.est)))
var.x<-var(x[v==1])
mu.x<-mean(x[v==1])

atten<-var.x/(var.x+var.u*p.est)
fixed<-gamma1/atten

atten1<-(var.w-var.u*p.est)/var.w
fixed1<-gamma1/atten1

t.est<-w[v==1]-x[v==1]
var.t<-ifelse(sum(s*v)==0,0,
ifelse(sum(s*v)==1,t.est^2,var(t.est)))
atten2<-var.x/(var.x+var.t)
fixed2<-mod.atten$coeff[2]/atten2

check<-sum(s*v)

if (check>0) {
A<-matrix(0,6,6)
A[1,1]<-1
A[1,2]<-mean(w)
A[2,1]<-mean(w)
A[2,2]<-sum(w^2)/N
A[3,3]<-sum(v)/N
A[4,3]<-2*sum(v*(x-mu.x))/N
A[4,4]<-sum(v)/N
A[5,5]<-sum(v)/N
A[6,6]<-sum(s*v)/N

psi<-cbind((y-mod.atten$coeff[1]-gamma1*w),
(y-mod.atten$coeff[1]-gamma1*w)*w,
(x-mu.x)*v,
((x-mu.x)^2-var.x)*v,
(s-p.est)*v,
((w-x)^2-var.u)*s*v)
B<-t(psi)%*%psi/N
var.theta<-solve(A)%*%B%*%solve(A)/N

dg<-c(0, 1/atten, 0, -gamma1*var.u*p.est/var.x^2, gamma1*var.u/var.x, gamma1*p.est/var.x )

var.fixed<-dg%*%var.theta%*%dg


A1<-matrix(0,6,6)
A1[1,1]<-1
A1[1,2]<-mean(w)
A1[2,1]<-mean(w)
A1[2,2]<-sum(w^2)/N
A1[3,3]<-1
A1[4,3]<-2*sum(w-mu.w)/N
A1[4,4]<-1
A1[5,5]<-sum(v)/N
A1[6,6]<-sum(s*v)/N

psi1<-cbind((y-mod.atten$coeff[1]-gamma1*w),
(y-mod.atten$coeff[1]-gamma1*w)*w,
(w-mu.w),
((w-mu.w)^2-var.w),
(s-p.est)*v,
((w-x)^2-var.u)*s*v)
B1<-t(psi1)%*%psi1/N
var.theta1<-solve(A1)%*%B1%*%solve(A1)/N

dg1<-c(0, var.w/(var.w-var.u*p.est), 0, -gamma1*var.u*p.est/(var.w-var.u*p.est)^2,
gamma1*var.w*var.u/(var.w-var.u*p.est)^2,gamma1*var.w*p.est/(var.w-var.u*p.est)^2)

var.fixed1<-dg1%*%var.theta1%*%dg1


A2<-matrix(0,5,5)
A2[1,1]<-1
A2[1,2]<-mean(w)
A2[2,1]<-mean(w)
A2[2,2]<-sum(w^2)/N
A2[3,3]<-sum(v)/N
A2[4,3]<-2*sum(v*(x-mu.x))/N
A2[4,4]<-sum(v)/N
A2[5,5]<-sum(v)/N

psi2<-cbind((y-mod.atten$coeff[1]-gamma1*w),
(y-mod.atten$coeff[1]-gamma1*w)*w,
(x-mu.x)*v,
((x-mu.x)^2-var.x)*v,
((w-x)^2-var.t)*v)
B2<-t(psi2)%*%psi2/N
var.theta2<-solve(A2)%*%B2%*%solve(A2)/N

dg2<-c(0, 1/atten2, 0, -gamma1*var.t/var.x^2, gamma1/var.x)

var.fixed2<-dg2%*%var.theta2%*%dg2

} else {
var.fixed<-var.fixed1<-var.fixed2<-var.gamma1
}

return(list("CorrectEst1"=fixed,"VarCorrectEst1"=var.fixed,
"CorrectEst2"=fixed1,"VarCorrectEst2"=var.fixed1,
"CorrectEst3"=fixed2,"VarCorrectEst3"=var.fixed2))

}





##### Simulation code when errors sometimes in predictor variable

N<-1000
p<-.05
pv<-.05
sigma<-.5
sigma.x<-50
sigma.u<-50
mean.x<-200
beta0<-6
beta1<- -.01

error.probs<-c(.05,.1,.2,.3)
num.audit<-c(25,50,100,200,300)

Nrep<-5000

perc.bias.audit1<-mean.sd.audit1<-cov.prob.audit1<-mse.audit1<-emp.sd.audit1<-matrix(NA,length(error.probs),length(num.audit))
perc.bias.audit2<-mean.sd.audit2<-cov.prob.audit2<-mse.audit2<-emp.sd.audit2<-matrix(NA,length(error.probs),length(num.audit))
perc.bias.audit3<-mean.sd.audit3<-cov.prob.audit3<-mse.audit3<-emp.sd.audit3<-matrix(NA,length(error.probs),length(num.audit))

perc.bias.att<-mean.sd.att<-cov.prob.att<-mse.att<-emp.sd.att<-NULL
perc.bias.complete<-mean.sd.complete<-cov.prob.complete<-mse.complete<-emp.sd.complete<-NULL


for (ii in 1:length(error.probs)) {

complete<-attenuated<-NULL
var.attenuated<-var.complete<-NULL
var.audit1 <- audit1<-lower.audit1 <-upper.audit1 <-matrix(NA,Nrep,length(num.audit))
var.audit2<-audit2<-lower.audit2<-upper.audit2<-matrix(NA,Nrep,length(num.audit))
var.audit3<-audit3<-lower.audit3<-upper.audit3<-matrix(NA,Nrep,length(num.audit))
lower.atten <- upper.atten<- lower.complete<- upper.complete<-NULL

p<-error.probs[ii]

for (i in 1:Nrep) {

x<-rnorm(N,mean.x,sigma.x)
s<-rbinom(N,1,p)
u<-rnorm(N,0,sigma.u)
w<-ifelse(s==1,x+u,x)
e<-rnorm(N,0,sigma)
y<-beta0+beta1*x+e

mod.complete<-lm(y~x)
mod.atten<-lm(y~w)

attenuated[i]<-mod.atten$coeff[2]
complete[i]<-mod.complete$coeff[2]
var.attenuated[i]<-summary(mod.atten)$coeff[2,2]^2
var.complete[i]<-summary(mod.complete)$coeff[2,2]^2

lower.atten[i]<-attenuated[i]-1.96*sqrt(var.attenuated[i])
upper.atten[i]<-attenuated[i]+1.96*sqrt(var.attenuated[i])
lower.complete[i]<-complete[i]-1.96*sqrt(var.complete[i])
upper.complete[i]<-complete[i]+1.96*sqrt(var.complete[i])

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

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

ans<-audit.correct(y,w,v,x,s)

audit1[i,j] <-ans$CorrectEst1
audit2[i,j]<-ans$CorrectEst2
audit3[i,j]<-ans$CorrectEst3
var.audit1[i,j] <-ans$VarCorrectEst1
var.audit2[i,j]<-ans$VarCorrectEst2
var.audit3[i,j]<-ans$VarCorrectEst3

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

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

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


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

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

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

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

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

}

round(perc.bias.att*100,2)
round(mean.sd.att,3)
round(cov.prob.att,3)
round(mse.att*10^7,1)

round(perc.bias.complete*100,2)
round(mean.sd.complete,3)
round(cov.prob.complete,3)
round(mse.complete*10^7,1)

round(perc.bias.audit1*100,2)
round(mean.sd.audit1,3)
round(cov.prob.audit1,3)
round(mse.audit1*10^7,1)

round(perc.bias.audit2*100,2)
round(mean.sd.audit2,3)
round(cov.prob.audit2,3)
round(mse.audit2*10^7,1)

round(perc.bias.audit3*100,2)
round(mean.sd.audit3,3)
round(cov.prob.audit3,3)
round(mse.audit3*10^7,1)








## audit.correct2() takes information from an audit, and corrects naive estimates when there are errors in the
## predictor and outcome variables which are potentially correlated. This function employs the methods given in
## Section 2.3, and returns the estimate CorrectEst and its variance as outlined in Section 2.3.


audit.correct2<-function(y,w,v,x,s) {

mod.atten<-lm(ystar~w)
gamma1<-mod.atten$coeff[2]
var.gamma1<-summary(mod.atten)$coeff[2,2]^2

var.w<-var(w)
mu.w<-mean(w)

p.est<-mean(s[v==1])
u.est<-w[s*v==1]-x[s*v==1]
var.u<-ifelse(sum(s*v)==0,0,
ifelse(sum(s*v)==1,u.est^2,var(u.est)*((sum(s*v)-1)/sum(s*v))))
var.x<-var(x[v==1])*((sum(v)-1)/sum(v))
mu.x<-mean(x[v==1])

resid.y<-y[v*s==1]-ystar[v*s==1]
resid.x<-x[v*s==1]-w[v*s==1]
cov.est<-ifelse(sum(s*v)==0,0,
ifelse(sum(s*v)==1,resid.y*resid.x,sum(resid.y*resid.x)/(sum(s*v)-1)))
mu.star<-mean(resid.y)

fixed<-(gamma1*(var.x+var.u*p.est) - p.est*cov.est)/var.x

check<-sum(s*v)

if (check>0) {
A<-matrix(0,7,7)
A[1,1]<-1
A[1,2]<-mean(w)
A[2,1]<-mean(w)
A[2,2]<-sum(w^2)/N
A[3,3]<-sum(v)/N
A[4,3]<-2*sum(v*(x-mu.x))/N
A[4,4]<-sum(v)/N
A[5,5]<-sum(v)/N
A[6,6]<-sum(s*v)/N
A[7,7]<-sum(s*v)/N

psi<-cbind((ystar-mod.atten$coeff[1]-gamma1*w),
(ystar-mod.atten$coeff[1]-gamma1*w)*w,
(x-mu.x)*v,
((x-mu.x)^2-var.x)*v,
(s-p.est)*v,
((w-x)^2-var.u)*s*v,
((ystar-y)*(w-x)-cov.est)*s*v
)
B<-t(psi)%*%psi/N
var.theta<-solve(A)%*%B%*%solve(A)/N

dg<-c(0, (var.x+p*var.u)/var.x, 0, -(gamma1*var.u*p.est-p.est*cov.est)/var.x^2,
(gamma1*var.u-cov.est)/var.x, gamma1*p.est/var.x, -p.est/var.x )

var.fixed<-dg%*%var.theta%*%dg

} else {
var.fixed<-var.gamma1
}

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

}









##### Simulation code when errors sometimes in both predictor and outcome variables


library(mvtnorm)

N<-1000
py<-.2
sigma<-.5
sigma.x<-50
sigma.u<-50
mean.x<-200
beta0<-6
beta1<- -.01
sigma.uy<-1
cor.u.ustar<- -0.5 ## Also used 0 and 0.5
sigma.ustar<-.5
mean.ustar<-0

error.probs<-c(.05,.1,.2,.3)
num.audit<-c(25,50,100,200,300)

Nrep<-5000

perc.bias.fixed<-mean.sd.fixed<-cov.prob.fixed<-mse.fixed<-emp.sd.fixed<-matrix(NA,length(error.probs),length(num.audit))
perc.bias.att<-mean.sd.att<-cov.prob.att<-mse.att<-emp.sd.att<-NULL
perc.bias.complete<-mean.sd.complete<-cov.prob.complete<-mse.complete<-emp.sd.complete<-NULL

for (ii in 1:length(error.probs)) {

complete<-attenuated<-NULL
var.attenuated<-var.complete<-NULL
check<-matrix(NA,Nrep,length(num.audit))
var.fixed <- fixed<-lower.fixed <-upper.fixed <-matrix(NA,Nrep,length(num.audit))
lower.atten <- upper.atten<- lower.complete<- upper.complete<-NULL

p<-error.probs[ii]

for (i in 1:Nrep) {

x<-rnorm(N,mean.x,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+e
sy<-rbinom(N,1,py)
uy<-rnorm(N,0,sigma.uy)
ystar<-y+sy*uy+s*ustar

mod.complete<-lm(y~x)
mod.atten<-lm(ystar~w)

attenuated[i]<-mod.atten$coeff[2]
complete[i]<-mod.complete$coeff[2]
var.attenuated[i]<-summary(mod.atten)$coeff[2,2]^2
var.complete[i]<-summary(mod.complete)$coeff[2,2]^2

lower.atten[i]<-attenuated[i]-1.96*sqrt(var.attenuated[i])
upper.atten[i]<-attenuated[i]+1.96*sqrt(var.attenuated[i])
lower.complete[i]<-complete[i]-1.96*sqrt(var.complete[i])
upper.complete[i]<-complete[i]+1.96*sqrt(var.complete[i])

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

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

ans<-audit.correct2(y,w,v,x,s)

fixed[i,j]<-ans$CorrectEst
var.fixed[i,j]<-ans$VarCorrectEst

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

}

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

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

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

}


round(perc.bias.att*100,2)
round(mean.sd.att,3)
round(cov.prob.att,3)
round(mse.att*10^7,1)

round(perc.bias.complete*100,2)
round(mean.sd.complete,3)
round(cov.prob.complete,3)
round(mse.complete*10^7,1)

round(perc.bias.fixed*100,2)
round(mean.sd.fixed,3)
round(cov.prob.fixed,3)
round(mse.fixed*10^7,1)

#############################################################################################
########### Simulations for covariates


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)
}
tstar<-ifelse(v==1,ystar-y,NA) ### In practice we can't differentiate between tstar1 and tstar2, but that doesn't matter.

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")


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,t[,j]*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,t[,j]*tstar-sigma.ttstar[j],0)
# for (k in 1:dim.w) {
# psi.sigma.ttstar[,k+dim.w*(j-1)]<-v*(s[,j]*u[,j]*s[,k]*ustar[,k]-sigma.ttstar[j,k])
# }
}


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)

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

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-dim.w)) {
A[j,j]<-1
}
for (j in (dim.A-dim.w+1):dim.A) {
A[j,j]<-nv/N
}

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

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

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
C[(dim.w+1):(dim.w+dim.z),1:dim.w]<-t(sigma.wz)
C[(dim.w+1):(dim.w+dim.z),(dim.w+1):(dim.w+dim.z)]<-sigma.zz

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

C.inv<-solve(C)

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

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-13)]
}}}}}
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-13)]
}}}}}
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

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)

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

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

}







sim.correct.cov<-function(N,py,sigma,sigma.x,sigma.u,mean.x,sigma.z,mean.z,cor.xz,beta0,betax,betaz,sigma.uy=1,cor.u.ustar,sigma.ustar,mean.ustar,
error.probs,num.audit,Nrep,filename) {
library(mvtnorm)
perc.bias.fixed.x=mean.sd.fixed.x=cov.prob.fixed.x=mse.fixed.x=emp.sd.fixed.x=matrix(NA,length(error.probs),length(num.audit))
perc.bias.att.x=mean.sd.att.x=cov.prob.att.x=mse.att.x=emp.sd.att.x=NULL
perc.bias.complete.x=mean.sd.complete.x=cov.prob.complete.x=mse.complete.x=emp.sd.complete.x=NULL
perc.bias.fixed.z=mean.sd.fixed.z=cov.prob.fixed.z=mse.fixed.z=emp.sd.fixed.z=matrix(NA,length(error.probs),length(num.audit))
perc.bias.att.z=mean.sd.att.z=cov.prob.att.z=mse.att.z=emp.sd.att.z=NULL
perc.bias.complete.z=mean.sd.complete.z=cov.prob.complete.z=mse.complete.z=emp.sd.complete.z=NULL

for (ii in 1:length(error.probs)) {
complete.x=attenuated.x=NULL
var.attenuated.x=var.complete.x=NULL
complete.z=attenuated.z=NULL
var.attenuated.z=var.complete.z=NULL
check=matrix(NA,Nrep,length(num.audit))
var.fixed.x = fixed.x=lower.fixed.x upper.fixed.x =matrix(NA,Nrep,length(num.audit))
lower.atten.x = upper.atten.x lower.complete.x= upper.complete.x=NULL
var.fixed.z = fixed.z=lower.fixed.z upper.fixed.z =matrix(NA,Nrep,length(num.audit))
lower.atten.z = upper.atten.z lower.complete.z= upper.complete.z=NULL
p=error.probs[ii]
for (i in 1:Nrep) {

Xmat<-rmvnorm(N, mean=c(mean.x,mean.z),sigma=matrix(c(sigma.x^2,rep(cor.xz*sigma.x*sigma.z,2),sigma.z^2),nrow=2))
x=Xmat[,1]
z=Xmat[,2]
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+betax*x+betaz*z+e
sy=rbinom(N,1,py)
uy=rnorm(N,0,sigma.uy)
ystar=y+sy*uy+s*ustar
mod.complete=lm(y~x+z)
mod.atten=lm(ystar~w+z)

attenuated.x[i]=mod.atten$coeff[2]
complete.x[i]=mod.complete$coeff[2]
var.attenuated.x[i]=summary(mod.atten)$coeff[2,2]^2
var.complete.x[i]=summary(mod.complete)$coeff[2,2]^2

attenuated.z[i]=mod.atten$coeff[3]
complete.z[i]=mod.complete$coeff[3]
var.attenuated.z[i]=summary(mod.atten)$coeff[3,2]^2
var.complete.z[i]=summary(mod.complete)$coeff[3,2]^2

lower.atten.x[i]=attenuated.x[i]-1.96*sqrt(var.attenuated.x[i])
upper.atten.x[i]=attenuated.x[i]+1.96*sqrt(var.attenuated.x[i])
lower.complete.x[i]=complete.x[i]-1.96*sqrt(var.complete.x[i])
upper.complete.x[i]=complete.x[i]+1.96*sqrt(var.complete.x[i])

lower.atten.z[i]=attenuated.z[i]-1.96*sqrt(var.attenuated.z[i])
upper.atten.z[i]=attenuated.z[i]+1.96*sqrt(var.attenuated.z[i])
lower.complete.z[i]=complete.z[i]-1.96*sqrt(var.complete.z[i])
upper.complete.z[i]=complete.z[i]+1.96*sqrt(var.complete.z[i])

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

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

ans=audit.correct.cov(ystar,w,z,v,x,y)

fixed.x[i,j]=ans$CorrectEst[1]
var.fixed.x[i,j]=ans$VarCorrectEst[1,1]

fixed.z[i,j]=ans$CorrectEst[2]
var.fixed.z[i,j]=ans$VarCorrectEst[2,2]

lower.fixed.x[i,j]=fixed.x[i,j]-1.96*sqrt(var.fixed.x[i,j])
upper.fixed.x[i,j]=fixed.x[i,j]+1.96*sqrt(var.fixed.x[i,j])

lower.fixed.z[i,j]=fixed.z[i,j]-1.96*sqrt(var.fixed.z[i,j])
upper.fixed.z[i,j]=fixed.z[i,j]+1.96*sqrt(var.fixed.z[i,j])
}

}

perc.bias.att.x[ii]=(mean(attenuated.x)-betax)/betax
mean.sd.att.x[ii]=mean((upper.atten.x-lower.atten.x)/(2*1.96))*100
cov.prob.att.x[ii]=sum(lower.atten.x<betax&upper.atten.x>betax)/Nrep
mse.att.x[ii]=mean((attenuated.x-betax)^2)

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

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

perc.bias.att.z[ii]=(mean(attenuated.z)-betaz)/betaz
mean.sd.att.z[ii]=mean((upper.atten.z-lower.atten.z)/(2*1.96))
cov.prob.att.z[ii]=sum(lower.atten.z<betaz&upper.atten.z>betaz)/Nrep
mse.att.z[ii]=mean((attenuated.z-betaz)^2)

perc.bias.fixed.z[ii,]=(apply(fixed.z,2,mean)-betaz)/betaz
mean.sd.fixed.z[ii,]=apply(sqrt(var.fixed.z),2,mean)
cov.prob.fixed.z[ii,]=apply(lower.fixed.z<betaz & upper.fixed.z>betaz,2,mean)
mse.fixed.z[ii,]=apply((fixed.z-betaz)^2,2,mean)
emp.sd.fixed.z[ii,]=sqrt(apply(fixed.z,2,var))

perc.bias.complete.z[ii]=(mean(complete.z)-betaz)/betaz
mean.sd.complete.z[ii]=mean((upper.complete.z-lower.complete.z)/(2*1.96))
cov.prob.complete.z[ii]=sum(lower.complete.z<betaz & upper.complete.z>betaz)/Nrep
mse.complete.z[ii]=mean((complete.z-betaz)^2)

}

save(perc.bias.att.x,mean.sd.att.x,cov.prob.att.x,mse.att.x,
perc.bias.complete.x,mean.sd.complete.x,cov.prob.complete.x,mse.complete.x,
perc.bias.fixed.x,mean.sd.fixed.x,cov.prob.fixed.x,mse.fixed.x,
perc.bias.att.z,mean.sd.att.z,cov.prob.att.z,mse.att.z,
perc.bias.complete.z,mean.sd.complete.z,cov.prob.complete.z,mse.complete.z,
perc.bias.fixed.z,mean.sd.fixed.z,cov.prob.fixed.z,mse.fixed.z,
file=filename)
}


sim.correct.cov(N=1000,py=.2,sigma=.5,sigma.x=50,sigma.u=50,mean.x=200,sigma.z=1,mean.z=0,cor.xz=.3,
beta0=6,betax=-0.01,betaz=1,sigma.uy=1,cor.u.ustar=0.5,
sigma.ustar=.5,mean.ustar=0,error.probs=c(.05,.1,.2,.3),num.audit=c(25,50,100,200,300),Nrep=5000,
filename="sims-pred-errors-orig-PosCor-cov.rda")

sim.correct.cov(N=1000,py=.2,sigma=.5,sigma.x=50,sigma.u=20,mean.x=200,sigma.z=1,mean.z=0,cor.xz=.3,
beta0=6,betax=-0.01,betaz=1,sigma.uy=1,cor.u.ustar=0.5,
sigma.ustar=.5,mean.ustar=0,error.probs=c(.05,.1,.2,.3),num.audit=c(25,50,100,200,300),Nrep=5000,
filename="sims-pred-errors-su20-1000PosCor-cov.rda")

sim.correct.cov(N=1000,py=.2,sigma=.5,sigma.x=50,sigma.u=50,mean.x=200,sigma.z=1,mean.z=0,cor.xz=.3,
beta0=6,betax=-0.01,betaz=1,sigma.uy=1,cor.u.ustar=0,
sigma.ustar=.5,mean.ustar=0,error.probs=c(.05,.1,.2,.3),num.audit=c(25,50,100,200,300),Nrep=5000,
filename="sims-pred-errors-su50-1000ZerCor-cov.rda")

sim.correct.cov(N=1000,py=.2,sigma=.5,sigma.x=50,sigma.u=20,mean.x=200,sigma.z=1,mean.z=0,cor.xz=.3,
beta0=6,betax=-0.01,betaz=1,sigma.uy=1,cor.u.ustar=0,
sigma.ustar=.5,mean.ustar=0,error.probs=c(.05,.1,.2,.3),num.audit=c(25,50,100,200,300),Nrep=5000,
filename="sims-pred-errors-su20-1000ZerCor-cov.rda")


sim.correct.cov(N=1000,py=.2,sigma=.5,sigma.x=50,sigma.u=50,mean.x=200,sigma.z=1,mean.z=0,cor.xz=.3,
beta0=6,betax=-0.01,betaz=1,sigma.uy=1,cor.u.ustar=-0.5,
sigma.ustar=.5,mean.ustar=0,error.probs=c(.05,.1,.2,.3),num.audit=c(25,50,100,200,300),Nrep=5000,
filename="sims-pred-errors-su50-1000NegCor-cov.rda")

sim.correct.cov(N=1000,py=.2,sigma=.5,sigma.x=50,sigma.u=20,mean.x=200,sigma.z=1,mean.z=0,cor.xz=.3,
beta0=6,betax=-0.01,betaz=1,sigma.uy=1,cor.u.ustar=-0.5,
sigma.ustar=.5,mean.ustar=0,error.probs=c(.05,.1,.2,.3),num.audit=c(25,50,100,200,300),Nrep=5000,
filename="sims-pred-errors-su20-1000NegCor-cov.rda")


sim.correct.cov(N=500,py=.2,sigma=.5,sigma.x=50,sigma.u=50,mean.x=200,sigma.z=1,mean.z=0,cor.xz=.3,
beta0=6,betax=-0.01,betaz=1,sigma.uy=1,cor.u.ustar=0.5,
sigma.ustar=.5,mean.ustar=0,error.probs=c(.05,.1,.2,.3),num.audit=c(25,50,100,200,300),Nrep=5000,
filename="sims-pred-errors-su50-500PosCor-cov.rda")

sim.correct.cov(N=500,py=.2,sigma=.5,sigma.x=50,sigma.u=20,mean.x=200,sigma.z=1,mean.z=0,cor.xz=.3,
beta0=6,betax=-0.01,betaz=1,sigma.uy=1,cor.u.ustar=0.5,
sigma.ustar=.5,mean.ustar=0,error.probs=c(.05,.1,.2,.3),num.audit=c(25,50,100,200,300),Nrep=5000,
filename="sims-pred-errors-su20-500PosCor-cov.rda")


sim.correct.cov(N=500,py=.2,sigma=.5,sigma.x=50,sigma.u=50,mean.x=200,sigma.z=1,mean.z=0,cor.xz=.3,
beta0=6,betax=-0.01,betaz=1,sigma.uy=1,cor.u.ustar=0,
sigma.ustar=.5,mean.ustar=0,error.probs=c(.05,.1,.2,.3),num.audit=c(25,50,100,200,300),Nrep=5000,
filename="sims-pred-errors-su50-500ZerCor-cov.rda")

sim.correct.cov(N=500,py=.2,sigma=.5,sigma.x=50,sigma.u=20,mean.x=200,sigma.z=1,mean.z=0,cor.xz=.3,
beta0=6,betax=-0.01,betaz=1,sigma.uy=1,cor.u.ustar=0,
sigma.ustar=.5,mean.ustar=0,error.probs=c(.05,.1,.2,.3),num.audit=c(25,50,100,200,300),Nrep=5000,
filename="sims-pred-errors-su20-500ZerCor-cov.rda")


sim.correct.cov(N=500,py=.2,sigma=.5,sigma.x=50,sigma.u=50,mean.x=200,sigma.z=1,mean.z=0,cor.xz=.3,
beta0=6,betax=-0.01,betaz=1,sigma.uy=1,cor.u.ustar=-0.5,
sigma.ustar=.5,mean.ustar=0,error.probs=c(.05,.1,.2,.3),num.audit=c(25,50,100,200,300),Nrep=5000,
filename="sims-pred-errors-su50-500NegCor-cov.rda")

sim.correct.cov(N=500,py=.2,sigma=.5,sigma.x=50,sigma.u=20,mean.x=200,sigma.z=1,mean.z=0,cor.xz=.3,
beta0=6,betax=-0.01,betaz=1,sigma.uy=1,cor.u.ustar=-0.5,
sigma.ustar=.5,mean.ustar=0,error.probs=c(.05,.1,.2,.3),num.audit=c(25,50,100,200,300),Nrep=5000,
filename="sims-pred-errors-su20-500NegCor-cov.rda")

###################################################################################
###### Analyses in the data example
###################################################################################

rm(list=ls())
load("analysis-data-CCASAnet-audit.rda")
source("audit-correct-cov.R")

### The source file: audit-correct-cov.R

z1<-ifelse(site=="argentina",1,0)
z2<-ifelse(site=="brazil",1,0)
z3<-ifelse(site=="chile",1,0)
z4<-ifelse(site=="honduras",1,0)
z5<-ifelse(site=="mexico",1,0)
z<-cbind(z1,z2,z3,z4,z5)

cat<-ifelse(site=="argentina","a",
ifelse(site=="brazil","b",
ifelse(site=="chile","c",
ifelse(site=="honduras","h",
ifelse(site=="mexico","m",
ifelse(site=="peru","p",NA))))))
poss.cat<-c("a","b","c","h","m","p")
p.cat<-table(cat)/length(cat)

mod.naive.full<-lm(ystar~w+w.age+z+male)
summary(mod.naive.full)


ans.full<-audit.correct.cov.vdependent.2plus(ystar,cbind(w,w.age),cbind(z,male),v,cbind(x,x.age),y,p.cat,cat,poss.cat)

cbind(round(mod.naive.full$coeff[-1],2),
round(mod.naive.full$coeff[-1]-1.96*summary(mod.naive.full)$coeff[-1,2],2),
round(mod.naive.full$coeff[-1]+1.96*summary(mod.naive.full)$coeff[-1,2],2))

cbind(round(ans.full$CorrectEst,2),
round(ans.full$CorrectEst-1.96*sqrt(diag(ans.full$VarCorrectEst)),2),
round(ans.full$CorrectEst+1.96*sqrt(diag(ans.full$VarCorrectEst)),2))


mod.naive.site<-lm(ystar~w+z)
summary(mod.naive.site)

ans.site<-audit.correct.cov.vdependent.2plus(ystar,w,z,v,x,y,p.cat,cat,poss.cat) ### These aren't exactly the same, even though they should be.
ans.site1<-audit.correct.cov.vdependent(ystar,w,z,v,x,y,p.cat,cat,poss.cat) ### However, they are very close to the same.
ans.site.ignore<-audit.correct.cov(ystar,w,z,v,x,y) ### Not correct because it ignores the audit-specific sampling,
### although results are similar

cbind(round(mod.naive.site$coeff[-1],2),
round(mod.naive.site$coeff[-1]-1.96*summary(mod.naive.site)$coeff[-1,2],2),
round(mod.naive.site$coeff[-1]+1.96*summary(mod.naive.site)$coeff[-1,2],2))

cbind(round(ans.site1$CorrectEst,2),
round(ans.site1$CorrectEst-1.96*sqrt(diag(ans.site1$VarCorrectEst)),2),
round(ans.site1$CorrectEst+1.96*sqrt(diag(ans.site1$VarCorrectEst)),2))

cbind(round(ans.site.ignore$CorrectEst,2),
round(ans.site.ignore$CorrectEst-1.96*sqrt(diag(ans.site.ignore$VarCorrectEst)),2),
round(ans.site.ignore$CorrectEst+1.96*sqrt(diag(ans.site.ignore$VarCorrectEst)),2))


##### Without covariates
N<-length(s)
mod.atten<-lm(ystar~w)
p.est<-mean(s[v==1])
u.est<-w[s*v==1]-x[s*v==1]
var.u<-var(u.est,na.rm=TRUE)#*(sum(s*v,na.rm=TRUE)-1)/sum(s*v,na.rm=TRUE)
var.x<-var(x[v==1],na.rm=TRUE)#*((sum(v)-1)/sum(v))
var.w<-var(w,na.rm=TRUE)
resid.y<-y[v*s==1]-ystar[v*s==1]
resid.x<-x[v*s==1]-w[v*s==1]
cov.est<-cor(resid.y[!is.na(resid.y)&!is.na(resid.x)],resid.x[!is.na(resid.y)&!is.na(resid.x)])*
sqrt(var(resid.y,na.rm=TRUE)*var(resid.x,na.rm=TRUE))
gamma1<-mod.atten$coeff[2]
fixed<-(gamma1*(var.x+var.u*p.est) - p.est*cov.est)/var.x
mu.x<-mean(x[v==1])
mu.w<-mean(w)
mu.star<-mean(resid.y,na.rm=TRUE)
var.betaN<-summary(mod.atten)$coeff[2,2]^2
A<-matrix(0,8,8)
A[1,1]<-1
A[1,2]<-mean(w)
A[2,1]<-mean(w)
A[2,2]<-sum(w^2)/N
A[3,3]<-sum(v)/N
A[4,3]<-2*sum(v*(x-mu.x))/N
A[4,4]<-sum(v)/N
A[5,5]<-sum(v)/N
A[6,6]<-sum(s*v)/N
A[7,7]<-sum(s*v)/N
A[8,7]<-sum(s*v*(w-x))/N
A[8,8]<-sum(s*v)/N
psi<-cbind((ystar-mod.atten$coeff[1]-gamma1*w),
(ystar-mod.atten$coeff[1]-gamma1*w)*w,
(x-mu.x)*v,
((x-mu.x)^2-var.x)*v,
(s-p.est)*v,
((w-x)^2-var.u)*s*v,
(ystar-y-mu.star)*s*v,
((ystar-y-mu.star)*(w-x)-cov.est)*s*v
)
B<-t(psi)%*%psi/N
var.theta<-solve(A)%*%B%*%solve(A)/N
dg<-c(0, (var.x+p.est*var.u)/var.x, 0, -(gamma1*var.u*p.est-p.est*cov.est)/var.x^2,
(gamma1*var.u-cov.est)/var.x, gamma1*p.est/var.x, 0, -p.est/var.x )
var.fixed<-dg%*%var.theta%*%dg
lower.fixed<-fixed-1.96*sqrt(var.fixed)
upper.fixed<-fixed+1.96*sqrt(var.fixed)
lower.atten<-gamma1-1.96*sqrt(var.betaN)
upper.atten<-gamma1+1.96*sqrt(var.betaN)

fixed
lower.fixed
upper.fixed

gamma1
lower.atten
upper.atten