Nothing
RCreliability.in<-function(r,z,W=NULL,Y){
n <- length(r)
t<-length(z)
meanz<-sapply(z, function(x) colMeans(x, na.rm = T))[1,]
sdz<-sapply(z, function(x) apply(x, 2, function(y) sd(y,na.rm = T)))[1,]
z.new <- sapply(1:t, function(x) (z[[x]]-meanz[x])/sdz[x], simplify = F)
if(is.null(W)){
zbar<-sapply(z.new,function(y) rowMeans(y,na.rm = T))#mean prevalence=0.503
#1. naive analysis
fit1<-glm(Y~zbar,family = "binomial")
beta.fit1<-fit1$coefficients #naive estimator
var1<-vcov(fit1) #naive covariance matrix
tab1<-summary(fit1)$coefficients
tab1[,1:2] <- tab1[,1:2]/c(1,sdz)
CI.low<-tab1[,1]-1.96*tab1[,2]
CI.high<-tab1[,1]+1.96*tab1[,2]
tab1<-cbind(tab1,exp(cbind(OR = tab1[, 1],CI.low,CI.high)))
#2. on xhat
muz<-as.numeric((r%*%zbar)/sum(r))
v<-sum(r)-sum(r^2)/sum(r)
dif<-sapply(1:n, function(x) zbar[x,]-muz)
if(t == 1){
sigmazstar<-t(dif)%*%(dif*r)/v
}else{
sigmazstar<-dif%*%t(dif*r)/v
}
diff<-as.matrix(na.omit(sapply(1:t,function(x) (z.new[[x]]-zbar[,x]))))
if(t==1){
sigma<-sum(sapply(1:nrow(diff),function(x) diff[x,]%*%t(diff[x,])))/sum(r-1)
}else{
sigma<-matrix(rowSums(sapply(1:nrow(diff),function(x) diff[x,]%*%t(diff[x,]))),t)/sum(r-1)
}
v12star<-sigmazstar-(n)*sigma/v
sigmax<-sigmazstar-(n)*sigma/v
sigmawithin<-sigma
icc<-sigmax%*%solve(sigmazstar)
if(t==1){
sigmazhat<-t(sapply(r,function(x) sigmazstar-(n)*sigma/v+sigma/x))
}else{
sigmazhat<-sapply(r,function(x) sigmazstar-(n)*sigma/v+sigma/x)
}
if(t==1){
xhat<-sapply(1:n,function(i) (v12star%*%solve(matrix(sigmazhat[,i],ncol=t))%*%zbar[i]))
}else{
xhat<-t(sapply(1:n,function(i) (v12star%*%solve(matrix(sigmazhat[,i],ncol=t))%*%zbar[i,])))
}
fit2<-glm(Y~xhat,family = "binomial")
beta.fit2<-fit2$coefficients
var2<-sandwich(fit2)
tab2<-summary(fit2)$coefficients
tab2[,2]<-sqrt(diag(var2))
tab2[,1:2] <- tab2[,1:2]/c(1,sdz)
CI.low<-tab2[,1]-1.96*tab2[,2]
CI.high<-tab2[,1]+1.96*tab2[,2]
tab2<-cbind(tab2,exp(cbind(OR = tab2[, 1],CI.low,CI.high)))
p<-as.vector(exp(beta.fit2 %*% t(cbind(1,xhat)))/(1+exp(beta.fit2 %*% t(cbind(1,xhat)))))
sigmaz.inv<-solve(sigmazstar)
Z<-t(cbind(zbar))
m<-matrix(0,nrow = t,ncol = t)
c<-matrix(0,nrow = t,ncol = t)
### sigmax
###diag
ddv.x<-sapply(1:t,function(x){
a<-rep(0,t)
a[x]<-a[x]+1
diag(a)
},simplify = F)
ddm.x<-sapply(1:t,function(x){
a<-rep(0,t)
a[x]<-a[x]+1
diag(a)
},simplify = F)
db.x<-sapply(1:t,function(x) t(sapply(1:n,function(y) (ddv.x[[x]]%*%solve(matrix(sigmazhat[,y],ncol=t))-
v12star%*%solve(matrix(sigmazhat[,y],ncol=t))%*%ddm.x[[x]]%*%solve(matrix(sigmazhat[,y],ncol=t)))%*%Z[,y])),simplify = F)
###off-diag
if(t>1){odv.x<-sapply(1:(t-1),function(x) sapply(min((x+1),t):t,function(y){
c[x,y]<-c[x,y]+1
c[y,x]<-c[y,x]+1
c
},simplify = F),simplify = F)
odm.x<-sapply(1:(t-1),function(x) sapply(min((x+1),t):t,function(y){
m[x,y]<-m[x,y]+1
m[y,x]<-m[y,x]+1
m
},simplify = F),simplify = F)
ob.x<-sapply(1:(t-1),function(x) sapply(1:(t-x),function(y)
t(sapply(1:n,function(u) (odv.x[[x]][[y]]%*%solve(matrix(sigmazhat[,u],ncol=t))-
v12star%*%solve(matrix(sigmazhat[,u],ncol=t))%*%odm.x[[x]][[y]]%*%solve(matrix(sigmazhat[,u],ncol=t)))%*%Z[,u])),simplify = F),simplify = F)
}
### sigma
###diag
db.0<-sapply(1:t,function(x) t(sapply(1:n, function(u) t((-ddv.x[[x]]%*%solve(matrix(sigmazhat[,u],ncol=t)))%*%(Z/r)[,u]))),simplify = F)
###off-diag
if(t>1){ob.0<-sapply(1:(t-1),function(x) sapply(1:(t-x),function(y)
t(sapply(1:n,function(u) t((-odv.x[[x]][[y]]%*%solve(matrix(sigmazhat[,u],ncol=t)))%*%(Z/r)[,u]))),simplify = F),simplify = F)}
m<-(t)*(t+1)/2+(t*(t+1))/2 #number of the covariance estimates
s<-t+1 #number of beta estimates
b<-if(t>1){
list(db.x,ob.x,db.0,ob.0)
}else list(db.x,db.0)
b<-sapply(1:m,function(x) matrix(unlist(b),nrow=n)[,(t*x-t+1):(t*x)],simplify = F)
d<-as.matrix(p*(1-p)%*%t(beta.fit2[2:(t+1)]))/n
bstar<-sapply(b,function(x) -rowSums(x*d))
a<-matrix(NA,ncol=m,nrow=s)
a[1,]<-colSums(bstar)
a[2:(t+1),]<-t(apply(as.matrix(xhat), 2, function(x) colSums(x * bstar)))
w = diag(p * (1 - p))
Astar <- t(cbind(1,xhat)) %*% w %*% cbind(1,xhat) /n
A<-rbind(cbind(diag(rep(-1,m)),
matrix(0,nrow=m,ncol=s)),cbind(a,Astar))
A<-rbind(cbind(matrix(diag(rep(-1,t),nrow=t),t),matrix(0,nrow=t,ncol=m+s)),
cbind(matrix(0,nrow=m+s,ncol=t),A))
dmgi<-sapply(1:(t),function(x) r*(((Z[x,])-muz[x]))/sum(r))
###sigmas
###diag
dgi<-r*t(((Z-muz)^2-diag(sigmazstar))/v)
###off-diag
if(t>1){ogi<-sapply(1:(t-1),function(x) sapply((x+1):(t), function(y)
(r*((Z[x,]-muz[x])*((Z[y,]-muz[y])-sigmazstar[x,y]))/v)))}
###within variance
###diag
dgi.0<-if(t==1){
sapply(1:t,function(y) (rowSums((z.new[[y]]-zbar[,y])^2,na.rm = T)-sigma*(r-1))/sum(r-1))
}else sapply(1:t,function(y) (rowSums((z.new[[y]]-zbar[,y])^2,na.rm = T)-diag(sigma)[y]*(r-1))/sum(r-1))
###off-diag
if(t>1){ogi.0<-sapply(1:(t-1),function(x) sapply((x+1):t, function(y)
(rowSums((z.new[[x]]-zbar[,x])*(z.new[[y]]-zbar[,y]),na.rm = T)-sigma[x,y]*(r-1))/sum(r-1)))
}
###betas
gi.beta<-cbind(1,xhat)*(Y-p) / n
gi.1<-if(t==1){
cbind(dmgi,dgi,gi.beta)
}else cbind(dmgi,dgi,matrix(unlist(ogi),nrow=n),gi.beta)
B1<-t(gi.1)%*%gi.1
gi.2<-if(t==1){
dgi.0
}else cbind(dgi.0,matrix(unlist(ogi.0),nrow=n))
B2<-t(gi.2)%*%gi.2
m1<-(t)*(t+1)/2
m2<-t*(t+1)/2
B<-rbind(cbind(B1[1:(t+m1),1:(t+m1)],matrix(0,t+m1,m2),B1[1:(t+m1),(t+m1+1):(t+m1+s)]),
cbind(matrix(0,m2,t+m1),B2,matrix(0,m2,s)),
cbind(B1[(t+m1+1):(t+m1+s),1:(t+m1)],matrix(0,s,m2),B1[(t+m1+1):(t+m1+s),(t+m1+1):(t+m1+s)]))
cov2<-solve(A)%*%B%*%t(solve(A))
tab3<-summary(fit2)$coefficients
tab3[,2]<-sqrt(diag(cov2)[(t+m+1):(t+m+s)])
tab3[,1:2] <- tab3[,1:2]/c(1,sdz)
CI.low<-tab3[,1]-1.96*tab3[,2]
CI.high<-tab3[,1]+1.96*tab3[,2]
tab3[,3]<-tab3[,1]/tab3[,2]
tab3[,4]<-2*pnorm(tab3[,3],lower.tail=FALSE)
tab3<-cbind(tab3,exp(cbind(OR = tab3[, 1],CI.low,CI.high)))
}
else{
W<-matrix(W,n)
t<-length(z)
q<-ncol(W)
meanw<-colMeans(W)
sdw<-apply(W, 2, function(y) sd(y))
W.new <- sapply(1:q, function(x) (W[,x]-meanw[x])/sdw[x])
zbar<-sapply(z.new,function(y) rowMeans(y,na.rm = T))#mean prevalence=0.503
#1. naive analysis
fit1<-glm(Y~zbar+W.new,family = "binomial")
beta.fit1<-fit1$coefficients #naive estimator
var1<-vcov(fit1) #naive covariance matrix
tab1<-summary(fit1)$coefficients
tab1[,1:2] <- tab1[,1:2]/c(1,sdz,sdw)
CI.low<-tab1[,1]-1.96*tab1[,2]
CI.high<-tab1[,1]+1.96*tab1[,2]
tab1<-cbind(tab1,exp(cbind(OR = tab1[, 1],CI.low,CI.high)))
#2. on xhat
muz<-as.numeric((r%*%zbar)/sum(r))
muw<-as.numeric(colSums(W.new)/n)
v<-sum(r)-sum(r^2)/sum(r)
dif<-rbind(sapply(1:n, function(x) zbar[x,]-muz),
sapply(1:n, function(x) W.new[x,]-muw))
sigmazstar<-dif%*%t(dif*r)/v
diff<-as.matrix(na.omit(sapply(1:t,function(x) (z.new[[x]]-zbar[,x]))))
if(t==1){
sigma<-sum(sapply(1:nrow(diff),function(x) diff[x,]%*%t(diff[x,])))/sum(r-1)
}else{
sigma<-matrix(rowSums(sapply(1:nrow(diff),function(x) diff[x,]%*%t(diff[x,]))),t)/sum(r-1)
}
v12star<-sigmazstar[1:t,]-cbind((n)*sigma/v,matrix(0,ncol = q,nrow=t))
sigmax<-sigmazstar-rbind(cbind((n)*sigma/v,matrix(0,ncol = q,nrow=t)),matrix(0,ncol=t+q,nrow=q))
sigmawithin<-rbind(cbind(sigma,matrix(0,ncol = q,nrow=t)),matrix(0,ncol=t+q,nrow=q))
icc<-sigmax%*%solve(sigmazstar)
colnames(sigmax)<-colnames(var1)[-1]
rownames(sigmax)<-rownames(var1)[-1]
colnames(sigmawithin)<-colnames(var1)[-1]
rownames(sigmawithin)<-rownames(var1)[-1]
colnames(icc)<-colnames(var1)[-1]
rownames(icc)<-rownames(var1)[-1]
sigmazhat<-sapply(r,function(x) sigmax+rbind(cbind(sigma/x,matrix(0,ncol = q,nrow=t)),matrix(0,ncol=t+q,nrow=q)))
if(t==1){
xhat<-sapply(1:n,function(i) (v12star%*%solve(matrix(sigmazhat[,i],ncol=t+q))%*%cbind(zbar,W.new)[i,]))
}else{
xhat<-t(sapply(1:n,function(i) (v12star%*%solve(matrix(sigmazhat[,i],ncol=t+q))%*%cbind(zbar,W.new)[i,])))
}
fit2<-glm(Y~xhat+W.new,family = "binomial")
beta.fit2<-fit2$coefficients
var2<-sandwich(fit2)
tab2<-summary(fit2)$coefficients
tab2[,2]<-sqrt(diag(var2))
tab2[,1:2] <- tab2[,1:2]/c(1,sdz,sdw)
CI.low<-tab2[,1]-1.96*tab2[,2]
CI.high<-tab2[,1]+1.96*tab2[,2]
tab2<-cbind(tab2,exp(cbind(OR = tab2[, 1],CI.low,CI.high)))
p<-as.vector(exp(beta.fit2 %*% t(cbind(1,xhat,W.new)))/(1+exp(beta.fit2 %*% t(cbind(1,xhat,W.new)))))
sigmaz.inv<-solve(sigmazstar)
Z<-t(cbind(zbar,W.new))
m<-matrix(0,nrow = t+q,ncol = t+q)
c<-matrix(0,nrow = t,ncol = t+q)
### sigmax
###diag
ddv.x<-sapply(1:t,function(x){
a<-rep(0,t)
a[x]<-a[x]+1
cbind(diag(a),matrix(0,nrow = t,ncol = q))
},simplify = F)
ddm.x<-sapply(1:t,function(x){
a<-rep(0,t+q)
a[x]<-a[x]+1
diag(a)
},simplify = F)
db.x<-sapply(1:t,function(x) t(sapply(1:n,function(y) (ddv.x[[x]]%*%solve(matrix(sigmazhat[,y],ncol=t+q))-
v12star%*%solve(matrix(sigmazhat[,y],ncol=t+q))%*%ddm.x[[x]]%*%solve(matrix(sigmazhat[,y],ncol=t+q)))%*%Z[,y])),simplify = F)
###off-diag
if(t>1){odv.x<-sapply(1:(t-1),function(x) sapply(min((x+1),t):t,function(y){
c[x,y]<-c[x,y]+1
c[y,x]<-c[y,x]+1
c
},simplify = F),simplify = F)
odm.x<-sapply(1:(t-1),function(x) sapply(min((x+1),t):t,function(y){
m[x,y]<-m[x,y]+1
m[y,x]<-m[y,x]+1
m
},simplify = F),simplify = F)
ob.x<-sapply(1:(t-1),function(x) sapply(1:(t-x),function(y)
t(sapply(1:n,function(u) (odv.x[[x]][[y]]%*%solve(matrix(sigmazhat[,u],ncol=t+q))-
v12star%*%solve(matrix(sigmazhat[,u],ncol=t+q))%*%odm.x[[x]][[y]]%*%solve(matrix(sigmazhat[,u],ncol=t+q)))%*%Z[,u])),simplify = F),simplify = F)
}
### sigmaxW
### all off-diag
odv.xw<-sapply(1:t,function(x) sapply(max((x+1),t+1):(t+q),function(y){
c[x,y]<-c[x,y]+1
c
},simplify = F),simplify = F)
odm.xw<-sapply(1:t,function(x) sapply(max((x+1),t+1):(t+q),function(y){
m[x,y]<-m[x,y]+1
m[y,x]<-m[y,x]+1
m
},simplify = F),simplify = F)
ob.xw<-sapply(1:t,function(x) sapply(1:q,function(y)
t(sapply(1:n,function(u) t((odv.xw[[x]][[y]]%*%solve(matrix(sigmazhat[,u],ncol=t+q))-
v12star%*%solve(matrix(sigmazhat[,u],ncol=t+q))%*%odm.xw[[x]][[y]]%*%solve(matrix(sigmazhat[,u],ncol=t+q)))%*%Z[,u]))),simplify = F),simplify = F)
### sigmaW
###diag
ddm.w<-sapply((t+1):(t+q),function(x){
a<-rep(0,t+q)
a[x]<-a[x]+1
diag(a)
},simplify = F)
db.w<-sapply(1:q,function(x)
t(sapply(1:n,function(u) t((-v12star%*%solve(matrix(sigmazhat[,u],ncol=t+q))%*%ddm.w[[x]]%*%solve(matrix(sigmazhat[,u],ncol=t+q)))%*%Z[,u]))),simplify = F)
###off-diag
if(q>1){odm.w<-sapply((t+1):(t+q-1),function(x) sapply(min((x+1),t+q):(t+q),function(y){
m[x,y]<-m[x,y]+1
m[y,x]<-m[y,x]+1
m
},simplify = F),simplify = F)
ob.w<-sapply(1:(q-1),function(x) sapply(1:(q-x),function(y)
t(sapply(1:n,function(u) (-v12star%*%solve(matrix(sigmazhat[,u],ncol=t+q))%*%odm.w[[x]][[y]]%*%solve(matrix(sigmazhat[,u],ncol=t+q)))%*%Z[,u])),simplify = F),simplify = F)
}
### sigma
###diag
db.0<-sapply(1:t,function(x) t(sapply(1:n, function(u) t((-ddv.x[[x]]%*%solve(matrix(sigmazhat[,u],ncol=t+q)))%*%(Z/r)[,u]))),simplify = F)
###off-diag
if(t>1){ob.0<-sapply(1:(t-1),function(x) sapply(1:(t-x),function(y)
t(sapply(1:n,function(u) t((-odv.x[[x]][[y]]%*%solve(matrix(sigmazhat[,u],ncol=t+q)))%*%(Z/r)[,u]))),simplify = F),simplify = F)}
m<-(t+q)*(t+q+1)/2+(t*(t+1))/2 #number of the covariance estimates
s<-t+q+1 #number of beta estimates
b<-if(t>1 & q>1){
list(db.x,db.w,ob.x,ob.xw,ob.w,db.0,ob.0)
}else if(q>1){
list(db.x,db.w,ob.xw,ob.w,db.0)
}else if(t>1){
list(db.x,db.w,ob.x,ob.xw,db.0,ob.0)
}else list(db.x,db.w,ob.xw,db.0)
b<-sapply(1:m,function(x) matrix(unlist(b),nrow=n)[,(t*x-t+1):(t*x)],simplify = F)
d<-as.matrix(p*(1-p)%*%t(beta.fit2[2:(t+1)]))/n
bstar<-sapply(b,function(x) -rowSums(x*d))
a<-matrix(NA,ncol=m,nrow=s)
a[1,]<-colSums(bstar)
a[2:(t+1),]<-t(apply(as.matrix(xhat), 2, function(x) colSums(x * bstar)))
a[(t+2):(t+q+1),]<-t(apply(as.matrix(W.new), 2, function(x) colSums(x * bstar)))
w <- diag(p * (1 - p))
Astar <- t(cbind(1,xhat, W.new)) %*% w %*% cbind(1,xhat, W.new)/n
A<-rbind(cbind(diag(rep(-1,m)),
matrix(0,nrow=m,ncol=s)),cbind(a,Astar))
A<-rbind(cbind(diag(rep(-1,t+q),nrow=t+q),matrix(0,nrow=t+q,ncol=m+s)),
cbind(matrix(0,nrow=m+s,ncol=t+q),A))
###sigmas
###diag
mu<-c(muz,muw)
dmgi<-sapply(1:(t+q),function(x) r*(Z[x,]-mu[x])/sum(r))
dgi<-sapply(1:(t+q),function(x) (r*((Z[x,]-mu[x])^2-diag(sigmazstar)[x]))/v)
###off-diag
ogi<-sapply(1:(t+q-1),function(x) sapply((x+1):(t+q), function(y)
(r*((Z[x,]-mu[x])*(Z[y,]-mu[y])-sigmazstar[x,y]))/v))
###within variance
###diag
dgi.0<-if(t==1){
sapply(1:t,function(y) (rowSums((z.new[[y]]-zbar[,y])^2,na.rm = T)-sigma*(r-1))/sum(r-1))
}else sapply(1:t,function(y) (rowSums((z.new[[y]]-zbar[,y])^2,na.rm = T)-diag(sigma)[y]*(r-1))/sum(r-1))
###off-diag
if(t>1){ogi.0<-sapply(1:(t-1),function(x) sapply((x+1):t, function(y)
(rowSums((z.new[[x]]-zbar[,x])*(z.new[[y]]-zbar[,y]),na.rm = T)-sigma[x,y]*(r-1))/sum(r-1)))
}
###betas
gi.beta<-cbind(1,xhat,W.new)*(Y-p)/n
gi.1<-cbind(dmgi,dgi,matrix(unlist(ogi),nrow=n),gi.beta)
B1<-crossprod(gi.1)
gi.2<-if(t==1){
dgi.0
}else cbind(dgi.0,matrix(unlist(ogi.0),nrow=n))
B2<-t(gi.2)%*%gi.2
m1<-(t+q)*(t+q+1)/2
m2<-t*(t+1)/2
B<-rbind(cbind(B1[1:(t+q+m1),1:(t+q+m1)],matrix(0,t+q+m1,m2),B1[1:(t+q+m1),(t+q+m1+1):(t+q+m1+s)]),
cbind(matrix(0,m2,t+q+m1),B2,matrix(0,m2,s)),
cbind(B1[(t+q+m1+1):(t+q+m1+s),1:(t+q+m1)],matrix(0,s,m2),B1[(t+q+m1+1):(t+q+m1+s),(t+q+m1+1):(t+q+m1+s)]))
cov2<-solve(A)%*%B%*%t(solve(A))
CI.low<-beta.fit2-1.96*sqrt(diag(cov2)[(t+q+m+1):(t+q+m+s)])
CI.high<-beta.fit2+1.96*sqrt(diag(cov2)[(t+q+m+1):(t+q+m+s)])
coverage.xhat0<-log(1.5)<=CI.high & log(1.5)>=CI.low
tab3<-summary(fit2)$coefficients
tab3[,2]<-sqrt(diag(cov2)[(t+q+m+1):(t+q+m+s)])
tab3[,1:2] <- tab3[,1:2]/c(1,sdz,sdw)
CI.low<-tab3[,1]-1.96*tab3[,2]
CI.high<-tab3[,1]+1.96*tab3[,2]
tab3[,3]<-tab3[,1]/tab3[,2]
tab3[,4]<-2*pnorm(tab3[,3],lower.tail=FALSE)
tab3<-cbind(tab3,exp(cbind(OR = tab3[, 1],CI.low,CI.high)))
}
list("Naive estimates"=tab1[-1,],
"Corrected estimates"=tab2[-1,],
"Corrected estimates, taking into account the extra variation due to estimating the parameters in the measurement error model"=tab3[-1,])
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.