Nothing
comparebcSens <-function (marker1, marker2, D, atSpec, alpha, plots){
if ((length(marker1) != length(D)) | (length(marker2) != length(D))) {
stop("ERROR: The length of the 'marker' and 'D' inputs must be equal.")
} else if (min(D) != 0 | max(D) != 1) {
stop("ERROR: Controls must be assigned a value of 0; cases must be assigned a value of 1. Both controls and cases should be included in the dataset.")
} else if (alpha <= 0 | alpha >= 1) {
stop("ERROR: The level of significance, alpha, should be set between 0 and 1. A common choice is 0.05.")
} else if (sum(is.na(marker1)) > 0 | sum(is.na(marker2)) > 0 | sum(is.na(D)) > 0 | sum(is.na(atSpec)) > 0) {
stop("ERROR: Please remove all missing data before running this function.")
} else if ((!is.numeric(atSpec)) | (length(atSpec) != 1)) {
stop("ERROR: 'atSpec' must be a single numeric value.")
} else if ((sum(marker1 < 0) > 0) | (sum(marker2 < 2)) > 0) {
stop("ERROR: To use the Box-Cox transformation, all marker values must be positive.")
} else {
erf <- function (x) 2 * pnorm(x * sqrt(2)) - 1
erfc <- function (x) 2 * pnorm(x * sqrt(2), lower.tail = FALSE)
erfinv <- function (x) qnorm((1 + x)/2)/sqrt(2)
erfcinv <- function (x) qnorm(x/2, lower.tail = FALSE)/sqrt(2)
pmalpha=qnorm(1-alpha/2);
if (plots!="on"){plots="off"}
W1a=marker1[D==0]
W1b=marker1[D==1]
W2a=marker2[D==0]
W2b=marker2[D==1]
Wa=cbind(W1a,W2a);
Wb=cbind(W1b,W2b);
na=length(W1a)
nb=length(W1b)
likbox2D<-function(W1a,W1b,W2a,W2b,lam){
na=length(W1a);
nb=length(W1b);
out=c();
#for (i in 1:length(h)){
W1alam= (W1a^lam[1]-1)/lam[1];
W2alam= (W2a^lam[2]-1)/lam[2];
W1blam= (W1b^lam[1]-1)/lam[1];
W2blam= (W2b^lam[2]-1)/lam[2];
Walam=cbind(W1alam,W2alam)
Wblam=cbind(W1blam,W2blam)
cova= cov(Walam)*(na-1)/na;
covb= cov(Wblam)*(nb-1)/nb;
mualam=c(mean(W1alam), mean(W2alam))
mublam=c(mean(W1blam), mean(W2blam))
out= -sum(log(dmvnorm(Walam,mualam,cova)))-sum(log(dmvnorm(Wblam,mublam,covb))) -(lam[1]-1)*sum(log(W1a))-(lam[2]-1)*sum(log(W2a))-(lam[1]-1)*sum(log(W1b))-(lam[2]-1)*sum(log(W2b));
return(out)
}
lam=c(2,2)
likbox2D(W1a,W1b,W2a,W2b,lam)
logL<-function(h){
likbox2D(W1a,W1b,W2a,W2b,h)
}
# init=c(0.8,0.8)
# lam=fminsearch(logL,init)
# lam=c(lam$optbase$xopt)
lam<- optim(c(1,1),logL,gr=NULL,method="BFGS", control=list(maxit=10000))
lam=c(lam$par)
lam
out <- optim(c(1,1), logL, method = "Nelder-Mead")
lam = out$par
###################################
W1alam= (W1a^lam[1]-1)/lam[1];
W2alam= (W2a^lam[2]-1)/lam[2];
W1blam= (W1b^lam[1]-1)/lam[1];
W2blam= (W2b^lam[2]-1)/lam[2];
Walam=cbind(W1alam,W2alam)
Wblam=cbind(W1blam,W2blam)
m1ahat=mean(W1alam)
m1bhat=mean(W1blam)
m2ahat=mean(W2alam)
m2bhat=mean(W2blam)
s1ahat=sqrt( var(W1alam)*(na-1)/na )
s1bhat=sqrt( var(W1blam)*(nb-1)/nb )
s2ahat=sqrt( var(W2alam)*(na-1)/na )
s2bhat=sqrt( var(W2blam)*(nb-1)/nb )
cova= cov(Walam)*(na-1)/na;cova=cova[1,2];
covb= cov(Wblam)*(nb-1)/nb;covb=covb[1,2];
lam1=lam[1]
lam2=lam[2]
##################################
h1_1=sum(-2/(2*s1ahat^2 - (2*cova^2)/s2ahat^2))
h1_2=sum(-(2*s1ahat*s2ahat^2*(cova*lam1 - lam2*s2ahat^2 - W2a^lam2*cova*lam1 + W1a^lam1*lam2*s2ahat^2 + cova*lam1*lam2*m2ahat - lam1*lam2*m1ahat*s2ahat^2))/(lam1*lam2*(s1ahat^2*s2ahat^2 - cova^2)^2))
h1_3=0
h1_4=0;
h1_5=sum(cova/(- cova^2 + s1ahat^2*s2ahat^2));
h1_6=sum((s2ahat*(2*cova^2*(lam2 - W1a^lam1*lam2 + lam1*lam2*m1ahat) - 2*cova*s1ahat^2*(lam1 - W2a^lam2*lam1 + lam1*lam2*m2ahat)))/(lam1*lam2*(s1ahat^2*s2ahat^2 - cova^2)^2))
h1_7=0
h1_8=0
h1_9=sum((cova^2*(lam2*m2ahat - W2a^lam2 + 1) + s2ahat^2*((lam2*m2ahat - W2a^lam2 + 1)*s1ahat^2 - 2*cova*lam2*m1ahat))/(lam2*(s1ahat^2*s2ahat^2 - cova^2)^2) - (cova*s2ahat^2*(2*lam2 - 2*W1a^lam1*lam2))/(lam1*lam2*(s1ahat^2*s2ahat^2 - cova^2)^2))
h1_10=0
h1_11=sum((2*s2ahat^2*(W1a^lam1*lam1*log(W1a) - W1a^lam1 + 1))/(lam1^2*(- 2*cova^2 + 2*s1ahat^2*s2ahat^2)))
h1_12=sum(-(cova*(W2a^lam2*lam2*log(W2a) - W2a^lam2 + 1))/(lam2^2*(- cova^2 + s1ahat^2*s2ahat^2)))
h2_1=sum(-(2*s1ahat*s2ahat^2*(cova*lam1 - lam2*s2ahat^2 - W2a^lam2*cova*lam1 + W1a^lam1*lam2*s2ahat^2 + cova*lam1*lam2*m2ahat - lam1*lam2*m1ahat*s2ahat^2))/(lam1*lam2*(s1ahat^2*s2ahat^2 - cova^2)^2))
h2_2=sum(- (3*s1ahat^2*s2ahat^6 + cova^4*(lam1^2*m2ahat^2 + lam1^2*s2ahat^2) - cova^3*(2*lam1*m2ahat*s2ahat^2 - 2*W1a^lam1*lam1*m2ahat*s2ahat^2 + 2*lam1^2*m1ahat*m2ahat*s2ahat^2) - cova*(6*lam1*m2ahat*s1ahat^2*s2ahat^4 - 6*W1a^lam1*lam1*m2ahat*s1ahat^2*s2ahat^4 + 6*lam1^2*m1ahat*m2ahat*s1ahat^2*s2ahat^4) + cova^2*(W1a^(2*lam1)*s2ahat^4 - 2*W1a^lam1*s2ahat^4 + s2ahat^4 + 2*lam1*m1ahat*s2ahat^4 + lam1^2*m1ahat^2*s2ahat^4 - 2*W1a^lam1*lam1*m1ahat*s2ahat^4 + 3*lam1^2*m2ahat^2*s1ahat^2*s2ahat^2) - 6*W1a^lam1*s1ahat^2*s2ahat^6 + 3*W1a^(2*lam1)*s1ahat^2*s2ahat^6 - lam1^2*s1ahat^4*s2ahat^6 + 3*lam1^2*m1ahat^2*s1ahat^2*s2ahat^6 + 6*lam1*m1ahat*s1ahat^2*s2ahat^6 - 6*W1a^lam1*lam1*m1ahat*s1ahat^2*s2ahat^6)/(lam1^2*(s1ahat^2*s2ahat^2 - cova^2)^3) - (cova^4*(W2a^(2*lam2)*lam1^2 - 2*W2a^lam2*lam1^2 + lam1^2) + cova^2*(3*lam1^2*s1ahat^2*s2ahat^2 + 3*W2a^(2*lam2)*lam1^2*s1ahat^2*s2ahat^2 - 6*W2a^lam2*lam1^2*s1ahat^2*s2ahat^2) - lam2*((2*W2a^lam2*lam1^2*m2ahat - 2*lam1^2*m2ahat)*cova^4 + (2*lam1*s2ahat^2 + 2*lam1^2*m1ahat*s2ahat^2 - 2*W1a^lam1*lam1*s2ahat^2 - 2*W2a^lam2*lam1*s2ahat^2 - 2*W2a^lam2*lam1^2*m1ahat*s2ahat^2 + 2*W1a^lam1*W2a^lam2*lam1*s2ahat^2)*cova^3 + (6*W2a^lam2*lam1^2*m2ahat*s1ahat^2*s2ahat^2 - 6*lam1^2*m2ahat*s1ahat^2*s2ahat^2)*cova^2 + (6*lam1*s1ahat^2*s2ahat^4 + 6*lam1^2*m1ahat*s1ahat^2*s2ahat^4 - 6*W1a^lam1*lam1*s1ahat^2*s2ahat^4 - 6*W2a^lam2*lam1*s1ahat^2*s2ahat^4 - 6*W2a^lam2*lam1^2*m1ahat*s1ahat^2*s2ahat^4 + 6*W1a^lam1*W2a^lam2*lam1*s1ahat^2*s2ahat^4)*cova))/(lam1^2*lam2^2*(s1ahat^2*s2ahat^2 - cova^2)^3))
h2_3=sum(0)
h2_4=sum(0)
h2_5=sum((s1ahat*(2*cova^2*(lam1 - W2a^lam2*lam1 + lam1*lam2*m2ahat) - 2*cova*s2ahat^2*(lam2 - W1a^lam1*lam2 + lam1*lam2*m1ahat)))/(lam1*lam2*(s1ahat^2*s2ahat^2 - cova^2)^2))
h2_6=sum((lam2*((2*cova*s2ahat^3*(2*lam1 + 2*lam1^2*m1ahat - 2*W1a^lam1*lam1 - 2*W2a^lam2*lam1 + 2*W1a^lam1*W2a^lam2*lam1 - 2*W2a^lam2*lam1^2*m1ahat) - 2*cova*s2ahat*(4*cova*lam1^2*m2ahat - 4*W2a^lam2*cova*lam1^2*m2ahat))*s1ahat^3 + 2*cova*s2ahat*(2*cova^2*lam1 + 2*cova^2*lam1^2*m1ahat - 2*W1a^lam1*cova^2*lam1 - 2*W2a^lam2*cova^2*lam1 - 2*W2a^lam2*cova^2*lam1^2*m1ahat + 2*W1a^lam1*W2a^lam2*cova^2*lam1)*s1ahat) - 2*cova*s1ahat^3*s2ahat*(2*cova*lam1^2 + 2*W2a^(2*lam2)*cova*lam1^2 - 4*W2a^lam2*cova*lam1^2))/(lam1^2*lam2^2*(s1ahat^2*s2ahat^2 - cova^2)^3) - ((4*cova^2*lam1^2*m2ahat^2*s2ahat - 2*cova*s2ahat^3*(2*lam1*m2ahat + cova*lam1^2 + 2*lam1^2*m1ahat*m2ahat - 2*W1a^lam1*lam1*m2ahat))*s1ahat^3 + (2*cova*(2*cova + 2*W1a^(2*lam1)*cova - 4*W1a^lam1*cova + 2*cova*lam1^2*m1ahat^2 + 4*cova*lam1*m1ahat - 4*W1a^lam1*cova*lam1*m1ahat)*s2ahat^3 + 2*cova*(cova^3*lam1^2 - 2*cova^2*lam1*m2ahat + 2*W1a^lam1*cova^2*lam1*m2ahat - 2*cova^2*lam1^2*m1ahat*m2ahat)*s2ahat)*s1ahat)/(lam1^2*(s1ahat^2*s2ahat^2 - cova^2)^3))
h2_7=sum(0)
h2_8=sum(0)
h2_9=sum(- (cova*(s2ahat^4*(- 4*lam2^2*m1ahat^2*s1ahat + 2*lam2^2*s1ahat^3) - s1ahat^3*s2ahat^2*(4*lam2*m2ahat - 4*W2a^lam2 + 2*W2a^(2*lam2) + 2*lam2^2*m2ahat^2 - 4*W2a^lam2*lam2*m2ahat + 2)) - cova^3*(s1ahat*(4*lam2*m2ahat - 4*W2a^lam2 + 2*W2a^(2*lam2) + 2*lam2^2*m2ahat^2 - 4*W2a^lam2*lam2*m2ahat + 2) + 2*lam2^2*s1ahat*s2ahat^2) + s1ahat^3*s2ahat^4*(2*lam2*m1ahat + 2*lam2^2*m1ahat*m2ahat - 2*W2a^lam2*lam2*m1ahat) + cova^2*s1ahat*s2ahat^2*(6*lam2*m1ahat + 6*lam2^2*m1ahat*m2ahat - 6*W2a^lam2*lam2*m1ahat))/(lam2^2*(s1ahat^2*s2ahat^2 - cova^2)^3) - (lam1*((6*lam2 + 6*lam2^2*m2ahat - 6*W1a^lam1*lam2 - 6*W2a^lam2*lam2 + 6*W1a^lam1*W2a^lam2*lam2 - 6*W1a^lam1*lam2^2*m2ahat)*cova^2*s1ahat*s2ahat^2 + (8*W1a^lam1*lam2^2*m1ahat - 8*lam2^2*m1ahat)*cova*s1ahat*s2ahat^4 + (2*lam2 + 2*lam2^2*m2ahat - 2*W1a^lam1*lam2 - 2*W2a^lam2*lam2 + 2*W1a^lam1*W2a^lam2*lam2 - 2*W1a^lam1*lam2^2*m2ahat)*s1ahat^3*s2ahat^4) - cova*s1ahat*s2ahat^4*(4*W1a^(2*lam1)*lam2^2 - 8*W1a^lam1*lam2^2 + 4*lam2^2))/(lam1^2*lam2^2*(s1ahat^2*s2ahat^2 - cova^2)^3))
h2_10=sum(0)
h2_11=sum((2*s1ahat*s2ahat^2*(W1a^lam1*lam1*log(W1a) - W1a^lam1 + 1)*(cova*lam1 - lam2*s2ahat^2 - W2a^lam2*cova*lam1 + W1a^lam1*lam2*s2ahat^2 + cova*lam1*lam2*m2ahat - lam1*lam2*m1ahat*s2ahat^2))/(lam1^3*lam2*(s1ahat^2*s2ahat^2 - cova^2)^2))
h2_12=sum(-(2*cova*s1ahat*(W2a^lam2*lam2*log(W2a) - W2a^lam2 + 1)*(cova*lam1 - lam2*s2ahat^2 - W2a^lam2*cova*lam1 + W1a^lam1*lam2*s2ahat^2 + cova*lam1*lam2*m2ahat - lam1*lam2*m1ahat*s2ahat^2))/(lam1*lam2^3*(s1ahat^2*s2ahat^2 - cova^2)^2))
h3_1=sum(0)
h3_2=sum(0)
h3_3=sum(-2/(2*s1bhat^2 - (2*covb^2)/s2bhat^2))
h3_4=sum(-(2*s1bhat*s2bhat^2*(covb*lam1 - lam2*s2bhat^2 - W2b^lam2*covb*lam1 + W1b^lam1*lam2*s2bhat^2 + covb*lam1*lam2*m2bhat - lam1*lam2*m1bhat*s2bhat^2))/(lam1*lam2*(s1bhat^2*s2bhat^2 - covb^2)^2))
h3_5=sum(0)
h3_6=sum(0)
h3_7=sum(covb/(- covb^2 + s1bhat^2*s2bhat^2))
h3_8=sum((s2bhat*(2*covb^2*(lam2 - W1b^lam1*lam2 + lam1*lam2*m1bhat) - 2*covb*s1bhat^2*(lam1 - W2b^lam2*lam1 + lam1*lam2*m2bhat)))/(lam1*lam2*(s1bhat^2*s2bhat^2 - covb^2)^2))
h3_9=sum(0)
h3_10=sum((covb^2*(lam2*m2bhat - W2b^lam2 + 1) + s2bhat^2*((lam2*m2bhat - W2b^lam2 + 1)*s1bhat^2 - 2*covb*lam2*m1bhat))/(lam2*(s1bhat^2*s2bhat^2 - covb^2)^2) - (covb*s2bhat^2*(2*lam2 - 2*W1b^lam1*lam2))/(lam1*lam2*(s1bhat^2*s2bhat^2 - covb^2)^2))
h3_11=sum((2*s2bhat^2*(W1b^lam1*lam1*log(W1b) - W1b^lam1 + 1))/(lam1^2*(- 2*covb^2 + 2*s1bhat^2*s2bhat^2)))
h3_12=sum(-(covb*(W2b^lam2*lam2*log(W2b) - W2b^lam2 + 1))/(lam2^2*(- covb^2 + s1bhat^2*s2bhat^2)))
h4_1=sum(0)
h4_2=sum(0)
h4_3=sum(-(2*s1bhat*s2bhat^2*(covb*lam1 - lam2*s2bhat^2 - W2b^lam2*covb*lam1 + W1b^lam1*lam2*s2bhat^2 + covb*lam1*lam2*m2bhat - lam1*lam2*m1bhat*s2bhat^2))/(lam1*lam2*(s1bhat^2*s2bhat^2 - covb^2)^2))
h4_4=sum(- (3*s1bhat^2*s2bhat^6 + covb^4*(lam1^2*m2bhat^2 + lam1^2*s2bhat^2) - covb^3*(2*lam1*m2bhat*s2bhat^2 - 2*W1b^lam1*lam1*m2bhat*s2bhat^2 + 2*lam1^2*m1bhat*m2bhat*s2bhat^2) - covb*(6*lam1*m2bhat*s1bhat^2*s2bhat^4 - 6*W1b^lam1*lam1*m2bhat*s1bhat^2*s2bhat^4 + 6*lam1^2*m1bhat*m2bhat*s1bhat^2*s2bhat^4) + covb^2*(W1b^(2*lam1)*s2bhat^4 - 2*W1b^lam1*s2bhat^4 + s2bhat^4 + 2*lam1*m1bhat*s2bhat^4 + lam1^2*m1bhat^2*s2bhat^4 - 2*W1b^lam1*lam1*m1bhat*s2bhat^4 + 3*lam1^2*m2bhat^2*s1bhat^2*s2bhat^2) - 6*W1b^lam1*s1bhat^2*s2bhat^6 + 3*W1b^(2*lam1)*s1bhat^2*s2bhat^6 - lam1^2*s1bhat^4*s2bhat^6 + 3*lam1^2*m1bhat^2*s1bhat^2*s2bhat^6 + 6*lam1*m1bhat*s1bhat^2*s2bhat^6 - 6*W1b^lam1*lam1*m1bhat*s1bhat^2*s2bhat^6)/(lam1^2*(s1bhat^2*s2bhat^2 - covb^2)^3) - (covb^4*(W2b^(2*lam2)*lam1^2 - 2*W2b^lam2*lam1^2 + lam1^2) + covb^2*(3*lam1^2*s1bhat^2*s2bhat^2 + 3*W2b^(2*lam2)*lam1^2*s1bhat^2*s2bhat^2 - 6*W2b^lam2*lam1^2*s1bhat^2*s2bhat^2) - lam2*((2*W2b^lam2*lam1^2*m2bhat - 2*lam1^2*m2bhat)*covb^4 + (2*lam1*s2bhat^2 + 2*lam1^2*m1bhat*s2bhat^2 - 2*W1b^lam1*lam1*s2bhat^2 - 2*W2b^lam2*lam1*s2bhat^2 - 2*W2b^lam2*lam1^2*m1bhat*s2bhat^2 + 2*W1b^lam1*W2b^lam2*lam1*s2bhat^2)*covb^3 + (6*W2b^lam2*lam1^2*m2bhat*s1bhat^2*s2bhat^2 - 6*lam1^2*m2bhat*s1bhat^2*s2bhat^2)*covb^2 + (6*lam1*s1bhat^2*s2bhat^4 + 6*lam1^2*m1bhat*s1bhat^2*s2bhat^4 - 6*W1b^lam1*lam1*s1bhat^2*s2bhat^4 - 6*W2b^lam2*lam1*s1bhat^2*s2bhat^4 - 6*W2b^lam2*lam1^2*m1bhat*s1bhat^2*s2bhat^4 + 6*W1b^lam1*W2b^lam2*lam1*s1bhat^2*s2bhat^4)*covb))/(lam1^2*lam2^2*(s1bhat^2*s2bhat^2 - covb^2)^3))
h4_5=sum(0)
h4_6=sum(0)
h4_7=sum((s1bhat*(2*covb^2*(lam1 - W2b^lam2*lam1 + lam1*lam2*m2bhat) - 2*covb*s2bhat^2*(lam2 - W1b^lam1*lam2 + lam1*lam2*m1bhat)))/(lam1*lam2*(s1bhat^2*s2bhat^2 - covb^2)^2))
h4_8=sum((lam2*((2*covb*s2bhat^3*(2*lam1 + 2*lam1^2*m1bhat - 2*W1b^lam1*lam1 - 2*W2b^lam2*lam1 + 2*W1b^lam1*W2b^lam2*lam1 - 2*W2b^lam2*lam1^2*m1bhat) - 2*covb*s2bhat*(4*covb*lam1^2*m2bhat - 4*W2b^lam2*covb*lam1^2*m2bhat))*s1bhat^3 + 2*covb*s2bhat*(2*covb^2*lam1 + 2*covb^2*lam1^2*m1bhat - 2*W1b^lam1*covb^2*lam1 - 2*W2b^lam2*covb^2*lam1 - 2*W2b^lam2*covb^2*lam1^2*m1bhat + 2*W1b^lam1*W2b^lam2*covb^2*lam1)*s1bhat) - 2*covb*s1bhat^3*s2bhat*(2*covb*lam1^2 + 2*W2b^(2*lam2)*covb*lam1^2 - 4*W2b^lam2*covb*lam1^2))/(lam1^2*lam2^2*(s1bhat^2*s2bhat^2 - covb^2)^3) - ((4*covb^2*lam1^2*m2bhat^2*s2bhat - 2*covb*s2bhat^3*(2*lam1*m2bhat + covb*lam1^2 + 2*lam1^2*m1bhat*m2bhat - 2*W1b^lam1*lam1*m2bhat))*s1bhat^3 + (2*covb*(2*covb + 2*W1b^(2*lam1)*covb - 4*W1b^lam1*covb + 2*covb*lam1^2*m1bhat^2 + 4*covb*lam1*m1bhat - 4*W1b^lam1*covb*lam1*m1bhat)*s2bhat^3 + 2*covb*(covb^3*lam1^2 - 2*covb^2*lam1*m2bhat + 2*W1b^lam1*covb^2*lam1*m2bhat - 2*covb^2*lam1^2*m1bhat*m2bhat)*s2bhat)*s1bhat)/(lam1^2*(s1bhat^2*s2bhat^2 - covb^2)^3))
h4_9=sum(0)
h4_10=sum(- (covb*(s2bhat^4*(- 4*lam2^2*m1bhat^2*s1bhat + 2*lam2^2*s1bhat^3) - s1bhat^3*s2bhat^2*(4*lam2*m2bhat - 4*W2b^lam2 + 2*W2b^(2*lam2) + 2*lam2^2*m2bhat^2 - 4*W2b^lam2*lam2*m2bhat + 2)) - covb^3*(s1bhat*(4*lam2*m2bhat - 4*W2b^lam2 + 2*W2b^(2*lam2) + 2*lam2^2*m2bhat^2 - 4*W2b^lam2*lam2*m2bhat + 2) + 2*lam2^2*s1bhat*s2bhat^2) + s1bhat^3*s2bhat^4*(2*lam2*m1bhat + 2*lam2^2*m1bhat*m2bhat - 2*W2b^lam2*lam2*m1bhat) + covb^2*s1bhat*s2bhat^2*(6*lam2*m1bhat + 6*lam2^2*m1bhat*m2bhat - 6*W2b^lam2*lam2*m1bhat))/(lam2^2*(s1bhat^2*s2bhat^2 - covb^2)^3) - (lam1*((6*lam2 + 6*lam2^2*m2bhat - 6*W1b^lam1*lam2 - 6*W2b^lam2*lam2 + 6*W1b^lam1*W2b^lam2*lam2 - 6*W1b^lam1*lam2^2*m2bhat)*covb^2*s1bhat*s2bhat^2 + (8*W1b^lam1*lam2^2*m1bhat - 8*lam2^2*m1bhat)*covb*s1bhat*s2bhat^4 + (2*lam2 + 2*lam2^2*m2bhat - 2*W1b^lam1*lam2 - 2*W2b^lam2*lam2 + 2*W1b^lam1*W2b^lam2*lam2 - 2*W1b^lam1*lam2^2*m2bhat)*s1bhat^3*s2bhat^4) - covb*s1bhat*s2bhat^4*(4*W1b^(2*lam1)*lam2^2 - 8*W1b^lam1*lam2^2 + 4*lam2^2))/(lam1^2*lam2^2*(s1bhat^2*s2bhat^2 - covb^2)^3))
h4_11=sum((2*s1bhat*s2bhat^2*(W1b^lam1*lam1*log(W1b) - W1b^lam1 + 1)*(covb*lam1 - lam2*s2bhat^2 - W2b^lam2*covb*lam1 + W1b^lam1*lam2*s2bhat^2 + covb*lam1*lam2*m2bhat - lam1*lam2*m1bhat*s2bhat^2))/(lam1^3*lam2*(s1bhat^2*s2bhat^2 - covb^2)^2))
h4_12=sum(-(2*covb*s1bhat*(W2b^lam2*lam2*log(W2b) - W2b^lam2 + 1)*(covb*lam1 - lam2*s2bhat^2 - W2b^lam2*covb*lam1 + W1b^lam1*lam2*s2bhat^2 + covb*lam1*lam2*m2bhat - lam1*lam2*m1bhat*s2bhat^2))/(lam1*lam2^3*(s1bhat^2*s2bhat^2 - covb^2)^2))
h5_1=sum(cova/(- cova^2 + s1ahat^2*s2ahat^2))
h5_2=sum((s1ahat*(2*cova^2*(lam1 - W2a^lam2*lam1 + lam1*lam2*m2ahat) - 2*cova*s2ahat^2*(lam2 - W1a^lam1*lam2 + lam1*lam2*m1ahat)))/(lam1*lam2*(s1ahat^2*s2ahat^2 - cova^2)^2))
h5_3=0
h5_4=0
h5_5=sum(-2/(2*s2ahat^2 - (2*cova^2)/s1ahat^2))
h5_6=sum(-(2*s1ahat^2*s2ahat*(cova*lam2 - lam1*s1ahat^2 - W1a^lam1*cova*lam2 + W2a^lam2*lam1*s1ahat^2 + cova*lam1*lam2*m1ahat - lam1*lam2*m2ahat*s1ahat^2))/(lam1*lam2*(s1ahat^2*s2ahat^2 - cova^2)^2))
h5_7=0
h5_8=0
h5_9=sum((cova^2*(lam1*m1ahat - W1a^lam1 + 1) + s1ahat^2*((lam1*m1ahat - W1a^lam1 + 1)*s2ahat^2 - 2*cova*lam1*m2ahat))/(lam1*(s1ahat^2*s2ahat^2 - cova^2)^2) - (cova*s1ahat^2*(2*lam1 - 2*W2a^lam2*lam1))/(lam1*lam2*(s1ahat^2*s2ahat^2 - cova^2)^2))
h5_10=0
h5_11=sum(-(cova*(W1a^lam1*lam1*log(W1a) - W1a^lam1 + 1))/(lam1^2*(- cova^2 + s1ahat^2*s2ahat^2)))
h5_12=sum((2*s1ahat^2*(W2a^lam2*lam2*log(W2a) - W2a^lam2 + 1))/(lam2^2*(- 2*cova^2 + 2*s1ahat^2*s2ahat^2)))
h6_1=sum((s2ahat*(2*cova^2*(lam2 - W1a^lam1*lam2 + lam1*lam2*m1ahat) - 2*cova*s1ahat^2*(lam1 - W2a^lam2*lam1 + lam1*lam2*m2ahat)))/(lam1*lam2*(s1ahat^2*s2ahat^2 - cova^2)^2))
h6_2=sum((lam2*((2*cova*s2ahat^3*(2*lam1 + 2*lam1^2*m1ahat - 2*W1a^lam1*lam1 - 2*W2a^lam2*lam1 + 2*W1a^lam1*W2a^lam2*lam1 - 2*W2a^lam2*lam1^2*m1ahat) - 2*cova*s2ahat*(4*cova*lam1^2*m2ahat - 4*W2a^lam2*cova*lam1^2*m2ahat))*s1ahat^3 + 2*cova*s2ahat*(2*cova^2*lam1 + 2*cova^2*lam1^2*m1ahat - 2*W1a^lam1*cova^2*lam1 - 2*W2a^lam2*cova^2*lam1 - 2*W2a^lam2*cova^2*lam1^2*m1ahat + 2*W1a^lam1*W2a^lam2*cova^2*lam1)*s1ahat) - 2*cova*s1ahat^3*s2ahat*(2*cova*lam1^2 + 2*W2a^(2*lam2)*cova*lam1^2 - 4*W2a^lam2*cova*lam1^2))/(lam1^2*lam2^2*(s1ahat^2*s2ahat^2 - cova^2)^3) - ((4*cova^2*lam1^2*m2ahat^2*s2ahat - 2*cova*s2ahat^3*(2*lam1*m2ahat + cova*lam1^2 + 2*lam1^2*m1ahat*m2ahat - 2*W1a^lam1*lam1*m2ahat))*s1ahat^3 + (2*cova*(2*cova + 2*W1a^(2*lam1)*cova - 4*W1a^lam1*cova + 2*cova*lam1^2*m1ahat^2 + 4*cova*lam1*m1ahat - 4*W1a^lam1*cova*lam1*m1ahat)*s2ahat^3 + 2*cova*(cova^3*lam1^2 - 2*cova^2*lam1*m2ahat + 2*W1a^lam1*cova^2*lam1*m2ahat - 2*cova^2*lam1^2*m1ahat*m2ahat)*s2ahat)*s1ahat)/(lam1^2*(s1ahat^2*s2ahat^2 - cova^2)^3))
h6_3=0
h6_4=0
h6_5=sum(-(2*s1ahat^2*s2ahat*(cova*lam2 - lam1*s1ahat^2 - W1a^lam1*cova*lam2 + W2a^lam2*lam1*s1ahat^2 + cova*lam1*lam2*m1ahat - lam1*lam2*m2ahat*s1ahat^2))/(lam1*lam2*(s1ahat^2*s2ahat^2 - cova^2)^2))
h6_6=sum(- (3*s1ahat^6*s2ahat^2 + cova^4*(lam2^2*m1ahat^2 + lam2^2*s1ahat^2) - cova^3*(2*lam2*m1ahat*s1ahat^2 - 2*W2a^lam2*lam2*m1ahat*s1ahat^2 + 2*lam2^2*m1ahat*m2ahat*s1ahat^2) - cova*(6*lam2*m1ahat*s1ahat^4*s2ahat^2 - 6*W2a^lam2*lam2*m1ahat*s1ahat^4*s2ahat^2 + 6*lam2^2*m1ahat*m2ahat*s1ahat^4*s2ahat^2) + cova^2*(W2a^(2*lam2)*s1ahat^4 - 2*W2a^lam2*s1ahat^4 + s1ahat^4 + 2*lam2*m2ahat*s1ahat^4 + lam2^2*m2ahat^2*s1ahat^4 - 2*W2a^lam2*lam2*m2ahat*s1ahat^4 + 3*lam2^2*m1ahat^2*s1ahat^2*s2ahat^2) - 6*W2a^lam2*s1ahat^6*s2ahat^2 + 3*W2a^(2*lam2)*s1ahat^6*s2ahat^2 - lam2^2*s1ahat^6*s2ahat^4 + 3*lam2^2*m2ahat^2*s1ahat^6*s2ahat^2 + 6*lam2*m2ahat*s1ahat^6*s2ahat^2 - 6*W2a^lam2*lam2*m2ahat*s1ahat^6*s2ahat^2)/(lam2^2*(s1ahat^2*s2ahat^2 - cova^2)^3) - (cova^4*(W1a^(2*lam1)*lam2^2 - 2*W1a^lam1*lam2^2 + lam2^2) + cova^2*(3*lam2^2*s1ahat^2*s2ahat^2 + 3*W1a^(2*lam1)*lam2^2*s1ahat^2*s2ahat^2 - 6*W1a^lam1*lam2^2*s1ahat^2*s2ahat^2) - lam1*((2*W1a^lam1*lam2^2*m1ahat - 2*lam2^2*m1ahat)*cova^4 + (2*lam2*s1ahat^2 + 2*lam2^2*m2ahat*s1ahat^2 - 2*W1a^lam1*lam2*s1ahat^2 - 2*W2a^lam2*lam2*s1ahat^2 - 2*W1a^lam1*lam2^2*m2ahat*s1ahat^2 + 2*W1a^lam1*W2a^lam2*lam2*s1ahat^2)*cova^3 + (6*W1a^lam1*lam2^2*m1ahat*s1ahat^2*s2ahat^2 - 6*lam2^2*m1ahat*s1ahat^2*s2ahat^2)*cova^2 + (6*lam2*s1ahat^4*s2ahat^2 + 6*lam2^2*m2ahat*s1ahat^4*s2ahat^2 - 6*W1a^lam1*lam2*s1ahat^4*s2ahat^2 - 6*W2a^lam2*lam2*s1ahat^4*s2ahat^2 - 6*W1a^lam1*lam2^2*m2ahat*s1ahat^4*s2ahat^2 + 6*W1a^lam1*W2a^lam2*lam2*s1ahat^4*s2ahat^2)*cova))/(lam1^2*lam2^2*(s1ahat^2*s2ahat^2 - cova^2)^3))
h6_7=0
h6_8=0
h6_9=sum(- (cova*(s1ahat^4*(- 4*lam1^2*m2ahat^2*s2ahat + 2*lam1^2*s2ahat^3) - s1ahat^2*s2ahat^3*(4*lam1*m1ahat - 4*W1a^lam1 + 2*W1a^(2*lam1) + 2*lam1^2*m1ahat^2 - 4*W1a^lam1*lam1*m1ahat + 2)) - cova^3*(s2ahat*(4*lam1*m1ahat - 4*W1a^lam1 + 2*W1a^(2*lam1) + 2*lam1^2*m1ahat^2 - 4*W1a^lam1*lam1*m1ahat + 2) + 2*lam1^2*s1ahat^2*s2ahat) + s1ahat^4*s2ahat^3*(2*lam1*m2ahat + 2*lam1^2*m1ahat*m2ahat - 2*W1a^lam1*lam1*m2ahat) + cova^2*s1ahat^2*s2ahat*(6*lam1*m2ahat + 6*lam1^2*m1ahat*m2ahat - 6*W1a^lam1*lam1*m2ahat))/(lam1^2*(s1ahat^2*s2ahat^2 - cova^2)^3) - (lam2*((6*lam1 + 6*lam1^2*m1ahat - 6*W1a^lam1*lam1 - 6*W2a^lam2*lam1 + 6*W1a^lam1*W2a^lam2*lam1 - 6*W2a^lam2*lam1^2*m1ahat)*cova^2*s1ahat^2*s2ahat + (8*W2a^lam2*lam1^2*m2ahat - 8*lam1^2*m2ahat)*cova*s1ahat^4*s2ahat + (2*lam1 + 2*lam1^2*m1ahat - 2*W1a^lam1*lam1 - 2*W2a^lam2*lam1 + 2*W1a^lam1*W2a^lam2*lam1 - 2*W2a^lam2*lam1^2*m1ahat)*s1ahat^4*s2ahat^3) - cova*s1ahat^4*s2ahat*(4*W2a^(2*lam2)*lam1^2 - 8*W2a^lam2*lam1^2 + 4*lam1^2))/(lam1^2*lam2^2*(s1ahat^2*s2ahat^2 - cova^2)^3))
h6_10=0
h6_11=sum(-(2*cova*s2ahat*(W1a^lam1*lam1*log(W1a) - W1a^lam1 + 1)*(cova*lam2 - lam1*s1ahat^2 - W1a^lam1*cova*lam2 + W2a^lam2*lam1*s1ahat^2 + cova*lam1*lam2*m1ahat - lam1*lam2*m2ahat*s1ahat^2))/(lam1^3*lam2*(s1ahat^2*s2ahat^2 - cova^2)^2))
h6_12=sum((2*s1ahat^2*s2ahat*(W2a^lam2*lam2*log(W2a) - W2a^lam2 + 1)*(cova*lam2 - lam1*s1ahat^2 - W1a^lam1*cova*lam2 + W2a^lam2*lam1*s1ahat^2 + cova*lam1*lam2*m1ahat - lam1*lam2*m2ahat*s1ahat^2))/(lam1*lam2^3*(s1ahat^2*s2ahat^2 - cova^2)^2))
h7_1=0
h7_2=0
h7_3=sum(covb/(- covb^2 + s1bhat^2*s2bhat^2))
h7_4=sum((s1bhat*(2*covb^2*(lam1 - W2b^lam2*lam1 + lam1*lam2*m2bhat) - 2*covb*s2bhat^2*(lam2 - W1b^lam1*lam2 + lam1*lam2*m1bhat)))/(lam1*lam2*(s1bhat^2*s2bhat^2 - covb^2)^2))
h7_5=0
h7_6=0
h7_7=sum(-2/(2*s2bhat^2 - (2*covb^2)/s1bhat^2))
h7_8=sum(-(2*s1bhat^2*s2bhat*(covb*lam2 - lam1*s1bhat^2 - W1b^lam1*covb*lam2 + W2b^lam2*lam1*s1bhat^2 + covb*lam1*lam2*m1bhat - lam1*lam2*m2bhat*s1bhat^2))/(lam1*lam2*(s1bhat^2*s2bhat^2 - covb^2)^2))
h7_9=0
h7_10=sum((covb^2*(lam1*m1bhat - W1b^lam1 + 1) + s1bhat^2*((lam1*m1bhat - W1b^lam1 + 1)*s2bhat^2 - 2*covb*lam1*m2bhat))/(lam1*(s1bhat^2*s2bhat^2 - covb^2)^2) - (covb*s1bhat^2*(2*lam1 - 2*W2b^lam2*lam1))/(lam1*lam2*(s1bhat^2*s2bhat^2 - covb^2)^2))
h7_11=sum(-(covb*(W1b^lam1*lam1*log(W1b) - W1b^lam1 + 1))/(lam1^2*(- covb^2 + s1bhat^2*s2bhat^2)))
h7_12=sum((2*s1bhat^2*(W2b^lam2*lam2*log(W2b) - W2b^lam2 + 1))/(lam2^2*(- 2*covb^2 + 2*s1bhat^2*s2bhat^2)))
h8_1=0
h8_2=0
h8_3=sum((s2bhat*(2*covb^2*(lam2 - W1b^lam1*lam2 + lam1*lam2*m1bhat) - 2*covb*s1bhat^2*(lam1 - W2b^lam2*lam1 + lam1*lam2*m2bhat)))/(lam1*lam2*(s1bhat^2*s2bhat^2 - covb^2)^2))
h8_4=sum((lam2*((2*covb*s2bhat^3*(2*lam1 + 2*lam1^2*m1bhat - 2*W1b^lam1*lam1 - 2*W2b^lam2*lam1 + 2*W1b^lam1*W2b^lam2*lam1 - 2*W2b^lam2*lam1^2*m1bhat) - 2*covb*s2bhat*(4*covb*lam1^2*m2bhat - 4*W2b^lam2*covb*lam1^2*m2bhat))*s1bhat^3 + 2*covb*s2bhat*(2*covb^2*lam1 + 2*covb^2*lam1^2*m1bhat - 2*W1b^lam1*covb^2*lam1 - 2*W2b^lam2*covb^2*lam1 - 2*W2b^lam2*covb^2*lam1^2*m1bhat + 2*W1b^lam1*W2b^lam2*covb^2*lam1)*s1bhat) - 2*covb*s1bhat^3*s2bhat*(2*covb*lam1^2 + 2*W2b^(2*lam2)*covb*lam1^2 - 4*W2b^lam2*covb*lam1^2))/(lam1^2*lam2^2*(s1bhat^2*s2bhat^2 - covb^2)^3) - ((4*covb^2*lam1^2*m2bhat^2*s2bhat - 2*covb*s2bhat^3*(2*lam1*m2bhat + covb*lam1^2 + 2*lam1^2*m1bhat*m2bhat - 2*W1b^lam1*lam1*m2bhat))*s1bhat^3 + (2*covb*(2*covb + 2*W1b^(2*lam1)*covb - 4*W1b^lam1*covb + 2*covb*lam1^2*m1bhat^2 + 4*covb*lam1*m1bhat - 4*W1b^lam1*covb*lam1*m1bhat)*s2bhat^3 + 2*covb*(covb^3*lam1^2 - 2*covb^2*lam1*m2bhat + 2*W1b^lam1*covb^2*lam1*m2bhat - 2*covb^2*lam1^2*m1bhat*m2bhat)*s2bhat)*s1bhat)/(lam1^2*(s1bhat^2*s2bhat^2 - covb^2)^3))
h8_5=0
h8_6=0
h8_7=sum(-(2*s1bhat^2*s2bhat*(covb*lam2 - lam1*s1bhat^2 - W1b^lam1*covb*lam2 + W2b^lam2*lam1*s1bhat^2 + covb*lam1*lam2*m1bhat - lam1*lam2*m2bhat*s1bhat^2))/(lam1*lam2*(s1bhat^2*s2bhat^2 - covb^2)^2))
h8_8=sum(- (3*s1bhat^6*s2bhat^2 + covb^4*(lam2^2*m1bhat^2 + lam2^2*s1bhat^2) - covb^3*(2*lam2*m1bhat*s1bhat^2 - 2*W2b^lam2*lam2*m1bhat*s1bhat^2 + 2*lam2^2*m1bhat*m2bhat*s1bhat^2) - covb*(6*lam2*m1bhat*s1bhat^4*s2bhat^2 - 6*W2b^lam2*lam2*m1bhat*s1bhat^4*s2bhat^2 + 6*lam2^2*m1bhat*m2bhat*s1bhat^4*s2bhat^2) + covb^2*(W2b^(2*lam2)*s1bhat^4 - 2*W2b^lam2*s1bhat^4 + s1bhat^4 + 2*lam2*m2bhat*s1bhat^4 + lam2^2*m2bhat^2*s1bhat^4 - 2*W2b^lam2*lam2*m2bhat*s1bhat^4 + 3*lam2^2*m1bhat^2*s1bhat^2*s2bhat^2) - 6*W2b^lam2*s1bhat^6*s2bhat^2 + 3*W2b^(2*lam2)*s1bhat^6*s2bhat^2 - lam2^2*s1bhat^6*s2bhat^4 + 3*lam2^2*m2bhat^2*s1bhat^6*s2bhat^2 + 6*lam2*m2bhat*s1bhat^6*s2bhat^2 - 6*W2b^lam2*lam2*m2bhat*s1bhat^6*s2bhat^2)/(lam2^2*(s1bhat^2*s2bhat^2 - covb^2)^3) - (covb^4*(W1b^(2*lam1)*lam2^2 - 2*W1b^lam1*lam2^2 + lam2^2) + covb^2*(3*lam2^2*s1bhat^2*s2bhat^2 + 3*W1b^(2*lam1)*lam2^2*s1bhat^2*s2bhat^2 - 6*W1b^lam1*lam2^2*s1bhat^2*s2bhat^2) - lam1*((2*W1b^lam1*lam2^2*m1bhat - 2*lam2^2*m1bhat)*covb^4 + (2*lam2*s1bhat^2 + 2*lam2^2*m2bhat*s1bhat^2 - 2*W1b^lam1*lam2*s1bhat^2 - 2*W2b^lam2*lam2*s1bhat^2 - 2*W1b^lam1*lam2^2*m2bhat*s1bhat^2 + 2*W1b^lam1*W2b^lam2*lam2*s1bhat^2)*covb^3 + (6*W1b^lam1*lam2^2*m1bhat*s1bhat^2*s2bhat^2 - 6*lam2^2*m1bhat*s1bhat^2*s2bhat^2)*covb^2 + (6*lam2*s1bhat^4*s2bhat^2 + 6*lam2^2*m2bhat*s1bhat^4*s2bhat^2 - 6*W1b^lam1*lam2*s1bhat^4*s2bhat^2 - 6*W2b^lam2*lam2*s1bhat^4*s2bhat^2 - 6*W1b^lam1*lam2^2*m2bhat*s1bhat^4*s2bhat^2 + 6*W1b^lam1*W2b^lam2*lam2*s1bhat^4*s2bhat^2)*covb))/(lam1^2*lam2^2*(s1bhat^2*s2bhat^2 - covb^2)^3))
h8_9=0
h8_10=sum(- (covb*(s1bhat^4*(- 4*lam1^2*m2bhat^2*s2bhat + 2*lam1^2*s2bhat^3) - s1bhat^2*s2bhat^3*(4*lam1*m1bhat - 4*W1b^lam1 + 2*W1b^(2*lam1) + 2*lam1^2*m1bhat^2 - 4*W1b^lam1*lam1*m1bhat + 2)) - covb^3*(s2bhat*(4*lam1*m1bhat - 4*W1b^lam1 + 2*W1b^(2*lam1) + 2*lam1^2*m1bhat^2 - 4*W1b^lam1*lam1*m1bhat + 2) + 2*lam1^2*s1bhat^2*s2bhat) + s1bhat^4*s2bhat^3*(2*lam1*m2bhat + 2*lam1^2*m1bhat*m2bhat - 2*W1b^lam1*lam1*m2bhat) + covb^2*s1bhat^2*s2bhat*(6*lam1*m2bhat + 6*lam1^2*m1bhat*m2bhat - 6*W1b^lam1*lam1*m2bhat))/(lam1^2*(s1bhat^2*s2bhat^2 - covb^2)^3) - (lam2*((6*lam1 + 6*lam1^2*m1bhat - 6*W1b^lam1*lam1 - 6*W2b^lam2*lam1 + 6*W1b^lam1*W2b^lam2*lam1 - 6*W2b^lam2*lam1^2*m1bhat)*covb^2*s1bhat^2*s2bhat + (8*W2b^lam2*lam1^2*m2bhat - 8*lam1^2*m2bhat)*covb*s1bhat^4*s2bhat + (2*lam1 + 2*lam1^2*m1bhat - 2*W1b^lam1*lam1 - 2*W2b^lam2*lam1 + 2*W1b^lam1*W2b^lam2*lam1 - 2*W2b^lam2*lam1^2*m1bhat)*s1bhat^4*s2bhat^3) - covb*s1bhat^4*s2bhat*(4*W2b^(2*lam2)*lam1^2 - 8*W2b^lam2*lam1^2 + 4*lam1^2))/(lam1^2*lam2^2*(s1bhat^2*s2bhat^2 - covb^2)^3))
h8_11=sum(-(2*covb*s2bhat*(W1b^lam1*lam1*log(W1b) - W1b^lam1 + 1)*(covb*lam2 - lam1*s1bhat^2 - W1b^lam1*covb*lam2 + W2b^lam2*lam1*s1bhat^2 + covb*lam1*lam2*m1bhat - lam1*lam2*m2bhat*s1bhat^2))/(lam1^3*lam2*(s1bhat^2*s2bhat^2 - covb^2)^2))
h8_12=sum((2*s1bhat^2*s2bhat*(W2b^lam2*lam2*log(W2b) - W2b^lam2 + 1)*(covb*lam2 - lam1*s1bhat^2 - W1b^lam1*covb*lam2 + W2b^lam2*lam1*s1bhat^2 + covb*lam1*lam2*m1bhat - lam1*lam2*m2bhat*s1bhat^2))/(lam1*lam2^3*(s1bhat^2*s2bhat^2 - covb^2)^2))
h9_1=sum((cova^2*(lam2*m2ahat - W2a^lam2 + 1) + s2ahat^2*((lam2*m2ahat - W2a^lam2 + 1)*s1ahat^2 - 2*cova*lam2*m1ahat))/(lam2*(s1ahat^2*s2ahat^2 - cova^2)^2) - (cova*s2ahat^2*(2*lam2 - 2*W1a^lam1*lam2))/(lam1*lam2*(s1ahat^2*s2ahat^2 - cova^2)^2))
h9_2=sum(- (cova*(s2ahat^4*(- 4*lam2^2*m1ahat^2*s1ahat + 2*lam2^2*s1ahat^3) - s1ahat^3*s2ahat^2*(4*lam2*m2ahat - 4*W2a^lam2 + 2*W2a^(2*lam2) + 2*lam2^2*m2ahat^2 - 4*W2a^lam2*lam2*m2ahat + 2)) - cova^3*(s1ahat*(4*lam2*m2ahat - 4*W2a^lam2 + 2*W2a^(2*lam2) + 2*lam2^2*m2ahat^2 - 4*W2a^lam2*lam2*m2ahat + 2) + 2*lam2^2*s1ahat*s2ahat^2) + s1ahat^3*s2ahat^4*(2*lam2*m1ahat + 2*lam2^2*m1ahat*m2ahat - 2*W2a^lam2*lam2*m1ahat) + cova^2*s1ahat*s2ahat^2*(6*lam2*m1ahat + 6*lam2^2*m1ahat*m2ahat - 6*W2a^lam2*lam2*m1ahat))/(lam2^2*(s1ahat^2*s2ahat^2 - cova^2)^3) - (lam1*((6*lam2 + 6*lam2^2*m2ahat - 6*W1a^lam1*lam2 - 6*W2a^lam2*lam2 + 6*W1a^lam1*W2a^lam2*lam2 - 6*W1a^lam1*lam2^2*m2ahat)*cova^2*s1ahat*s2ahat^2 + (8*W1a^lam1*lam2^2*m1ahat - 8*lam2^2*m1ahat)*cova*s1ahat*s2ahat^4 + (2*lam2 + 2*lam2^2*m2ahat - 2*W1a^lam1*lam2 - 2*W2a^lam2*lam2 + 2*W1a^lam1*W2a^lam2*lam2 - 2*W1a^lam1*lam2^2*m2ahat)*s1ahat^3*s2ahat^4) - cova*s1ahat*s2ahat^4*(4*W1a^(2*lam1)*lam2^2 - 8*W1a^lam1*lam2^2 + 4*lam2^2))/(lam1^2*lam2^2*(s1ahat^2*s2ahat^2 - cova^2)^3))
h9_3=0
h9_4=0
h9_5=sum((cova^2*(lam1*m1ahat - W1a^lam1 + 1) + s1ahat^2*((lam1*m1ahat - W1a^lam1 + 1)*s2ahat^2 - 2*cova*lam1*m2ahat))/(lam1*(s1ahat^2*s2ahat^2 - cova^2)^2) - (cova*s1ahat^2*(2*lam1 - 2*W2a^lam2*lam1))/(lam1*lam2*(s1ahat^2*s2ahat^2 - cova^2)^2))
h9_6=sum(- (cova*(s1ahat^4*(- 4*lam1^2*m2ahat^2*s2ahat + 2*lam1^2*s2ahat^3) - s1ahat^2*s2ahat^3*(4*lam1*m1ahat - 4*W1a^lam1 + 2*W1a^(2*lam1) + 2*lam1^2*m1ahat^2 - 4*W1a^lam1*lam1*m1ahat + 2)) - cova^3*(s2ahat*(4*lam1*m1ahat - 4*W1a^lam1 + 2*W1a^(2*lam1) + 2*lam1^2*m1ahat^2 - 4*W1a^lam1*lam1*m1ahat + 2) + 2*lam1^2*s1ahat^2*s2ahat) + s1ahat^4*s2ahat^3*(2*lam1*m2ahat + 2*lam1^2*m1ahat*m2ahat - 2*W1a^lam1*lam1*m2ahat) + cova^2*s1ahat^2*s2ahat*(6*lam1*m2ahat + 6*lam1^2*m1ahat*m2ahat - 6*W1a^lam1*lam1*m2ahat))/(lam1^2*(s1ahat^2*s2ahat^2 - cova^2)^3) - (lam2*((6*lam1 + 6*lam1^2*m1ahat - 6*W1a^lam1*lam1 - 6*W2a^lam2*lam1 + 6*W1a^lam1*W2a^lam2*lam1 - 6*W2a^lam2*lam1^2*m1ahat)*cova^2*s1ahat^2*s2ahat + (8*W2a^lam2*lam1^2*m2ahat - 8*lam1^2*m2ahat)*cova*s1ahat^4*s2ahat + (2*lam1 + 2*lam1^2*m1ahat - 2*W1a^lam1*lam1 - 2*W2a^lam2*lam1 + 2*W1a^lam1*W2a^lam2*lam1 - 2*W2a^lam2*lam1^2*m1ahat)*s1ahat^4*s2ahat^3) - cova*s1ahat^4*s2ahat*(4*W2a^(2*lam2)*lam1^2 - 8*W2a^lam2*lam1^2 + 4*lam1^2))/(lam1^2*lam2^2*(s1ahat^2*s2ahat^2 - cova^2)^3))
h9_7=0
h9_8=0
h9_9=sum(- (s1ahat^2*(cova^2*(6*lam2*m2ahat - 6*W2a^lam2 + 3*W2a^(2*lam2) + 3*lam2^2*m2ahat^2 - 6*W2a^lam2*lam2*m2ahat + 3) - cova*(6*lam2*m1ahat*s2ahat^2 - 6*W2a^lam2*lam2*m1ahat*s2ahat^2 + 6*lam2^2*m1ahat*m2ahat*s2ahat^2) + lam2^2*m1ahat^2*s2ahat^4) - cova^3*(2*lam2*m1ahat + 2*lam2^2*m1ahat*m2ahat - 2*W2a^lam2*lam2*m1ahat) + s1ahat^4*(W2a^(2*lam2)*s2ahat^2 - 2*W2a^lam2*s2ahat^2 + s2ahat^2 - lam2^2*s2ahat^4 + 2*lam2*m2ahat*s2ahat^2 + lam2^2*m2ahat^2*s2ahat^2 - 2*W2a^lam2*lam2*m2ahat*s2ahat^2) + cova^4*lam2^2 + 3*cova^2*lam2^2*m1ahat^2*s2ahat^2)/(lam2^2*(s1ahat^2*s2ahat^2 - cova^2)^3) - (cova^2*(3*lam2^2*s2ahat^2 - 6*W1a^lam1*lam2^2*s2ahat^2 + 3*W1a^(2*lam1)*lam2^2*s2ahat^2) + s1ahat^2*(lam2^2*s2ahat^4 - 2*W1a^lam1*lam2^2*s2ahat^4 + W1a^(2*lam1)*lam2^2*s2ahat^4) - lam1*(cova^3*(2*lam2 + 2*lam2^2*m2ahat - 2*W1a^lam1*lam2 - 2*W2a^lam2*lam2 + 2*W1a^lam1*W2a^lam2*lam2 - 2*W1a^lam1*lam2^2*m2ahat) - cova^2*(6*lam2^2*m1ahat*s2ahat^2 - 6*W1a^lam1*lam2^2*m1ahat*s2ahat^2) + s1ahat^2*(cova*(6*lam2*s2ahat^2 + 6*lam2^2*m2ahat*s2ahat^2 - 6*W1a^lam1*lam2*s2ahat^2 - 6*W2a^lam2*lam2*s2ahat^2 - 6*W1a^lam1*lam2^2*m2ahat*s2ahat^2 + 6*W1a^lam1*W2a^lam2*lam2*s2ahat^2) - 2*lam2^2*m1ahat*s2ahat^4 + 2*W1a^lam1*lam2^2*m1ahat*s2ahat^4)))/(lam1^2*lam2^2*(s1ahat^2*s2ahat^2 - cova^2)^3))
h9_10=0
h9_11=sum(-((W1a^lam1*lam1*log(W1a) - W1a^lam1 + 1)*(cova^2*lam1 - 2*cova*lam2*s2ahat^2 + lam1*s1ahat^2*s2ahat^2 - W2a^lam2*cova^2*lam1 + 2*W1a^lam1*cova*lam2*s2ahat^2 + cova^2*lam1*lam2*m2ahat - W2a^lam2*lam1*s1ahat^2*s2ahat^2 + lam1*lam2*m2ahat*s1ahat^2*s2ahat^2 - 2*cova*lam1*lam2*m1ahat*s2ahat^2))/(lam1^3*lam2*(s1ahat^2*s2ahat^2 - cova^2)^2))
h9_12=sum(-((W2a^lam2*lam2*log(W2a) - W2a^lam2 + 1)*(cova^2*lam2 - 2*cova*lam1*s1ahat^2 + lam2*s1ahat^2*s2ahat^2 - W1a^lam1*cova^2*lam2 + 2*W2a^lam2*cova*lam1*s1ahat^2 + cova^2*lam1*lam2*m1ahat - W1a^lam1*lam2*s1ahat^2*s2ahat^2 + lam1*lam2*m1ahat*s1ahat^2*s2ahat^2 - 2*cova*lam1*lam2*m2ahat*s1ahat^2))/(lam1*lam2^3*(s1ahat^2*s2ahat^2 - cova^2)^2))
h10_1=0
h10_2=0
h10_3=sum((covb^2*(lam2*m2bhat - W2b^lam2 + 1) + s2bhat^2*((lam2*m2bhat - W2b^lam2 + 1)*s1bhat^2 - 2*covb*lam2*m1bhat))/(lam2*(s1bhat^2*s2bhat^2 - covb^2)^2) - (covb*s2bhat^2*(2*lam2 - 2*W1b^lam1*lam2))/(lam1*lam2*(s1bhat^2*s2bhat^2 - covb^2)^2))
h10_4=sum(- (covb*(s2bhat^4*(- 4*lam2^2*m1bhat^2*s1bhat + 2*lam2^2*s1bhat^3) - s1bhat^3*s2bhat^2*(4*lam2*m2bhat - 4*W2b^lam2 + 2*W2b^(2*lam2) + 2*lam2^2*m2bhat^2 - 4*W2b^lam2*lam2*m2bhat + 2)) - covb^3*(s1bhat*(4*lam2*m2bhat - 4*W2b^lam2 + 2*W2b^(2*lam2) + 2*lam2^2*m2bhat^2 - 4*W2b^lam2*lam2*m2bhat + 2) + 2*lam2^2*s1bhat*s2bhat^2) + s1bhat^3*s2bhat^4*(2*lam2*m1bhat + 2*lam2^2*m1bhat*m2bhat - 2*W2b^lam2*lam2*m1bhat) + covb^2*s1bhat*s2bhat^2*(6*lam2*m1bhat + 6*lam2^2*m1bhat*m2bhat - 6*W2b^lam2*lam2*m1bhat))/(lam2^2*(s1bhat^2*s2bhat^2 - covb^2)^3) - (lam1*((6*lam2 + 6*lam2^2*m2bhat - 6*W1b^lam1*lam2 - 6*W2b^lam2*lam2 + 6*W1b^lam1*W2b^lam2*lam2 - 6*W1b^lam1*lam2^2*m2bhat)*covb^2*s1bhat*s2bhat^2 + (8*W1b^lam1*lam2^2*m1bhat - 8*lam2^2*m1bhat)*covb*s1bhat*s2bhat^4 + (2*lam2 + 2*lam2^2*m2bhat - 2*W1b^lam1*lam2 - 2*W2b^lam2*lam2 + 2*W1b^lam1*W2b^lam2*lam2 - 2*W1b^lam1*lam2^2*m2bhat)*s1bhat^3*s2bhat^4) - covb*s1bhat*s2bhat^4*(4*W1b^(2*lam1)*lam2^2 - 8*W1b^lam1*lam2^2 + 4*lam2^2))/(lam1^2*lam2^2*(s1bhat^2*s2bhat^2 - covb^2)^3))
h10_5=0
h10_6=0
h10_7=sum((covb^2*(lam1*m1bhat - W1b^lam1 + 1) + s1bhat^2*((lam1*m1bhat - W1b^lam1 + 1)*s2bhat^2 - 2*covb*lam1*m2bhat))/(lam1*(s1bhat^2*s2bhat^2 - covb^2)^2) - (covb*s1bhat^2*(2*lam1 - 2*W2b^lam2*lam1))/(lam1*lam2*(s1bhat^2*s2bhat^2 - covb^2)^2))
h10_8=sum(- (covb*(s1bhat^4*(- 4*lam1^2*m2bhat^2*s2bhat + 2*lam1^2*s2bhat^3) - s1bhat^2*s2bhat^3*(4*lam1*m1bhat - 4*W1b^lam1 + 2*W1b^(2*lam1) + 2*lam1^2*m1bhat^2 - 4*W1b^lam1*lam1*m1bhat + 2)) - covb^3*(s2bhat*(4*lam1*m1bhat - 4*W1b^lam1 + 2*W1b^(2*lam1) + 2*lam1^2*m1bhat^2 - 4*W1b^lam1*lam1*m1bhat + 2) + 2*lam1^2*s1bhat^2*s2bhat) + s1bhat^4*s2bhat^3*(2*lam1*m2bhat + 2*lam1^2*m1bhat*m2bhat - 2*W1b^lam1*lam1*m2bhat) + covb^2*s1bhat^2*s2bhat*(6*lam1*m2bhat + 6*lam1^2*m1bhat*m2bhat - 6*W1b^lam1*lam1*m2bhat))/(lam1^2*(s1bhat^2*s2bhat^2 - covb^2)^3) - (lam2*((6*lam1 + 6*lam1^2*m1bhat - 6*W1b^lam1*lam1 - 6*W2b^lam2*lam1 + 6*W1b^lam1*W2b^lam2*lam1 - 6*W2b^lam2*lam1^2*m1bhat)*covb^2*s1bhat^2*s2bhat + (8*W2b^lam2*lam1^2*m2bhat - 8*lam1^2*m2bhat)*covb*s1bhat^4*s2bhat + (2*lam1 + 2*lam1^2*m1bhat - 2*W1b^lam1*lam1 - 2*W2b^lam2*lam1 + 2*W1b^lam1*W2b^lam2*lam1 - 2*W2b^lam2*lam1^2*m1bhat)*s1bhat^4*s2bhat^3) - covb*s1bhat^4*s2bhat*(4*W2b^(2*lam2)*lam1^2 - 8*W2b^lam2*lam1^2 + 4*lam1^2))/(lam1^2*lam2^2*(s1bhat^2*s2bhat^2 - covb^2)^3))
h10_9=0
h10_10=sum(- (s1bhat^2*(covb^2*(6*lam2*m2bhat - 6*W2b^lam2 + 3*W2b^(2*lam2) + 3*lam2^2*m2bhat^2 - 6*W2b^lam2*lam2*m2bhat + 3) - covb*(6*lam2*m1bhat*s2bhat^2 - 6*W2b^lam2*lam2*m1bhat*s2bhat^2 + 6*lam2^2*m1bhat*m2bhat*s2bhat^2) + lam2^2*m1bhat^2*s2bhat^4) - covb^3*(2*lam2*m1bhat + 2*lam2^2*m1bhat*m2bhat - 2*W2b^lam2*lam2*m1bhat) + s1bhat^4*(W2b^(2*lam2)*s2bhat^2 - 2*W2b^lam2*s2bhat^2 + s2bhat^2 - lam2^2*s2bhat^4 + 2*lam2*m2bhat*s2bhat^2 + lam2^2*m2bhat^2*s2bhat^2 - 2*W2b^lam2*lam2*m2bhat*s2bhat^2) + covb^4*lam2^2 + 3*covb^2*lam2^2*m1bhat^2*s2bhat^2)/(lam2^2*(s1bhat^2*s2bhat^2 - covb^2)^3) - (covb^2*(3*lam2^2*s2bhat^2 - 6*W1b^lam1*lam2^2*s2bhat^2 + 3*W1b^(2*lam1)*lam2^2*s2bhat^2) + s1bhat^2*(lam2^2*s2bhat^4 - 2*W1b^lam1*lam2^2*s2bhat^4 + W1b^(2*lam1)*lam2^2*s2bhat^4) - lam1*(covb^3*(2*lam2 + 2*lam2^2*m2bhat - 2*W1b^lam1*lam2 - 2*W2b^lam2*lam2 + 2*W1b^lam1*W2b^lam2*lam2 - 2*W1b^lam1*lam2^2*m2bhat) - covb^2*(6*lam2^2*m1bhat*s2bhat^2 - 6*W1b^lam1*lam2^2*m1bhat*s2bhat^2) + s1bhat^2*(covb*(6*lam2*s2bhat^2 + 6*lam2^2*m2bhat*s2bhat^2 - 6*W1b^lam1*lam2*s2bhat^2 - 6*W2b^lam2*lam2*s2bhat^2 - 6*W1b^lam1*lam2^2*m2bhat*s2bhat^2 + 6*W1b^lam1*W2b^lam2*lam2*s2bhat^2) - 2*lam2^2*m1bhat*s2bhat^4 + 2*W1b^lam1*lam2^2*m1bhat*s2bhat^4)))/(lam1^2*lam2^2*(s1bhat^2*s2bhat^2 - covb^2)^3))
h10_11=sum(-((W1b^lam1*lam1*log(W1b) - W1b^lam1 + 1)*(covb^2*lam1 - 2*covb*lam2*s2bhat^2 + lam1*s1bhat^2*s2bhat^2 - W2b^lam2*covb^2*lam1 + 2*W1b^lam1*covb*lam2*s2bhat^2 + covb^2*lam1*lam2*m2bhat - W2b^lam2*lam1*s1bhat^2*s2bhat^2 + lam1*lam2*m2bhat*s1bhat^2*s2bhat^2 - 2*covb*lam1*lam2*m1bhat*s2bhat^2))/(lam1^3*lam2*(s1bhat^2*s2bhat^2 - covb^2)^2))
h10_12=sum(-((W2b^lam2*lam2*log(W2b) - W2b^lam2 + 1)*(covb^2*lam2 - 2*covb*lam1*s1bhat^2 + lam2*s1bhat^2*s2bhat^2 - W1b^lam1*covb^2*lam2 + 2*W2b^lam2*covb*lam1*s1bhat^2 + covb^2*lam1*lam2*m1bhat - W1b^lam1*lam2*s1bhat^2*s2bhat^2 + lam1*lam2*m1bhat*s1bhat^2*s2bhat^2 - 2*covb*lam1*lam2*m2bhat*s1bhat^2))/(lam1*lam2^3*(s1bhat^2*s2bhat^2 - covb^2)^2))
h11_1=sum((2*s2ahat^2*(W1a^lam1*lam1*log(W1a) - W1a^lam1 + 1))/(lam1^2*(- 2*cova^2 + 2*s1ahat^2*s2ahat^2)))
h11_2=sum((2*s1ahat*s2ahat^2*(W1a^lam1*lam1*log(W1a) - W1a^lam1 + 1)*(cova*lam1 - lam2*s2ahat^2 - W2a^lam2*cova*lam1 + W1a^lam1*lam2*s2ahat^2 + cova*lam1*lam2*m2ahat - lam1*lam2*m1ahat*s2ahat^2))/(lam1^3*lam2*(s1ahat^2*s2ahat^2 - cova^2)^2))
h11_3=sum((2*s2bhat^2*(W1b^lam1*lam1*log(W1b) - W1b^lam1 + 1))/(lam1^2*(- 2*covb^2 + 2*s1bhat^2*s2bhat^2)))
h11_4=sum((2*s1bhat*s2bhat^2*(W1b^lam1*lam1*log(W1b) - W1b^lam1 + 1)*(covb*lam1 - lam2*s2bhat^2 - W2b^lam2*covb*lam1 + W1b^lam1*lam2*s2bhat^2 + covb*lam1*lam2*m2bhat - lam1*lam2*m1bhat*s2bhat^2))/(lam1^3*lam2*(s1bhat^2*s2bhat^2 - covb^2)^2))
h11_5=sum(-(cova*(W1a^lam1*lam1*log(W1a) - W1a^lam1 + 1))/(lam1^2*(- cova^2 + s1ahat^2*s2ahat^2)))
h11_6=sum(-(2*cova*s2ahat*(W1a^lam1*lam1*log(W1a) - W1a^lam1 + 1)*(cova*lam2 - lam1*s1ahat^2 - W1a^lam1*cova*lam2 + W2a^lam2*lam1*s1ahat^2 + cova*lam1*lam2*m1ahat - lam1*lam2*m2ahat*s1ahat^2))/(lam1^3*lam2*(s1ahat^2*s2ahat^2 - cova^2)^2))
h11_7=sum(-(covb*(W1b^lam1*lam1*log(W1b) - W1b^lam1 + 1))/(lam1^2*(- covb^2 + s1bhat^2*s2bhat^2)))
h11_8=sum(-(2*covb*s2bhat*(W1b^lam1*lam1*log(W1b) - W1b^lam1 + 1)*(covb*lam2 - lam1*s1bhat^2 - W1b^lam1*covb*lam2 + W2b^lam2*lam1*s1bhat^2 + covb*lam1*lam2*m1bhat - lam1*lam2*m2bhat*s1bhat^2))/(lam1^3*lam2*(s1bhat^2*s2bhat^2 - covb^2)^2))
h11_9=sum(-((W1a^lam1*lam1*log(W1a) - W1a^lam1 + 1)*(cova^2*lam1 - 2*cova*lam2*s2ahat^2 + lam1*s1ahat^2*s2ahat^2 - W2a^lam2*cova^2*lam1 + 2*W1a^lam1*cova*lam2*s2ahat^2 + cova^2*lam1*lam2*m2ahat - W2a^lam2*lam1*s1ahat^2*s2ahat^2 + lam1*lam2*m2ahat*s1ahat^2*s2ahat^2 - 2*cova*lam1*lam2*m1ahat*s2ahat^2))/(lam1^3*lam2*(s1ahat^2*s2ahat^2 - cova^2)^2))
h11_10=sum(-((W1b^lam1*lam1*log(W1b) - W1b^lam1 + 1)*(covb^2*lam1 - 2*covb*lam2*s2bhat^2 + lam1*s1bhat^2*s2bhat^2 - W2b^lam2*covb^2*lam1 + 2*W1b^lam1*covb*lam2*s2bhat^2 + covb^2*lam1*lam2*m2bhat - W2b^lam2*lam1*s1bhat^2*s2bhat^2 + lam1*lam2*m2bhat*s1bhat^2*s2bhat^2 - 2*covb*lam1*lam2*m1bhat*s2bhat^2))/(lam1^3*lam2*(s1bhat^2*s2bhat^2 - covb^2)^2))
h11_11=sum(((2*((W1a^lam1-1)/lam1^2-(W1a^lam1*log(W1a))/lam1)^2)/s1ahat^2-(2*(m1ahat-(W1a^lam1-1)/lam1)*((2*W1a^lam1-2)/lam1^3-(2*W1a^lam1*log(W1a))/lam1^2+(W1a^lam1*log(W1a)^2)/lam1))/s1ahat^2+(2*cova*(m2ahat-(W2a^lam2-1)/lam2)*((2*W1a^lam1-2)/lam1^3-(2*W1a^lam1*log(W1a))/lam1^2+(W1a^lam1*log(W1a)^2)/lam1))/(s1ahat^2*s2ahat^2))/((2*cova^2)/(s1ahat^2*s2ahat^2)-2))+sum(((2*((W1b^lam1-1)/lam1^2-(W1b^lam1*log(W1b))/lam1)^2)/s1bhat^2-(2*(m1bhat-(W1b^lam1-1)/lam1)*((2*W1b^lam1-2)/lam1^3-(2*W1b^lam1*log(W1b))/lam1^2+(W1b^lam1*log(W1b)^2)/lam1))/s1bhat^2+(2*covb*(m2bhat-(W2b^lam2-1)/lam2)*((2*W1b^lam1-2)/lam1^3-(2*W1b^lam1*log(W1b))/lam1^2+(W1b^lam1*log(W1b)^2)/lam1))/(s1bhat^2*s2bhat^2))/((2*covb^2)/(s1bhat^2*s2bhat^2)-2))
h11_12=sum((2*cova*(W1a^lam1*lam1*log(W1a)-W1a^lam1+1)*(W2a^lam2*lam2*log(W2a)-W2a^lam2+1))/(lam1^2*lam2^2*(-2*cova^2+2*s1ahat^2*s2ahat^2)))+sum((2*covb*(W1b^lam1*lam1*log(W1b)-W1b^lam1+1)*(W2b^lam2*lam2*log(W2b)-W2b^lam2+1))/(lam1^2*lam2^2*(-2*covb^2+2*s1bhat^2*s2bhat^2)))
h12_1=sum(-(cova*(W2a^lam2*lam2*log(W2a) - W2a^lam2 + 1))/(lam2^2*(- cova^2 + s1ahat^2*s2ahat^2)))
h12_2=sum(-(2*cova*s1ahat*(W2a^lam2*lam2*log(W2a) - W2a^lam2 + 1)*(cova*lam1 - lam2*s2ahat^2 - W2a^lam2*cova*lam1 + W1a^lam1*lam2*s2ahat^2 + cova*lam1*lam2*m2ahat - lam1*lam2*m1ahat*s2ahat^2))/(lam1*lam2^3*(s1ahat^2*s2ahat^2 - cova^2)^2))
h12_3=sum(-(covb*(W2b^lam2*lam2*log(W2b) - W2b^lam2 + 1))/(lam2^2*(- covb^2 + s1bhat^2*s2bhat^2)))
h12_4=sum(-(2*covb*s1bhat*(W2b^lam2*lam2*log(W2b) - W2b^lam2 + 1)*(covb*lam1 - lam2*s2bhat^2 - W2b^lam2*covb*lam1 + W1b^lam1*lam2*s2bhat^2 + covb*lam1*lam2*m2bhat - lam1*lam2*m1bhat*s2bhat^2))/(lam1*lam2^3*(s1bhat^2*s2bhat^2 - covb^2)^2))
h12_5=sum((2*s1ahat^2*(W2a^lam2*lam2*log(W2a) - W2a^lam2 + 1))/(lam2^2*(- 2*cova^2 + 2*s1ahat^2*s2ahat^2)))
h12_6=sum((2*s1ahat^2*s2ahat*(W2a^lam2*lam2*log(W2a) - W2a^lam2 + 1)*(cova*lam2 - lam1*s1ahat^2 - W1a^lam1*cova*lam2 + W2a^lam2*lam1*s1ahat^2 + cova*lam1*lam2*m1ahat - lam1*lam2*m2ahat*s1ahat^2))/(lam1*lam2^3*(s1ahat^2*s2ahat^2 - cova^2)^2))
h12_7=sum((2*s1bhat^2*(W2b^lam2*lam2*log(W2b) - W2b^lam2 + 1))/(lam2^2*(- 2*covb^2 + 2*s1bhat^2*s2bhat^2)))
h12_8=sum((2*s1bhat^2*s2bhat*(W2b^lam2*lam2*log(W2b) - W2b^lam2 + 1)*(covb*lam2 - lam1*s1bhat^2 - W1b^lam1*covb*lam2 + W2b^lam2*lam1*s1bhat^2 + covb*lam1*lam2*m1bhat - lam1*lam2*m2bhat*s1bhat^2))/(lam1*lam2^3*(s1bhat^2*s2bhat^2 - covb^2)^2))
h12_9=sum(-((W2a^lam2*lam2*log(W2a) - W2a^lam2 + 1)*(cova^2*lam2 - 2*cova*lam1*s1ahat^2 + lam2*s1ahat^2*s2ahat^2 - W1a^lam1*cova^2*lam2 + 2*W2a^lam2*cova*lam1*s1ahat^2 + cova^2*lam1*lam2*m1ahat - W1a^lam1*lam2*s1ahat^2*s2ahat^2 + lam1*lam2*m1ahat*s1ahat^2*s2ahat^2 - 2*cova*lam1*lam2*m2ahat*s1ahat^2))/(lam1*lam2^3*(s1ahat^2*s2ahat^2 - cova^2)^2))
h12_10=sum(-((W2b^lam2*lam2*log(W2b) - W2b^lam2 + 1)*(covb^2*lam2 - 2*covb*lam1*s1bhat^2 + lam2*s1bhat^2*s2bhat^2 - W1b^lam1*covb^2*lam2 + 2*W2b^lam2*covb*lam1*s1bhat^2 + covb^2*lam1*lam2*m1bhat - W1b^lam1*lam2*s1bhat^2*s2bhat^2 + lam1*lam2*m1bhat*s1bhat^2*s2bhat^2 - 2*covb*lam1*lam2*m2bhat*s1bhat^2))/(lam1*lam2^3*(s1bhat^2*s2bhat^2 - covb^2)^2))
h12_11=sum((2*cova*(W1a^lam1*lam1*log(W1a)-W1a^lam1+1)*(W2a^lam2*lam2*log(W2a)-W2a^lam2+1))/(lam1^2*lam2^2*(-2*cova^2+2*s1ahat^2*s2ahat^2)))+sum((2*covb*(W1b^lam1*lam1*log(W1b)-W1b^lam1+1)*(W2b^lam2*lam2*log(W2b)-W2b^lam2+1))/(lam1^2*lam2^2*(-2*covb^2+2*s1bhat^2*s2bhat^2)))
h12_12=sum(((2*((W2a^lam2-1)/lam2^2-(W2a^lam2*log(W2a))/lam2)^2)/s2ahat^2-(2*(m2ahat-(W2a^lam2-1)/lam2)*((2*W2a^lam2-2)/lam2^3-(2*W2a^lam2*log(W2a))/lam2^2+(W2a^lam2*log(W2a)^2)/lam2))/s2ahat^2+(2*cova*(m1ahat-(W1a^lam1-1)/lam1)*((2*W2a^lam2-2)/lam2^3-(2*W2a^lam2*log(W2a))/lam2^2+(W2a^lam2*log(W2a)^2)/lam2))/(s1ahat^2*s2ahat^2))/((2*cova^2)/(s1ahat^2*s2ahat^2)-2))+sum(((2*((W2b^lam2-1)/lam2^2-(W2b^lam2*log(W2b))/lam2)^2)/s2bhat^2-(2*(m2bhat-(W2b^lam2-1)/lam2)*((2*W2b^lam2-2)/lam2^3-(2*W2b^lam2*log(W2b))/lam2^2+(W2b^lam2*log(W2b)^2)/lam2))/s2bhat^2+(2*covb*(m1bhat-(W1b^lam1-1)/lam1)*((2*W2b^lam2-2)/lam2^3-(2*W2b^lam2*log(W2b))/lam2^2+(W2b^lam2*log(W2b)^2)/lam2))/(s1bhat^2*s2bhat^2))/((2*covb^2)/(s1bhat^2*s2bhat^2)-2))
h1r=c(h1_1,h1_2,h1_3,h1_4,h1_5,h1_6,h1_7,h1_8,h1_9,h1_10,h1_11,h1_12)
h2r=c(h2_1,h2_2,h2_3,h2_4,h2_5,h2_6,h2_7,h2_8,h2_9,h2_10,h2_11,h2_12)
h3r=c(h3_1,h3_2,h3_3,h3_4,h3_5,h3_6,h3_7,h3_8,h3_9,h3_10,h3_11,h3_12)
h4r=c(h4_1,h4_2,h4_3,h4_4,h4_5,h4_6,h4_7,h4_8,h4_9,h4_10,h4_11,h4_12)
h5r=c(h5_1,h5_2,h5_3,h5_4,h5_5,h5_6,h5_7,h5_8,h5_9,h5_10,h5_11,h5_12)
h6r=c(h6_1,h6_2,h6_3,h6_4,h6_5,h6_6,h6_7,h6_8,h6_9,h6_10,h6_11,h6_12)
h7r=c(h7_1,h7_2,h7_3,h7_4,h7_5,h7_6,h7_7,h7_8,h7_9,h7_10,h7_11,h7_12)
h8r=c(h8_1,h8_2,h8_3,h8_4,h8_5,h8_6,h8_7,h8_8,h8_9,h8_10,h8_11,h8_12)
h9r=c(h9_1,h9_2,h9_3,h9_4,h9_5,h9_6,h9_7,h9_8,h9_9,h9_10,h9_11,h9_12)
h10r=c(h10_1,h10_2,h10_3,h10_4,h10_5,h10_6,h10_7,h10_8,h10_9,h10_10,h10_11,h10_12)
h11r=c(h11_1,h11_2,h11_3,h11_4,h11_5,h11_6,h11_7,h11_8,h11_9,h11_10,h11_11,h11_12)
h12r=c(h12_1,h12_2,h12_3,h12_4,h12_5,h12_6,h12_7,h12_8,h12_9,h12_10,h12_11,h12_12)
HCF=rbind(h1r,h2r,h3r,h4r,h5r,h6r,h7r,h8r,h9r,h10r,h11r,h12r)
HCF[1,2]=0;HCF[1,6]=0;HCF[1,9]=0;
HCF[2,1]=0;HCF[2,5]=0;
HCF[3,4]=0;HCF[3,8]=0;HCF[3,10]=0;HCF[4,3]=0;
HCF[4,7]=0;
HCF[5,2]=0;HCF[5,6]=0;HCF[5,9]=0;
HCF[6,1]=0;HCF[6,5]=0;
HCF[7,4]=0;HCF[7,8]=0;HCF[7,10]=0;
HCF[8,3]=0;HCF[8,7]=0;
HCF[9,1]=0;HCF[9,5]=0;
HCF[10,3]=0;
HCF[10,7]=0;
HCF[1,1]=na*HCF[1,1];
HCF[1,5]=na*HCF[1,5];
HCF[5,1]=na*HCF[5,1];
HCF[3,3]=nb*HCF[3,3];
HCF[3,7]=nb*HCF[3,7];
HCF[7,3]=nb*HCF[7,3];
HCF[7,7]=nb*HCF[7,7];
HCF[5,5]=na*HCF[5,5];
VV=inv(-HCF)
att=1-atSpec;
t=att;
s1alam=sqrt( var(W1alam)*(na-1)/na )
s1blam=sqrt( var(W1blam)*(nb-1)/nb )
s2alam=sqrt( var(W2alam)*(na-1)/na )
s2blam=sqrt( var(W2blam)*(nb-1)/nb )
m1alam=mean(W1alam)
m1blam=mean(W1blam)
m2alam=mean(W2alam)
m2blam=mean(W2blam)
#############################################################################
#############################################################################
#==================AT A GIVEN T FOR THE ORIGINAL===================
ROC1or= (1-pnorm(qnorm(1-t,m1alam,s1alam),m1blam,s1blam));
ROC2or= (1-pnorm(qnorm(1-t,m2alam,s2alam),m2blam,s2blam));
dm1a = -(2^(1/2)*exp(-((m1ahat - m1bhat)/s1bhat + (2^(1/2)*s1ahat*erfcinv(2*t))/s1bhat)^2/2))/(2*s1bhat*pi^(1/2))
dm1b = (2^(1/2)*exp(-((m1ahat - m1bhat)/s1bhat + (2^(1/2)*s1ahat*erfcinv(2*t))/s1bhat)^2/2))/(2*s1bhat*pi^(1/2))
ds1a = -(erfcinv(2*t)*exp(-((m1ahat - m1bhat)/s1bhat + (2^(1/2)*s1ahat*erfcinv(2*t))/s1bhat)^2/2))/(s1bhat*pi^(1/2))
ds1b = (2^(1/2)*exp(-((m1ahat - m1bhat)/s1bhat + (2^(1/2)*s1ahat*erfcinv(2*t))/s1bhat)^2/2)*((m1ahat - m1bhat)/s1bhat^2 + (2^(1/2)*s1ahat*erfcinv(2*t))/s1bhat^2))/(2*pi^(1/2))
der1=c(dm1a, ds1a, dm1b, ds1b);
vROC1=c(der1)%*%VV[1:4,1:4]%*%t(t(c(der1)))
vROC1=c(vROC1)
dm2a = -(2^(1/2)*exp(-((m2ahat - m2bhat)/s2bhat + (2^(1/2)*s2ahat*erfcinv(2*t))/s2bhat)^2/2))/(2*s2bhat*pi^(1/2))
dm2b = (2^(1/2)*exp(-((m2ahat - m2bhat)/s2bhat + (2^(1/2)*s2ahat*erfcinv(2*t))/s2bhat)^2/2))/(2*s2bhat*pi^(1/2))
ds2a = -(erfcinv(2*t)*exp(-((m2ahat - m2bhat)/s2bhat + (2^(1/2)*s2ahat*erfcinv(2*t))/s2bhat)^2/2))/(s2bhat*pi^(1/2))
ds2b = (2^(1/2)*exp(-((m2ahat - m2bhat)/s2bhat + (2^(1/2)*s2ahat*erfcinv(2*t))/s2bhat)^2/2)*((m2ahat - m2bhat)/s2bhat^2 + (2^(1/2)*s2ahat*erfcinv(2*t))/s2bhat^2))/(2*pi^(1/2))
der2=c(dm2a, ds2a, dm2b, ds2b);
vROC2=c(der2)%*%VV[5:8,5:8]%*%t(t(c(der2)))
covROC=c(der1, 0, 0, 0, 0)%*%VV[1:8,1:8]%*%t(t(c(0, 0, 0, 0, der2)))
Z=(ROC2or-ROC1or)/(sqrt(vROC2+vROC1-2*covROC))
SEZ=sqrt(vROC1+vROC2-2*covROC)
pval2tZ=2*pnorm(-abs(Z));
CIoriginal=c((ROC2or-ROC1or)-pmalpha*SEZ, (ROC2or-ROC1or)+pmalpha*SEZ)
CIoriginal
#############################################################################
#############################################################################
#==================AT A GIVEN T FOR THE TRANSORMED===================
ROC1MC= qnorm(1-pnorm(qnorm(1-t,m1alam,s1alam),m1blam,s1blam));
ROC2MC= qnorm(1-pnorm(qnorm(1-t,m2alam,s2alam),m2blam,s2blam));
dm1a =-1/s1blam;
ds1a =-(2^(1/2)*erfcinv(2*t))/s1blam;
dm1b =1/s1blam;
ds1b =(m1alam - m1blam)/s1blam^2 + (2^(1/2)*s1alam*erfcinv(2*t))/s1blam^2;
der1=c(dm1a, ds1a, dm1b, ds1b);
vROC1=c(der1)%*%VV[1:4,1:4]%*%t(t(c(der1)))
vROC1=c(vROC1)
dm2a =-1/s2blam;
ds2a =-(2^(1/2)*erfcinv(2*t))/s2blam;
dm2b =1/s2blam;
ds2b =(m2alam - m2blam)/s2blam^2 + (2^(1/2)*s2alam*erfcinv(2*t))/s2blam^2;
der2=c(dm2a, ds2a, dm2b, ds2b);
vROC2=c(der2)%*%VV[5:8,5:8]%*%t(t(c(der2)))
vROC2=c(vROC2)
covROC=c(der1, 0, 0, 0, 0)%*%VV[1:8,1:8]%*%t(t(c(0, 0, 0, 0, der2)))
covROC=c(covROC)
Zroct=(ROC2MC-ROC1MC)/(sqrt(vROC2+vROC1-2*covROC))
SEroct=sqrt(vROC1+vROC2-2*covROC)
pval2t_roct=2*pnorm(-abs(Zroct));pval2t=pval2t_roct;
CIZstar=c((ROC2MC-ROC1MC)-pmalpha*SEroct, (ROC2MC-ROC1MC)+pmalpha*SEroct)
CIZstar
SE1=pnorm(ROC1MC)
SE2=pnorm(ROC2MC)
res <- data.frame(Sens1 = formattable(SE1, digits = 4, format = "f"),
Sens2 = formattable(SE2, digits = 4, format = "f"),
diff = formattable(SE2 - SE1, digits = 4, format = "f"),
p_value_probit = formattable(pval2t, digits = 4, format = "f"),
p_value = formattable(pval2tZ, digits = 4, format = "f"),
ci_ll = formattable(CIoriginal[1], digits = 4, format = "f"),
ci_ul = formattable(CIoriginal[2], digits = 4, format = "f"))
rownames(res) <- c("Estimates:")
colnames(res) <- c("Sens 1", "Sens 2", "Diff", "P-Val (Probit)", "P-Val", "CI (LL)", "CI (UL)")
res <- formattable(as.matrix(res), digits = 4, format = "f")
res
#====TWO ROC FUNCTIONS: ONE IS THE BOXCOX AND THE OTHER THE REFERENCE LINE================
roc1<-function(t){
na = length(W1alam)
nb = length(W1blam)
s1alam=sqrt( var(W1alam)*(na-1)/na )
s1blam=sqrt( var(W1blam)*(nb-1)/nb )
1-pnorm(qnorm(1-t,
mean=mean(W1alam),
sd=s1alam),
mean=mean(W1blam),
sd=s1blam)
}
roc2<-function(t){
na = length(W2alam)
nb = length(W2blam)
s2alam=sqrt( var(W2alam)*(na-1)/na )
s2blam=sqrt( var(W2blam)*(nb-1)/nb )
1-pnorm(qnorm(1-t,
mean=mean(W2alam),
sd=s2alam),
mean=mean(W2blam),
sd=s2blam)
}
rocuseless<-function(t){
1-pnorm(qnorm(1-t,mean=1,sd=1),mean=1,sd=1)
}
#================IF PLOTS ARE REQUESTED PLOT THE ROCS==
if (plots=="on") {
# x11()
plot(linspace(0,1,1000),roc1(linspace(0,1,1000)),main="Box-Cox Based ROCs",xlab="FPR = 1 - Specificity",ylab="TPR = Sensitivity",type="l",col="red")
lines(linspace(0,1,1000),roc2(linspace(0,1,1000)),main="Box-Cox Based ROCs",xlab="FPR = 1 - Specificity",ylab="TPR = Sensitivity",type="l",col="black")
lines(linspace(0,1,1000),linspace(0,1,1000),type="l", lty=2)
points(1-atSpec,SE1, col = "red")
points(1-atSpec,SE2, col = "black")
legend("bottomright", legend=c(paste("ROC for Marker 1 with Sens =", formattable(SE1, digits = 4, format = "f"), "at Spec =", formattable(atSpec, digits = 4, format = "f")),
paste("ROC for Marker 2 with Sens =", formattable(SE2, digits = 4, format = "f"), "at Spec =", formattable(atSpec, digits = 4, format = "f")),
paste("P-value for the Sens difference:", formattable(pval2t, digits = 4, format = "f"))),
col=c("red", "black", "white"), lty=c(1, 1, NA), pch = c(NA, NA, NA), cex=0.8)
}
return(list(resultstable=res,
Sens1=SE1,
Sens2=SE2,
pvalue_probit_difference= pval2t,
pvalue_difference= pval2tZ,
CI_difference= CIoriginal,
roc1=roc1,
roc2=roc2,
transx1=W1alam,
transy1=W1blam,
transformation.parameter.1 = lam[1],
transx2=W2alam,
transy2=W2blam,
transformation.parameter.2 = lam[2]))
}
}
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.