Subscribe for X-meta Updates!

© 2017 by XMETA. Created using Wix.com

Visits since 08-25-2017

## CODE FOR MMETA FUNCTION

Please choose your preferred code according to the programming language you use

## STATA CODE FOR MMOM METHOD IN MMETA FUNCTION

capture program drop mmom

program mmom, rclass

version 1.1

args ys1 ys2 vars1 vars2

** ys1 ys2 are the two measures from the same study

** vars1 and vars2 are the associated variance

tempvar w1 w2

gen `w1' = 1/vars1

gen `w2' = 1/vars2

tempvar w1_ys1 w2_ys2

gen `w1_ys1' = `w1'*ys1

summ (`w1_ys1'), meanonly

scalar sum_w1_ys1 = r(sum)

gen `w2_ys2' = `w2'*ys2

summ (`w2_ys2'), meanonly

scalar sum_w2_ys2 = r(sum)

summ (`w1'), meanonly

scalar sum_w1 = r(sum)

summ (`w2'), meanonly

scalar sum_w2 = r(sum)

scalar y1_weight = sum_w1_ys1/sum_w1

scalar y2_weight = sum_w2_ys2/sum_w2

gen large_vars1 = 1-1*(vars1 > 10^4)

gen large_vars2 = 1-1*(vars2 > 10^4)

summ (`large_vars1'), meanonly

scalar n1 = r(sum)

summ (`large_vars2'), meanonly

scalar n2 = r(sum)

gen inter1 = `w1'*(ys1-y1_weight)^2

gen inter2 = `w2'*(ys2-y2_weight)^2

summ (inter1), meanonly

scalar Q1 = r(sum)

summ (inter2), meanonly

scalar Q2 = r(sum)

tempvar w1_2 w2_2

gen `w1_2' = `w1'^2

gen `w2_2' = `w2'^2

summ (`w1_2'), meanonly

scalar sum_w1_2 = r(sum)

summ (`w2_2'), meanonly

scalar sum_w2_2 = r(sum)

scalar tau1_2_hat = max(0, (Q1-(n1-1))/(sum_w1-sum_w1_2/sum_w1))

scalar tau2_2_hat = max(0, (Q2-(n2-1))/(sum_w2-sum_w2_2/sum_w2))

* variance estimate:

gen w1_star = 1/(vars1 + tau1_2_hat)

gen w2_star = 1/(vars2 + tau2_2_hat)

summ (w1_star), meanonly

scalar sum_w1_star = r(sum)

summ (w2_star), meanonly

scalar sum_w2_star = r(sum)

gen w1_star_ys1 = w1_star * ys1

gen w2_star_ys2 = w2_star * ys2

summ (w1_star_ys1), meanonly

scalar sum_w1_star_ys1 = r(sum)

summ (w2_star_ys2), meanonly

scalar sum_w2_star_ys2 = r(sum)

scalar beta1_hat = sum_w1_star_ys1/sum_w1_star

scalar beta2_hat = sum_w2_star_ys2/sum_w2_star

scalar var_beta1_hat = 1/sum_w1_star

scalar var_beta2_hat = 1/sum_w2_star

gen mycov_beta_piece = (w1_star/sum_w1_star)*(w2_star/sum_w2_star)*(ys1 - beta1_hat)*(ys2 - beta2_hat)

summ (mycov_beta_piece), meanonly

scalar mycov_beta = r(sum)

matrix beta_hat = (beta1_hat , beta2_hat)

matrix sigma_hat = (var_beta1_hat, mycov_beta \ mycov_beta,var_beta2_hat)

matrix list beta_hat

matrix list sigma_hat

end

## R CODE FOR MMETA FUNCTION

mmeta <-

function(data, rhow, type, k, method) {

nn.reml <-

function(y1, s1, y2, s2, rhow) {

S=cbind(s1^2, rhow*s1*s2, s2^2)

myreml= mvmeta(cbind(y1, y2), S, method="reml")

coef = coef(myreml)[1,]

vcov = vcov(myreml)

rhob <- cov2cor(myreml\$Psi)[2,1]

myresults = list(coefficients = coef, vcov= vcov, rhob=rhob)

return(myresults)

}

nn.cl <-

function(y1, s1, y2, s2) {

m=length(y1)

estim.pseudo=rep(NA, 4)

my.pseudo1= mvmeta(y1,s1^2,method="reml")

my.pseudo2= mvmeta(y2,s2^2,method="reml")

estim.pseudo[c(1:2)]=c(coef(my.pseudo1), my.pseudo1\$Psi)

estim.pseudo[c(3:4)]=c(coef(my.pseudo2), my.pseudo2\$Psi)

estim.pseudo = matrix(estim.pseudo, nrow=1)

colnames(estim.pseudo)=c("beta1", "tau1.2", "beta2", "tau2.2")

Score1.beta = (y1-coef(my.pseudo1))/(s1^2+my.pseudo1\$Psi)

Score2.beta = (y2-coef(my.pseudo2))/(s2^2+my.pseudo2\$Psi)

Score1.tau2 = -1/(2*(s1^2+my.pseudo1\$Psi)) + (y1-coef(my.pseudo1))^2/(s1^2+my.pseudo1\$Psi)^2/2

Score2.tau2 = -1/(2*(s2^2+my.pseudo2\$Psi)) + (y2-coef(my.pseudo2))^2/(s2^2+my.pseudo2\$Psi)^2/2

Score1 = rbind(Score1.beta, Score1.tau2); Score2 = rbind(Score2.beta, Score2.tau2)

I11.hat = Score1%*%t(Score1)/m; I22.hat = Score2%*%t(Score2)/m

I12.hat = Score1%*%t(Score2)/m

myoff.diag = solve(I11.hat,tol=1e-100)%*%I12.hat%*%solve(I22.hat,tol=1e-100)/m

myupper = cbind(solve(m*I11.hat,tol=1e-100), (myoff.diag))

mylower = cbind(t(myoff.diag), solve(m*I22.hat,tol=1e-100))

myV=rbind(myupper,mylower)

colnames(myV)=c("beta1", "tau1.2", "beta2", "tau2.2")

rownames(myV)=c("beta1", "tau1.2", "beta2", "tau2.2")

myresults = list(coefficients = estim.pseudo, vcov= myV)

return(myresults)

}

nn.mom= function(y1, s1, y2, s2){

v1 = s1^2

v2 = s2^2

w1 = 1/(v1)

w2 = 1/(v2)

y1.weight = sum(w1*y1)/sum(w1)

y2.weight = sum(w2*y2)/sum(w2)

n1 = sum(1-1*(v1 > 10^4))

n2 = sum(1-1*(v2 > 10^4))

Q1 = sum(w1*(y1-y1.weight)^2)

Q2 = sum(w2*(y2-y2.weight)^2)

tau1.2.hat = max(0, (Q1-(n1-1))/(sum(w1)-sum(w1^2)/sum(w1)))

tau2.2.hat = max(0, (Q2-(n2-1))/(sum(w2)-sum(w2^2)/sum(w2)))

## variance estimate:

w1.star = 1/(v1 + tau1.2.hat)

w2.star = 1/(v2 + tau2.2.hat)

beta1.hat = sum(w1.star*y1)/sum(w1.star)

beta2.hat = sum(w2.star*y2)/sum(w2.star)

var.beta1.hat = 1/sum(w1.star)

var.beta2.hat = 1/sum(w2.star)

mycov.beta = sum((w1.star/sum(w1.star))*(w2.star/sum(w2.star))*(y1 - beta1.hat)*(y2 - beta2.hat))

beta.hat = cbind(beta1=beta1.hat,beta2=beta2.hat)

myV = matrix(c(var.beta1.hat,mycov.beta,mycov.beta,var.beta2.hat),nrow = 2, byrow = T)

colnames(myV)=c("beta1", "beta2")

rownames(myV)=c("beta1", "beta2")

myresults = list(coefficients = beta.hat, vcov= myV)

return(myresults)

}

log.Riley.Lik.reml.func = function(mygamma, y1, s1, y2, s2){

m = length(y1)

beta1 = mygamma[1]; beta2= mygamma[2]

log.psi1.2 = mygamma[3]; log.psi2.2 = mygamma[4]; omega = mygamma[5]

psi1.2 = exp(log.psi1.2); psi2.2 = exp(log.psi2.2)

rho = 2*plogis(omega)-1

s1.2 = s1^2

s2.2 = s2^2

Sigma = matrix(0, nrow=(2*m), ncol=(2*m))

for(i in 1:m){

cov.y1.y2 = rho*sqrt((psi1.2+s1.2[i])*(psi2.2+s2.2[i]))

Sigma.i = matrix(c(s1.2[i]+psi1.2, cov.y1.y2, cov.y1.y2, s2.2[i]+psi2.2), nrow=2)

Sigma[((i-1)*2+c(1:2)), ((i-1)*2+c(1:2))] = Sigma.i

}

y1.temp = rep(y1, each=2); y2.temp = rep(y2, each=2)

y1.temp2 = y1.temp%*%diag(rep(c(1,0), times=m))

y2.temp2 = y2.temp%*%diag(rep(c(0,1), times=m))

Y = as.vector(y1.temp2+y2.temp2)

x.design=cbind(rep(c(1,0), times=m), rep(c(0,1),times=m))

mybeta=x.design%*%c(beta1,beta2)

mylik = -0.5*((m-2)*log(2*pi)-log(det(t(x.design)%*%x.design)) +

log(det(Sigma))+log(det(t(x.design)%*%solve(Sigma,tol=1e-100)%*%x.design)) +

t(Y-mybeta)%*%solve(Sigma,tol=1e-100)%*%(Y-mybeta))

mylik = as.numeric(mylik)

return(mylik)

}

log.Riley.Lik.reml.i.func=function(mygamma, y1.i, s1.i, y2.i, s2.i){

result=log.Riley.Lik.reml.func(mygamma, y1.i, s1.i, y2.i, s2.i)

return(result)

}

robustV.ftn=function(mygamma, y1, s1, y2, s2){

m = length(y1)

SigmaB=array(NA, c(5,5,m))

SigmaA=array(NA, c(5,5,m))

for (i in 1:m){

y1.i=y1[i]

s1.i=s1[i]

y2.i=y2[i]

s2.i=s2[i]

myB=jacobian(log.Riley.Lik.reml.i.func, mygamma, y1.i=y1.i, s1.i=s1.i, y2.i=y2.i, s2.i=s2.i)

SigmaB[,,i]=t(myB)%*%myB

myA=hessian(log.Riley.Lik.reml.i.func, mygamma, y1.i=y1.i, s1.i=s1.i, y2.i=y2.i, s2.i=s2.i)

SigmaA[,,i]=-myA

}

E.A=apply(SigmaA, c(1:2), mean) # sensitivity matrix

E.B=apply(SigmaB, c(1:2), mean) # variability matrix

myV=solve(E.A)%*%E.B%*%solve(E.A)

return(myV)

}

nn.rs=function(y1, s1, y2, s2){

mygamma.init = c(0,0,0,0,0)

myresults = optim(mygamma.init, log.Riley.Lik.reml.func, y1=y1, s1=s1, y2=y2, s2=s2, hessian=TRUE,

control = list(fnscale=-1,maxit=1000))

beta1=myresults\$par[1]

tau1.2=exp(myresults\$par[3])

beta2=myresults\$par[2]

tau2.2=exp(myresults\$par[4])

rhos=2*plogis(myresults\$par[5])-1

par=cbind(beta1, tau1.2, beta2, tau2.2, rhos)

myV=robustV.ftn(myresults\$par, y1=y1, s1=s1, y2=y2, s2=s2)/length(y1)

colnames(myV)=c("beta1", "log.tau1.2", "beta2", "log.tau2.2", "omega")

rownames(myV)=c("beta1", "log.tau1.2", "beta2", "log.tau2.2", "omega")

return(list(coefficients=par, vcov=myV))

}

myLik.indep.log=function(mypar, mydat) {

a1.temp <- mypar[1]; b1.temp <- mypar[2]

a2.temp <- mypar[3]; b2.temp <- mypar[4]

a1 <- exp(a1.temp); b1 <- exp(b1.temp)

a2 <- exp(a2.temp); b2 <- exp(b2.temp)

temp1 <- (lgamma(a1+mydat\$y1) + lgamma(b1+mydat\$n1-mydat\$y1)

+ lgamma(a2+mydat\$y2) + lgamma(b2+mydat\$n2-mydat\$y2)

+ lgamma(a1+b1) + lgamma(a2+b2))

temp2 <- (lgamma(a1) + lgamma(b1) + lgamma(a2) + lgamma(b2)

+ lgamma(a1+b1+mydat\$n1) + lgamma(a2+b2+mydat\$n2))

myLogLik <- sum(temp1 - temp2)

return(myLogLik)

}

par.cal=function(mypar) {

a1 <- exp(mypar[1]); b1 <- exp(mypar[2])

a2 <- exp(mypar[3]); b2 <- exp(mypar[4])

eta <- mypar[5]

cc <- sqrt(a1*a2*b1*b2)/sqrt((a1+b1+1)*(a2+b2+1))

upper.bound <- cc/max(a1*b2, a2*b1)

lower.bound <- -cc/max(a1*a2, b1*b2)

expit.eta= exp(eta)/(1+exp(eta))

rho <- (upper.bound-lower.bound)*expit.eta + lower.bound

return(c(a1,b1,a2,b2,rho,eta))

}

bb.cl=function(y1,n1,y2,n2){

nstudy=length(y1)

initial.val.gen <- function(y, n) {

BBfit <- betabin(cbind(y, n-y)~1, ~1, data=data.frame(y=y,n=n))

expit.BB <- exp(as.numeric(BBfit@param[1]))/(1+exp(as.numeric(BBfit@param[1])))

a.ini <- expit.BB*(1/as.numeric(BBfit@param[2])-1)

b.ini <- (1/as.numeric(BBfit@param[2])-1)*(1-expit.BB)

return(list(a=a.ini, b=b.ini))

}

mle.CL <- function(y1=y1,n1=n1,y2=y2,n2=n2) {

init.val <- rep(0, 5)

fit1 <- initial.val.gen(y1, n1)

fit2 <- initial.val.gen(y2, n2)

init.val[1] <- log(fit1\$a);

init.val[2] <- log(fit1\$b)

init.val[3] <- log(fit2\$a);

init.val[4] <- log(fit2\$b)

MLE.inde.log <- optim(init.val[1:4], myLik.indep.log, method = "L-BFGS-B",

lower=rep(-20,4), upper=rep(20,4),

control = list(fnscale=-1,maxit=1000),

hessian = T, mydat=list(y1=y1,n1=n1,y2=y2,n2=n2))

mypar<- par.cal(MLE.inde.log\$par)

rho<-0;

eta<-NA;

hessian.log<-MLE.inde.log\$hessian

colnames(hessian.log)<-c("loga1","logb1","loga2","logb2")

rownames(hessian.log)<-c("loga1","logb1","loga2","logb2")

conv=MLE.inde.log\$convergence

a1 <- mypar[1]; b1 <- mypar[2];

a2 <- mypar[3]; b2 <- mypar[4];

prior.MLE<-c(a1, b1, a2, b2, rho ,eta)

return(list(MLE=prior.MLE,hessian.log=hessian.log,conv=conv))

}

myLik.indep.vector=function(mypar, mydat) {

a1.temp <- mypar[1]; b1.temp <- mypar[2]

a2.temp <- mypar[3]; b2.temp <- mypar[4]

a1 <- exp(a1.temp); b1 <- exp(b1.temp)

a2 <- exp(a2.temp); b2 <- exp(b2.temp)

temp1 <- (lgamma(a1+mydat\$y1) + lgamma(b1+mydat\$n1-mydat\$y1)

+ lgamma(a2+mydat\$y2) + lgamma(b2+mydat\$n2-mydat\$y2)

+ lgamma(a1+b1) + lgamma(a2+b2))

temp2 <- (lgamma(a1) + lgamma(b1) + lgamma(a2) + lgamma(b2)

+ lgamma(a1+b1+mydat\$n1) + lgamma(a2+b2+mydat\$n2))

myLogLik <- (temp1 - temp2)

return(myLogLik)

}

sandwich.var.cal=function(mypar, myhessian.indep, mydat){

npar=length(mypar)

delta=1e-8

myD = matrix(NA, nrow=npar, ncol=nstudy)

for(i in 1:npar){

par.for=par.back=mypar

par.for[i]=par.for[i]+delta

par.back[i]=par.back[i]-delta

myD[i,] = (myLik.indep.vector(par.for, mydat)-myLik.indep.vector(par.back,mydat))/(2*delta)

}

score.squared = myD%*%t(myD)

inv.hessian = solve(-myhessian.indep)

results=inv.hessian%*% score.squared %*%inv.hessian

return(results)

}

out=mle.CL(y1=y1,n1=n1,y2=y2,n2=n2)

mle=out\$MLE

hessian.log=out\$hessian.log

a1=mle[1];b1=mle[2]

a2=mle[3];b2=mle[4]

mypar=log(c(a1,b1,a2,b2))

mydat=list(y1=y1,n1=n1,y2=y2,n2=n2)

covar.sandwich=sandwich.var.cal(mypar, hessian.log, mydat)

myD <- matrix(c(-1, 1, 1, -1), nrow=1)

myVar.log <- as.numeric(myD %*% covar.sandwich %*% t(myD))

logOR <- log((a2/b2)/(a1/b1))

OR_CI=c(exp(logOR-1.96*sqrt(myVar.log)), exp(logOR+1.96*sqrt(myVar.log)))

#return(list(logOR=logOR, OR_CI=OR_CI, se=sqrt(myVar.log),hessian.log=hessian.log,conv=out\$conv,mle=mle))

return(list(logOR=logOR, OR_CI=OR_CI, se=sqrt(myVar.log),hessian.log=hessian.log,conv=out\$conv,mle=mle))

}

logit = function(p) log(p/(1-p))

expit = function(b) exp(b)/(1+exp(b))

dlogitnorm = function(p, mu, tau2){

(2*pi*tau2)^(-1/2)*exp(-(qlogis(p)-mu)^2/(2*tau2))/(p*(1-p))

}

integrand1<-function(n,y, mypar2, p.grid){

dbinom(y, size=n, prob=p.grid)*dlogitnorm(p.grid, mypar2[1], mypar2[2])

}

integrand2<-function(n,y,mypar2, p.grid){

dbinom(y, size=n, prob=p.grid)*dlogitnorm(p.grid, mypar2[1], mypar2[2])*2*(qlogis(p.grid)-mypar2[1])/(2*mypar2[2])

}

integrand3<-function(n,y,mypar2, p.grid){

dbinom(y, size=n, prob=p.grid)*dlogitnorm(p.grid, mypar2[1], mypar2[2])*(-1/(2*mypar2[2])+(qlogis(p.grid)-mypar2[1])^2/(2*mypar2[2]^2))

}

mystandize = function(MM){

diag(1/sqrt(diag(MM)))%*%MM%*%diag((1/sqrt(diag(MM))))

}

bn.cl=function(y1, n1, y2, n2){

id=c(1:length(y1))

mydat1 = data.frame(id,n1,y1,n2,y2)

names(mydat1) = c("id","Nd","SeY","Nn","SpY")

m = nrow(mydat1)

estim = estim.orig = mbse.orig = matrix(NA, nrow=4, ncol=1)

mySandwich = matrix(NA,nrow=4,ncol=4)

fit.SeY = glmmML(cbind(SeY, Nd-SeY)~1, data = mydat1, cluster = id, method= "ghq", n.points=8,     control=list(epsilon=1e-08 , maxit=50000,trace=FALSE))

fit.SpY = glmmML(cbind(SpY, Nn-SpY)~1, data = mydat1, cluster = id, method= "ghq", n.points=8, control=list(epsilon=1e-08 , maxit=50000,trace=FALSE))

estim[c(1: 2), 1] = c(as.numeric(fit.SeY\$coefficients), (fit.SeY\$sigma)^2)

estim[c(3: 4), 1] = c(as.numeric(fit.SpY\$coefficients), (fit.SpY\$sigma)^2)

estim.orig[c(1: 2), 1] = c(plogis(estim[1, 1]), estim[2, 1])

estim.orig[c(3: 4), 1] = c(plogis(estim[3, 1]), estim[4, 1])

rownames(estim)=c("SeY", "SeY.S.2", "SpY", "SpY.S.2")

rownames(estim.orig)=c("SeY", "SeY.S.2", "SpY", "SpY.S.2")

I11.hat = solve(m*fit.SeY\$variance)

I22.hat = solve(m*fit.SpY\$variance)

int1 = int2  = int3 = rep(NA, length=m)

est1 = estim[c(1: 2), 1] ##plug in the point estimate based on sensitivity

for(i in 1: m){

int1[i] = integrate(integrand1, lower=0, upper=1, n=mydat1\$Nd[i], mypar2=est1, y=mydat1\$SeY[i])\$value

int2[i] = integrate(integrand2, lower=0, upper=1, n=mydat1\$Nd[i], mypar2=est1, y=mydat1\$SeY[i])\$value

int3[i] = integrate(integrand3, lower=0, upper=1, n=mydat1\$Nd[i], mypar2=est1, y=mydat1\$SeY[i])\$value

}

deriveSeY.mu = int2/int1

deriveSeY.tau2 = int3/int1

B.SeY = cbind(deriveSeY.mu, deriveSeY.tau2)

int1 = int2  = int3 = rep(NA, length=m)

est2 = estim[c(3: 4), 1] ##plug in the point estimate based on specificity

for(i in 1: m){

int1[i] = integrate(integrand1, lower=0, upper=1, n=mydat1\$Nn[i], mypar2=est2, y=mydat1\$SpY[i])\$value

int2[i] = integrate(integrand2, lower=0, upper=1, n=mydat1\$Nn[i], mypar2=est2, y=mydat1\$SpY[i])\$value

int3[i] = integrate(integrand3, lower=0, upper=1, n=mydat1\$Nn[i], mypar2=est2, y=mydat1\$SpY[i])\$value

}

deriveSpY.mu = int2/int1

deriveSpY.tau2 = int3/int1

B.SpY = cbind(deriveSpY.mu, deriveSpY.tau2)

I12.hat = t(B.SeY)%*%B.SpY/m

myoff.diag = solve(I11.hat)%*%I12.hat%*%solve(I22.hat)/m

myupper = cbind(fit.SeY\$variance, myoff.diag)

mylower = cbind(t(myoff.diag), fit.SpY\$variance)

myV = rbind(myupper, mylower)

colnames(myV)=c("SeY", "SeY.S.2", "SpY", "SpY.S.2")

rownames(myV)=c("SeY", "SeY.S.2", "SpY", "SpY.S.2")

return(list(coefficients=estim, par.orig=estim.orig, vcov=myV))

}

score.Li.func = function(PiY,SeY,SpY,n,n1,n0,PiY.est,SeY.est,SpY.est){

beta0 = PiY.est[1]; beta1 = SeY.est[1]; beta2 = SpY.est[1]

mu0 = expit(beta0); mu1 = expit(beta1); mu2 = expit(beta2)

phi0 = PiY.est[2]; phi1 = SeY.est[2]; phi2 = SpY.est[2]

theta0 = phi0/(1-phi0); theta1 = phi1/(1-phi1); theta2 = phi2/(1-phi2)

k01 = c(0:(PiY-1)); k02 = c(0:(n-PiY-1)); k03 = c(0:(n-1))

k11 = c(0:(SeY-1)); k12 = c(0:(n1-SeY-1)); k13 = c(0:(n1-1))

k21 = c(0:(SpY-1)); k22 = c(0:(n0-SpY-1)); k23 = c(0:(n0-1))

# notice for outbound

if(PiY == 0) k01 = NULL; if(PiY == n) k02 = NULL

if(SeY == 0) k11 = NULL; if(SeY == n1) k12 = NULL

if(SpY == 0) k21 = NULL; if(SpY == n0) k22 = NULL

Dbeta0 = Dbeta1 = Dbeta2 = Dphi0 = Dphi1 = Dphi2 = 0

Dphi0 = (sum(k01/(mu0+k01*theta0))+sum(k02/(1-mu0+k02*theta0))- sum(k03/(1+k03*theta0)))/((1-phi0)^2)

Dphi1 = (sum(k11/(mu1+k11*theta1))+sum(k12/(1-mu1+k12*theta1))-sum(k13/(1+k13*theta1)))/((1-phi1)^2)

Dphi2 = (sum(k21/(mu2+k21*theta2))+sum(k22/(1-mu2+k22*theta2))-sum(k23/(1+k23*theta2)))/((1-phi2)^2)

Dbeta0 = sum(mu0*(1-mu0)/(mu0+k01*theta0))-sum(mu0*(1-mu0)/(1-mu0+k02*theta0))

Dbeta1 = sum(mu1*(1-mu1)/(mu1+k11*theta1))-sum(mu1*(1-mu1)/(1-mu1+k12*theta1))

Dbeta2 = sum(mu2*(1-mu2)/(mu2+k21*theta2))-sum(mu2*(1-mu2)/(1-mu2+k22*theta2))

score.PiY = as.vector(c(Dbeta0,Dphi0))

score.SeY = as.vector(c(Dbeta1,Dphi1))

score.SpY = as.vector(c(Dbeta2,Dphi2))

return(list(score.PiY=score.PiY, score.SeY=score.SeY, score.SpY=score.SpY))

}

score.Li.all.func= function(SeY,SpY,n1,n0,SeY.est,SpY.est){

beta1 = SeY.est[1] ; beta2 = SpY.est[1]

mu1 = expit(beta1) ; mu2 = expit(beta2)

phi1 = SeY.est[2] ; phi2 =SpY.est[2]

theta1 = phi1/(1-phi1) ; theta2 = phi2/(1-phi2)

k11 = c(0:(SeY-1)); k12 = c(0:(n1-SeY-1)); k13 = c(0:(n1-1))

k21 = c(0:(SpY-1)); k22 = c(0:(n0-SpY-1)); k23 = c(0:(n0-1))

# notice for outbound

if(SeY == 0) k11 = NULL ; if(SeY == n1) k12 = NULL

if(SpY == 0) k21 = NULL ; if(SpY == n0) k22 = NULL

Dbeta1 = Dbeta2 = Dphi1 = Dphi2 = 0

Dphi1 = (sum(k11/(mu1+k11*theta1))+sum(k12/(1-mu1+k12*theta1))- sum(k13/(1+k13*theta1)))/((1-phi1)^2)

Dphi2 = (sum(k21/(mu2+k21*theta2))+sum(k22/(1-mu2+k22*theta2))- sum(k23/(1+k23*theta2)))/((1-phi2)^2)

Dbeta1 = sum(mu1*(1-mu1)/(mu1+k11*theta1))- sum(mu1*(1-mu1)/(1-mu1+k12*theta1))

Dbeta2 = sum(mu2*(1-mu2)/(mu2+k21*theta2))-sum(mu2*(1-mu2)/(1-mu2+k22*theta2))

score.SeY = as.vector(c(Dbeta1,Dphi1))

score.SpY = as.vector(c(Dbeta2,Dphi2))

return(list(score.SeY=score.SeY, score.SpY=score.SpY))

}

Fisher.inf = function(PiY.del,SeY.del,SpY.del,n.del,n1.del, n0.del, SeY,SpY,n1,n0,PiY.est,SeY.est, SpY.est){

m = length(SeY) # the length of all study, including case control study & cohort study

m2 = length(PiY.del) # the length of the cohort study only

I00.array = I01.array = I02.array = array(NA,c(2,2,m2))

I11.array = I12.array = I22.array = array(NA,c(2,2,m))

for(i in 1:m2){

temp = score.Li.func (PiY.del[i], SeY.del[i],SpY.del[i],n.del[i], n1.del[i],n0.del[i],PiY.est, SeY.est, SpY.est)

score.PiY = temp\$score.PiY

score.SeY = temp\$score.SeY

score.SpY = temp\$score.SpY

I00.array[ , , i] = score.PiY%*%t(score.PiY)

I01.array[ , , i] = score.PiY%*%t(score.SeY)

I02.array[ , , i] = score.PiY%*%t(score.SpY)

}

for(i in 1:m){

temp = score.Li.all.func (SeY[i],SpY[i],n1[i],n0[i],SeY.est,SpY.est)

score.SeY = temp\$score.SeY

score.SpY = temp\$score.SpY

I11.array[ , , i] = score.SeY%*%t(score.SeY)

I12.array[ , , i] = score.SeY%*%t(score.SpY)

I22.array[ , , i] = score.SpY%*%t(score.SpY)

}

I00 = apply(I00.array, c(1,2), mean)

I01 = apply(I01.array, c(1,2), mean)

I02 = apply(I02.array, c(1,2), mean)

I11 = apply(I11.array, c(1,2), mean)

I12 = apply(I12.array, c(1,2), mean)

I22 = apply(I22.array, c(1,2), mean)

return(list(I00=I00, I01=I01, I02=I02, I11=I11, I12=I12, I22=I22))

}

tb.cl=function(n, PiY, SeY, n1, SpY, n2){

m = length(SeY) # the length of all study, including case control study & cohort study

m2 = sum(is.na(PiY)!=1) # the length of the cohort study only

estim = estim.orig = mbse.orig = matrix(NA,nrow=6,ncol=1)

mySandwich = matrix(NA,nrow=6,ncol=6)

mydat1 = data.frame(n,PiY,n1,SeY,n2, SpY)

names(mydat1) = c("n","PiY", "n1","SeY","n0","SpY")

mydat2 = mydat1[is.na(mydat1\$PiY)==FALSE, ] ## dataset used for estimating prevalence

fit.PiY = betabin(cbind(PiY, n-PiY)~1, ~1, hessian=T, data=mydat2, link="logit")

fit.SeY = betabin(cbind(SeY, n1-SeY)~1, ~1, hessian=T, data=mydat1, link="logit")

fit.SpY = betabin(cbind(SpY, n0-SpY)~1, ~1, hessian=T, data=mydat1, link="logit")

estim[c(1:2), 1] = c(as.numeric(fit.PiY@param))

estim[c(3:4), 1] = c(as.numeric(fit.SeY@param))

estim[c(5:6), 1] = c(as.numeric(fit.SpY@param))

estim.orig[c(1: 2), 1] = c(expit(estim[1, 1]), estim[2, 1])

estim.orig[c(3: 4), 1] = c(expit(estim[3, 1]), estim[4, 1])

estim.orig[c(5: 6), 1] = c(expit(estim[5, 1]), estim[6, 1])

rownames(estim)=c("PiY1", "PiY2", "SeY1", "SeY2", "SpY1", "SpY2")

rownames(estim.orig)=c("PiY1", "PiY2", "SeY1", "SeY2", "SpY1", "SpY2")

fit.PiY.var = matrix(c(fit.PiY@varparam[1,1], fit.PiY@varparam[1,2], fit.PiY@varparam[2,1],fit.PiY@varparam[2,2]), nrow=2,byrow=TRUE)

fit.SeY.var = matrix(c(fit.SeY@varparam[1,1], fit.SeY@varparam[1,2], fit.SeY@varparam[2,1],fit.SeY@varparam[2,2]), nrow=2,byrow=TRUE)

fit.SpY.var = matrix(c(fit.SpY@varparam[1,1], fit.SpY@varparam[1,2], fit.SpY@varparam[2,1],fit.SpY@varparam[2,2]), nrow=2,byrow=TRUE)

I00.hat = solve(m2*fit.PiY.var)

I11.hat = solve(m*fit.SeY.var)

I22.hat = solve(m*fit.SpY.var)

estim.PiY = estim[c(1:2), 1]

estim.SeY = estim[c(3:4), 1]

estim.SpY = estim[c(5:6), 1]

temp = Fisher.inf(mydat2\$PiY,mydat2\$SeY,mydat2\$SpY,mydat2\$n,mydat2\$n1,mydat2\$n0, mydat1\$SeY,mydat1\$SpY,mydat1\$n1,mydat1\$n0,estim.PiY, estim.SeY, estim.SpY)

I01.hat = temp\$I01 ; I02.hat = temp\$I02 ; I12.hat = temp\$I12

myoff.diag01 = solve(I00.hat)%*%I01.hat%*%solve(I11.hat)/m

myoff.diag02 = solve(I00.hat)%*%I02.hat%*%solve(I22.hat)/m

myoff.diag12 = solve(I11.hat)%*%I12.hat%*%solve(I22.hat)/m

myupper = cbind((fit.PiY.var), sqrt(m/m2)*(myoff.diag01), sqrt(m/m2)*(myoff.diag02))

mymiddle = cbind(sqrt(m/m2)*t(myoff.diag01), (fit.SeY.var), (myoff.diag12))

mylower = cbind(sqrt(m/m2)*t(myoff.diag02), t(myoff.diag12), (fit.SpY.var))

myV = rbind(myupper, mymiddle, mylower)

colnames(myV)=c("PiY1", "PiY2", "SeY1", "SeY2", "SpY1", "SpY2")

rownames(myV)=c("PiY1", "PiY2", "SeY1", "SeY2", "SpY1", "SpY2")

mu0.mbse = (estim.orig[1, 1]*(1-estim.orig[1, 1]))*sqrt(diag(myV)[1])

mu1.mbse = (estim.orig[3, 1]*(1-estim.orig[3, 1]))*sqrt(diag(myV)[3])

mu2.mbse = (estim.orig[5, 1]*(1-estim.orig[5, 1]))*sqrt(diag(myV)[5])

s1.2.mbse = sqrt(diag(myV)[2])

s2.2.mbse = sqrt(diag(myV)[4])

s3.2.mbse = sqrt(diag(myV)[6])

mbse.orig = c(mu0.mbse, s1.2.mbse, mu1.mbse, s2.2.mbse, mu2.mbse, s3.2.mbse)

return(list(coefficients=estim, estim.orig=estim.orig, vcov=myV))

}

tn.cl=function(n, PiY2, SeY1, SeY2, SpY1, SpY2, Nd, Nn){

estim = estim.orig = mbse.orig = matrix(NA, nrow=6, ncol=1)

m1=length(SeY1)

m2=length(SeY2)

m=m1+m2

PiY = c(rep(NA, length=m1),PiY2)

SeY = c(SeY1, SeY2)

SpY = c(SpY1, SpY2)

id=c(1:m)

mydat = data.frame(id, rep(n, length=m), PiY, Nd ,SeY,Nn, SpY)

names(mydat) = c("id", "Nt","PiY", "Nd","SeY","Nn","SpY") ## codebook for mydat: study id, number of total subjects, number of disease ppl among total subjects, number of diease ppl, number of being postive test among disease ppl, number of disease-free ppl, number of being negative test among disease-free ppl

mydat2=na.omit(mydat)

fit.PiY = glmmML(cbind(PiY, Nt-PiY)~1, data = mydat2, cluster = id, method= "ghq", n.points=8,     control=list(epsilon=1e-08 , maxit=50000,trace=FALSE))

fit.SeY = glmmML(cbind(SeY, Nd-SeY)~1, data = mydat, cluster = id, method= "ghq", n.points=8,     control=list(epsilon=1e-08 , maxit=50000,trace=FALSE))

fit.SpY = glmmML(cbind(SpY, Nn-SpY)~1, data = mydat, cluster = id, method= "ghq", n.points=8, control=list(epsilon=1e-08 , maxit=50000,trace=FALSE))

estim[c(1: 2),1] = c(as.numeric(fit.PiY\$coefficients), (fit.PiY\$sigma)^2)

estim[c(3: 4),1] = c(as.numeric(fit.SeY\$coefficients), (fit.SeY\$sigma)^2)

estim[c(5: 6),1] = c(as.numeric(fit.SpY\$coefficients), (fit.SpY\$sigma)^2)

estim.orig[c(1: 2),1] = c(plogis(estim[1,1]), estim[2,1])

estim.orig[c(3: 4),1] = c(plogis(estim[3,1]), estim[4,1])

estim.orig[c(5: 6),1] = c(plogis(estim[5,1]), estim[6,1])

rownames(estim)=c("PiY1", "PiY2", "SeY1", "SeY2", "SpY1", "SpY2")

rownames(estim.orig)=c("PiY1", "PiY2", "SeY1", "SeY2", "SpY1", "SpY2")

I00.hat = solve(m2*fit.PiY\$variance)

I11.hat = solve(m*fit.SeY\$variance)

I22.hat = solve(m*fit.SpY\$variance)

int1 = int2  = int3 = rep(NA, length=m2)

est0 = estim[c(1: 2),1] ##plug in the point estimate based on disease prevalence

for(i in 1: m2){

int1[i] = integrate(integrand1, lower=0, upper=1, n=mydat2\$Nt[i], mypar2=est0, y=mydat2\$PiY[i])\$value

int2[i] = integrate(integrand2, lower=0, upper=1, n=mydat2\$Nt[i], mypar2=est0, y=mydat2\$PiY[i])\$value

int3[i] = integrate(integrand3, lower=0, upper=1, n=mydat2\$Nt[i], mypar2=est0, y=mydat2\$PiY[i])\$value

}

##

derivePiY.mu = int2/int1

derivePiY.tau2 = int3/int1

B.PiY = cbind(derivePiY.mu, derivePiY.tau2)

int1 = int2  = int3 = rep(NA, length=m)

est1 = estim[c(3: 4),1] ##plug in the point estimate based on sensitivity

for(i in 1: m){

int1[i] = integrate(integrand1, lower=0, upper=1, n=mydat\$Nd[i], mypar2=est1, y=mydat\$SeY[i])\$value

int2[i] = integrate(integrand2, lower=0, upper=1, n=mydat\$Nd[i], mypar2=est1, y=mydat\$SeY[i])\$value

int3[i] = integrate(integrand3, lower=0, upper=1, n=mydat\$Nd[i], mypar2=est1, y=mydat\$SeY[i])\$value

}

deriveSeY.mu = int2/int1

deriveSeY.tau2 = int3/int1

B.SeY = cbind(deriveSeY.mu, deriveSeY.tau2)

B.SeY.del = B.SeY[-c(1:m1), ] ## delete the first m1 studies

int1 = int2  = int3 = rep(NA, length=m)

est2 = estim[c(5: 6),1] ##plug in the point estimate based on specificity

for(i in 1: m){

int1[i] = integrate(integrand1, lower=0, upper=1, n=mydat\$Nn[i], mypar2=est2, y=mydat\$SpY[i])\$value

int2[i] = integrate(integrand2, lower=0, upper=1, n=mydat\$Nn[i], mypar2=est2, y=mydat\$SpY[i])\$value

int3[i] = integrate(integrand3, lower=0, upper=1, n=mydat\$Nn[i], mypar2=est2, y=mydat\$SpY[i])\$value

}

deriveSpY.mu = int2/int1

deriveSpY.tau2 = int3/int1

B.SpY = cbind(deriveSpY.mu, deriveSpY.tau2)

B.SpY.del = B.SpY[-c(1:m1), ] ## delete the first m1 studies

I01.hat = t(B.PiY)%*%B.SeY.del/m2

I02.hat = t(B.PiY)%*%B.SpY.del/m2

I12.hat = t(B.SeY)%*%B.SpY/m

myoff.diag01 = solve(I00.hat)%*%I01.hat%*%solve(I11.hat)/m

myoff.diag02 = solve(I00.hat)%*%I02.hat%*%solve(I22.hat)/m

myoff.diag12 = solve(I11.hat)%*%I12.hat%*%solve(I22.hat)/m

myupper = cbind(fit.PiY\$variance, sqrt(m/m2)*(myoff.diag01), sqrt(m/m2)*(myoff.diag02))

mymiddle = cbind(sqrt(m/m2)*t(myoff.diag01), fit.SeY\$variance, (myoff.diag12))

mylower = cbind(sqrt(m/m2)*t(myoff.diag02), t(myoff.diag12), fit.SpY\$variance)

myV = rbind(myupper, mymiddle, mylower)

colnames(myV)=c("PiY1", "PiY2", "SeY1", "SeY2", "SpY1", "SpY2")

rownames(myV)=c("PiY1", "PiY2", "SeY1", "SeY2", "SpY1", "SpY2")

return(list(coefficients=estim, estim.orig=estim.orig, vcov=myV))

}

if (missing(data)){data=NULL}

if (is.null(data)){

stop("The dataset must be specified.")

}

if (missing(k)){k=NULL}

if (is.null(k)){

stop("The number of outcomes must be specified.")

}

if (missing(type)){type=NULL}

if (is.null(type)){

stop('The type of outcome must be specified. Please input "continuous" or "binary". ')

}

if (missing(method)){method=NULL}

if (is.null(method)){

stop("A method must be specified. ")

}

if (missing(rhow)){rhow=NULL}

if ((is.null(rhow))&(method=="nn.reml")){

stop('The within-study correlations are required for the input method "nn.reml".')

}

if (k>2){

stop('The method for MMA with more than 2 outcomes are currently under development. ')

}

if (type=="continuous"){

y1=data\$y1

s1=data\$s1

y2=data\$y2

s2=data\$s2

if(method=="nn.cl"){

fit=nn.cl(y1, s1, y2, s2)}

if(method=="nn.reml"){

fit=nn.reml(y1, s1, y2, s2, rhow)}

if(method=="nn.mom"){

fit=nn.mom(y1, s1, y2, s2)}

if(method=="nn.rs"){

fit=nn.rs(y1, s1, y2, s2)}

if ((method%in%c("nn.reml", "nn.cl", "nn.mom", "nn.rs"))!=1){

stop('The input method is not available for continuous outcomes. Please choose from

"nn.reml", "nn.cl", "nn.mom", or "nn.rs". ')

}

}

if(type=="binary"){

if (method=="bb.cl"){

y1=data\$y1

n1=data\$n1

y2=data\$y2

n2=data\$n2

fit=bb.cl(y1,n1,y2,n2)

}

if (method=="bn.cl"){

y1=data\$y1

n1=data\$n1

y2=data\$y2

n2=data\$n2

fit=bn.cl(y1,n1,y2,n2)

}

if (method=="tb.cl"){

n=data\$n

PiY=data\$PiY

SeY=data\$SeY

n1=data\$n1

SpY=data\$SpY

n2=data\$n2

fit=tb.cl(n, PiY, SeY, n1, SpY, n2)

}

if (method=="tn.cl"){

n=data\$n

PiY2=data\$PiY2

SeY1=data\$SeY1

SeY2=data\$SeY2

SpY1=data\$SpY1

SpY2=data\$SpY2

Nd=data\$Nd

Nn=data\$Nn

fit=tn.cl(n, PiY2, SeY1, SeY2, SpY1, SpY2, Nd, Nn)

}

if ((method%in%c("bb.cl", "bn.cl", "tb.cl", "tn.cl"))!=1){

stop('The input method is not available for binary outcomes. Please choose from

"bb.cl", "bn.cl", "tb.cl", or "tn.cl". ')

}

}

res=list(type=type, k=k, method=method, coefficients=fit\$coefficients, vcov=fit\$vcov)

if(method=="bb.cl"){

res=list(type=type, k=k, method=method, logOR=fit\$logOR, OR_CI=fit\$OR_CI, se=fit\$se,hessian.log=fit\$hessian.log,conv=fit\$conv,mle=fit\$mle)

}

class(res) = c("mmeta")

return(res)

}

This site was designed with the
.com
website builder. Create your website today.