Nothing
momeff.nonzeromean=function(tstat,n1,n2=n1,gamma2,lower.df=6.1,upper.df=100,approx=TRUE) {
#################### non-zero mean effect size case
### some pre-computed constants
v=n1+n2-2
vv=v*v
et=mean(tstat)
et2=mean(tstat^2)
et3=mean(tstat^3)
etet=et*et
et2et2=et2^2
gv3=gamma(v/2-3/2)
sgv3=sqrt(gv3)
gv1=gamma(v/2-1/2)
gv=gamma(v/2)
gvgv=gv*gv
gv1gv1=gv1*gv1
#### after eliminating mu and sigma2, need to solve the following quadratic function of pi0
a=(
-3*et*gv3*gv1gv1*v
+2*et3*gv1gv1*gv1
)
b=(
+3*et*et2*gv3*gv1gv1*v
+3*et*gv3*gv1gv1*v
-4*et3*gv1gv1*gv1
-6*et*et2*gv3*gv1gv1
)
cc=(
+2*et3*gv1gv1*gv1
+6*et*et2*gv3*gv1gv1
+4*et^3*gv3*gvgv
-3*et*et2*gv3*gv1gv1*v
)
discriminant=b*b-4*a*cc
ans=ncp.mean=ncp.var=numeric(0)
if(discriminant>=0){ ## having real solutions
pi01=(sqrt(discriminant)-b)/2/a ### solution 1
mu=et/(1-pi01)/sqrt(v/2)/gv1*gv
(sigma22=(et2-pi01*v/(v-2))/(1-pi01)/v*(v-2)-1-mu*mu)## this is based on 2nd moment eqn.
(sigma23=(et3/(1-pi01)*2/v/sqrt(v/2)*gv/gv3-mu^3)/3/mu-1)## this is based on 3rd moment eqn. ## agrees:)
if(pi01>=0 && pi01<=1 && sigma22>0 && abs(sigma22-sigma23)<1e-6){
ans=c(ans,pi01)
ncp.mean=c(ncp.mean,mu)
ncp.var=c(ncp.var,(sigma22+sigma23)/2)
}
pi01=(-sqrt(discriminant)-b)/2/a ### solution 2
mu=et/(1-pi01)/sqrt(v/2)/gv1*gv
(sigma22=(et2-pi01*v/(v-2))/(1-pi01)/v*(v-2)-1-mu*mu)## this is based on 2nd moment eqn.
(sigma23=(et3/(1-pi01)*2/v/sqrt(v/2)*gv/gv3-mu^3)/3/mu-1)## this is based on 3rd moment eqn. ## agrees:)
if(pi01>=0 && pi01<=1 && sigma22>0 && abs(sigma22-sigma23)<1e-6){
ans=c(ans,pi01)
ncp.mean=c(ncp.mean,mu)
ncp.var=c(ncp.var,(sigma22+sigma23)/2)
}
if(length(ans)==1){ ## having 1 feasible solution
eff.mean=ncp.mean*sqrt(1/n1+1/n2)
gamma2=ncp.var*(1/n1+1/n2)
attr(ans,'eff.mean')=eff.mean
attr(ans,'gamma2')=gamma2
attr(ans,'method')='non-zero mean, exact solution'
return(ans)
# cat('unique pi0 ans=',ans,ncp.mean,sigma22,sigma23,fill=TRUE)
}else if (length(ans)==2){ ## having 2 feasible solutions
# cat('double pi0 ans=',ans,ncp.mean,ncp.var,fill=TRUE)
quants=quantile(tstat,seq(0,1,length=500))
deltas=seq(quantile(tstat,.001),quantile(tstat,.999),length=201)
dd=suppressWarnings(outer(quants,(deltas[-1]+deltas[-length(deltas)])/2,
function(xx,yy)dt(xx,v,yy)))
ll=c(-sum(log(
ans[1]*dt(quants,v)+(1-ans[1])*dd%*%(pnorm(deltas[-1],ncp.mean[1],sqrt(ncp.var[1]))-
pnorm(deltas[-length(deltas)],ncp.mean[1],sqrt(ncp.var[1])))
)),-sum(log(
ans[2]*dt(quants,v)+(1-ans[2])*dd%*%(pnorm(deltas[-1],ncp.mean[2],sqrt(ncp.var[2]))-
pnorm(deltas[-length(deltas)],ncp.mean[2],sqrt(ncp.var[2])))
)))
ans=ans[which.min(ll)]
eff.mean=ncp.mean[which.min(ll)]*sqrt(1/n1+1/n2)
gamma2=ncp.var[which.min(ll)]*(1/n1+1/n2)
attr(ans,'eff.mean')=eff.mean
attr(ans,'gamma2')=gamma2
attr(ans,'method')='non-zero mean, exact solution'
return(ans)
# cat('mle pi0 ans=',ans,ncp.mean[which.min(ll)],ncp.var[which.min(ll)],fill=TRUE)
}
}
## if program runs thru here, neither solutions are feasible
if(!approx)return(NA)
opt.fit=try(optimize(function(pp)abs(a*pp*pp+b*pp+cc),
c(0,lastbin(2*pt(abs(tstat),v,lower.tail=FALSE))), ### this range could be changed to c(0,1)
tol=1e-12))
if(inherits(opt.fit,'try-error')) return(NA)
# print(opt.fit)
pi01=opt.fit$minimum
mu=et/(1-pi01)/sqrt(v/2)/gv1*gv
(sigma22=(et2-pi01*v/(v-2))/(1-pi01)/v*(v-2)-1-mu*mu)## this is based on 2nd moment eqn.
(sigma23=(et3/(1-pi01)*2/v/sqrt(v/2)*gv/gv3-mu^3)/3/mu-1)## this is based on 3rd moment eqn.
## BECAUSE THIS IS APPROX SOLUTION, THE ABOVE TWO WILL NOT AGREE.
# cat(pi01,mu,sigma22,sigma23,fill=TRUE)
if(sigma22>0 && sigma23>0 ){## two solutions, check mle
ans=rep(pi01,2)
ncp.mean=rep(mu,2)
ncp.var=c(sigma22,sigma23)
# cat('double approx pi0 ans=',ans,ncp.mean,ncp.var,fill=TRUE)
quants=quantile(tstat,seq(0,1,length=500))
deltas=seq(quantile(tstat,.001),quantile(tstat,.999),length=201)
dd=suppressWarnings(outer(quants,(deltas[-1]+deltas[-length(deltas)])/2,
function(xx,yy)dt(xx,v,yy)))
ll=c(-sum(log(
ans[1]*dt(quants,v)+(1-ans[1])*dd%*%(pnorm(deltas[-1],ncp.mean[1],sqrt(ncp.var[1]))-
pnorm(deltas[-length(deltas)],ncp.mean[1],sqrt(ncp.var[1])))
)),-sum(log(
ans[2]*dt(quants,v)+(1-ans[2])*dd%*%(pnorm(deltas[-1],ncp.mean[2],sqrt(ncp.var[2]))-
pnorm(deltas[-length(deltas)],ncp.mean[2],sqrt(ncp.var[2])))
)))
ans=ans[which.min(ll)]
eff.mean=ncp.mean[which.min(ll)]*sqrt(1/n1+1/n2)
gamma2=ncp.var[which.min(ll)]*(1/n1+1/n2)
attr(ans,'eff.mean')=eff.mean
attr(ans,'gamma2')=gamma2
attr(ans,'method')='non-zero mean, approximate solution'
return(ans)
# cat('mle approx pi0 ans=',ans[1],ncp.mean[1],ncp.var[which.min(ll)],fill=TRUE)
}else if (sigma22*sigma23<0){ ## one solution
ans=pi01
ncp.mean=mu
ncp.var=max(c(sigma22,sigma23))
eff.mean=ncp.mean*sqrt(1/n1+1/n2)
gamma2=ncp.var*(1/n1+1/n2)
attr(ans,'eff.mean')=eff.mean
attr(ans,'gamma2')=gamma2
attr(ans,'method')='non-zero mean, approximate solution'
return(ans)
# cat('unique approx pi0 ans=',ans,ncp.mean,sigma22,sigma23,fill=TRUE)
}else { ## no solution, report result, but giving a warning
ans=pi01
ncp.mean=mu
ncp.var=mean(c(sigma22,sigma23))
warning('gamma2<0')
eff.mean=ncp.mean*sqrt(1/n1+1/n2)
gamma2=ncp.var*(1/n1+1/n2)
attr(ans,'eff.mean')=eff.mean
attr(ans,'gamma2')=gamma2
attr(ans,'method')='non-zero mean, approximate solution'
return(ans)
# cat('bad approx pi0 ans=',ans,ncp.mean,sigma22,sigma23,fill=TRUE)
}
}## of non-zero mean model
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.