Nothing
ScoreTest_wSaddleApprox_Get_X1 = function(X1)
{
q1<-ncol(X1)
if(q1>=2)
{
if(sum(abs(X1[,1]-X1[,2]))==0)
{
X1=X1[,-2]
q1<-q1-1
}
}
qr1<-qr(X1)
if(qr1$rank < q1){
X1.svd<-svd(X1)
X1 = X1.svd$u[,1:qr1$rank]
}
return(X1)
}
ScoreTest_wSaddleApprox_NULL_Model <- function(formula, data=NULL)
{
X1<-model.matrix(formula,data=data)
X1<-ScoreTest_wSaddleApprox_Get_X1(X1)
glmfit= glm(formula, data=data, family = "binomial")
convflag=0
if(glmfit$converged)
{
mu = glmfit$fitted.values
if(mean(mu)/mean(glmfit$y)>0.001 & (1-mean(mu))/(1-mean(glmfit$y))>0.001) convflag<-1 #Check that the null model converged properly with glm
}
if(convflag==0)
{
firthfit=fast.logistf.fit(x=X1,y=glmfit$y)
eta<-X1%*%firthfit$beta
mu = as.vector(exp(eta)/(1+exp(eta)))
}
V = mu*(1-mu)
res = glmfit$y- mu
n1<-length(res)
XV = t(X1 * V)
XVX_inv= solve(t(X1)%*%(X1 * V))
XXVX_inv= X1 %*% XVX_inv
re<-list(y=glmfit$y, mu=mu, res=res, V=V, X1=X1, XV=XV, XXVX_inv =XXVX_inv)
class(re)<-"SA_NULL"
return(re)
}
add_logp<-function(p1,p2)
{
p1<- -abs(p1)
p2<- -abs(p2)
maxp<-max(p1,p2)
minp<-min(p1,p2)
return(maxp+log(1+exp(minp-maxp)))
}
Korg<-function(t, mu, g)
{
n.t<-length(t)
out<-rep(0,n.t)
for(i in 1:n.t){
t1<-t[i]
temp<-log(1 - mu + mu * exp(g* t1))
out[i]<-sum(temp)
}
return(out)
}
K1_adj<-function(t, mu, g, q)
{
n.t<-length(t)
out<-rep(0,n.t)
for(i in 1:n.t){
t1<-t[i]
temp1<-(1 - mu)* exp(-g * t1) + mu
temp2<-mu *g
out[i]<-sum(temp2/temp1)-q
}
return(out)
}
K2<-function(t, mu, g)
{
n.t<-length(t)
out<-rep(0,n.t)
for(i in 1:n.t){
t1<-t[i]
temp1<-((1 - mu)* exp(-g * t1) + mu)^2
temp2<-(1-mu) * mu * g^2 * exp(-g*t1)
out[i]<-sum(temp2/temp1,na.rm=TRUE)
}
return(out)
}
Get_Normal_Spline<-function(var,nodes)
{
y1<-var*nodes
y2<-rep(var,length(nodes))
return(cbind(y1,y2))
}
Get_Saddle_Spline<-function(mu,g,nodes)
{
m1<-sum(mu*g)
y1<-K1_adj(nodes,mu,g,0)-m1
y2<-K2(nodes,mu,g)
return(cbind(nodes,y1,y2)[which(nodes!=0),])
}
getroot_K1<-function(init,mu,g,q,m1,tol=.Machine$double.eps^0.25,maxiter=1000)
{
g.pos<-sum(g[which(g>0)])
g.neg<- sum(g[which(g<0)])
if(q>=g.pos || q<=g.neg)
{
return(list(root=Inf,n.iter=0,Is.converge=TRUE))
} else {
t<-init
K1_eval<-K1_adj(t,mu,g,q)
prevJump<- Inf
rep<-1
repeat
{
K2_eval<-K2(t,mu,g)
tnew<-t-K1_eval/K2_eval
if(is.na(tnew))
{
conv=FALSE
break
}
if(abs(tnew-t)<tol)
{
conv<-TRUE
break
}
if(rep==maxiter)
{
conv<-FALSE
break
}
newK1<-K1_adj(tnew,mu,g,q)
if(sign(K1_eval)!=sign(newK1))
{
if(abs(tnew-t)>prevJump-tol)
{
tnew<-t+sign(newK1-K1_eval)*prevJump/2
newK1<-K1_adj(tnew,mu,g,q)
prevJump<-prevJump/2
} else {
prevJump<-abs(tnew-t)
}
}
rep<-rep+1
t<-tnew
K1_eval<-newK1
}
return(list(root=t,n.iter=rep,Is.converge=conv))
}
}
Get_Saddle_Prob<-function(zeta, mu, g, q,log.p=FALSE)
{
k1<-Korg(zeta, mu, g)
k2<-K2(zeta, mu, g)
if(is.finite(k1) && is.finite(k2))
{
temp1<-zeta * q - k1
w<-sign(zeta) * (2 *temp1)^{1/2}
v<- zeta * (k2)^{1/2}
Z.test<-w + 1/w * log(v/w)
if(Z.test > 0){
pval<-pnorm(Z.test, lower.tail = FALSE,log.p=log.p)
} else {
pval= -pnorm(Z.test, lower.tail = TRUE,log.p=log.p)
}
} else {
if(log.p)
{
pval<- -Inf
} else {
pval<-0
}
}
return(pval)
}
Saddle_Prob<-function(q, mu, g, Cutoff=2,alpha,output="P",nodes.fixed,nodes.init,log.p=FALSE)
{
m1<-sum(mu * g)
var1<-sum(mu * (1-mu) * g^2)
p1=NULL
p2=NULL
if(output=="metaspline")
{
if(is.null(nodes.fixed)==TRUE)
{
nodes<-nodes.init
nodes<-getnodes(nodes,mu,g)$nodes
} else {
nodes<-unique(c(nodes.fixed,0))
nodes<-nodes[order(nodes)]
}
}
Score<- q-m1
#
qinv = -sign(q-m1) * abs(q-m1) + m1
# Noadj
pval.noadj<-pchisq((q - m1)^2/var1, lower.tail = FALSE, df=1,log.p=log.p)
Is.converge=TRUE
if(Cutoff=="BE"){
rho<-sum(((abs(g))^3)*mu*(1-mu)*(mu^2+(1-mu)^2))
B<-0.56*rho*var1^(-3/2)
p<-B+alpha/2
Cutoff=ifelse(p>=0.496,0.01,qnorm(p,lower.tail=F))
} else if(Cutoff < 10^-1){
Cutoff=10^-1
}
#
if(output=="metaspline")
{
splfun<-Get_Saddle_Spline(mu,g,nodes)
}
if(abs(q - m1)/sqrt(var1) < Cutoff ){
pval=pval.noadj
} else {
out.uni1<-getroot_K1(0, mu=mu, g=g, q=q)
out.uni2<-getroot_K1(0, mu=mu, g=g, q=qinv)
if(out.uni1$Is.converge==TRUE && out.uni2$Is.converge==TRUE)
{
p1<-tryCatch(Get_Saddle_Prob(out.uni1$root, mu, g, q,log.p=log.p),error=function(e) {
if(log.p) return(pval.noadj-log(2))
else return(pval.noadj/2)})
p2<-tryCatch(Get_Saddle_Prob(out.uni2$root, mu, g, qinv,log.p=log.p),error=function(e) {
if(log.p) return(pval.noadj-log(2))
else return(pval.noadj/2)})
if(log.p)
{
pval = add_logp(p1,p2)
} else {
pval = abs(p1)+abs(p2)
}
Is.converge=TRUE
} else {
print("Error_Converge")
pval<-pval.noadj
Is.converge=FALSE
}
}
if(pval!=0 && pval.noadj/pval>10^3)
{
return(Saddle_Prob(q, mu, g, Cutoff=Cutoff*2,alpha,output,nodes.fixed,nodes.init,log.p=log.p))
} else if(output=="metaspline")
{
return(list(p.value=pval, p.value.NA=pval.noadj, Is.converge=Is.converge,Score=Score,splfun=splfun,var=var1))
} else {
return(list(p.value=pval, p.value.NA=pval.noadj, Is.converge=Is.converge, Score=Score))
}
}
TestSPA<-function(G, obj.null, Cutoff=2,alpha,output,nodes.fixed=NULL,nodes.init,log.p=FALSE)
{
if(class(obj.null) != "SA_NULL"){
stop("obj.null should be a returned object from ScoreTest_wSaddleApprox_NULL_Model")
}
y = obj.null$y
mu = obj.null$mu
res = obj.null$res
n.g<-sum(G)
if(n.g/(2*length(G))>0.5)
{
G<-2-G
n.g<-sum(G)
}
G1<-G - obj.null$XXVX_inv %*% (obj.null$XV %*% G)
q<-sum(G1 * y)
out<-Saddle_Prob(q, mu=mu, g=G1, Cutoff=Cutoff,alpha=alpha,output=output,nodes.fixed=nodes.fixed,nodes.init=nodes.init,log.p=log.p)
return(out)
}
Korg_fast<-function(t, mu, g,gNA,gNB,muNA,muNB,NAmu,NAsigma)
{
n.t<-length(t)
out<-rep(0,n.t)
for(i in 1:n.t){
t1<-t[i]
temp<-log(1 - muNB + muNB * exp(gNB* t1))
out[i]<-sum(temp)+NAmu*t1+0.5*NAsigma*t1^2
}
return(out)
}
K1_adj_fast<-function(t, mu, g, q,gNA,gNB,muNA,muNB,NAmu,NAsigma)
{
n.t<-length(t)
out<-rep(0,n.t)
for(i in 1:n.t){
t1<-t[i]
temp1<-(1 - muNB)* exp(-gNB * t1) + muNB
temp2<-muNB *gNB
temp3<-NAmu+NAsigma*t1
out[i]<-sum(temp2/temp1)+temp3-q
}
return(out)
}
K2_fast<-function(t, mu, g,gNA,gNB,muNA,muNB,NAmu,NAsigma)
{
n.t<-length(t)
out<-rep(0,n.t)
for(i in 1:n.t){
t1<-t[i]
temp1<-((1 - muNB)* exp(-gNB * t1) + muNB)^2
temp2<-(1-muNB) * muNB * gNB^2 * exp(-gNB*t1)
out[i]<-sum(temp2/temp1,na.rm=TRUE)+NAsigma
}
return(out)
}
Get_Saddle_Spline_fast<-function(nodes,mu,g,gNA,gNB,muNA,muNB,NAmu,NAsigma,m1)
{
y1<-K1_adj_fast(nodes,mu, g, 0,gNA,gNB,muNA,muNB,NAmu,NAsigma)-m1
y2<-K2_fast(nodes,mu, g,gNA,gNB,muNA,muNB,NAmu,NAsigma)
return(cbind(nodes,y1,y2)[which(nodes!=0),])
}
getroot_K1_fast<-function(init,mu,g,q,m1,gNA,gNB,muNA,muNB,NAmu,NAsigma,tol=.Machine$double.eps^0.25,maxiter=1000)
{
g.pos<-sum(g[which(g>0)])
g.neg<- sum(g[which(g<0)])
if(q>=g.pos || q<=g.neg)
{
return(list(root=Inf,n.iter=0,Is.converge=TRUE))
} else {
t<-init
K1_eval<-K1_adj_fast(t,mu,g,q,gNA,gNB,muNA,muNB,NAmu,NAsigma)
prevJump<- Inf
rep<-1
repeat
{
K2_eval<-K2_fast(t,mu,g,gNA,gNB,muNA,muNB,NAmu,NAsigma)
tnew<-t-K1_eval/K2_eval
if(is.na(tnew))
{
conv=FALSE
break
}
if(abs(tnew-t)<tol)
{
conv<-TRUE
break
}
if(rep==maxiter)
{
conv<-FALSE
break
}
newK1<-K1_adj_fast(tnew,mu,g,q,gNA,gNB,muNA,muNB,NAmu,NAsigma)
if(sign(K1_eval)!=sign(newK1))
{
if(abs(tnew-t)>prevJump-tol)
{
tnew<-t+sign(newK1-K1_eval)*prevJump/2
newK1<-K1_adj_fast(tnew,mu,g,q,gNA,gNB,muNA,muNB,NAmu,NAsigma)
prevJump<-prevJump/2
} else {
prevJump<-abs(tnew-t)
}
}
rep<-rep+1
t<-tnew
K1_eval<-newK1
}
return(list(root=t,n.iter=rep,Is.converge=conv))
}
}
Get_Saddle_Prob_fast<-function(zeta, mu, g, q,gNA,gNB,muNA,muNB,NAmu,NAsigma,log.p=FALSE)
{
k1<-Korg_fast(zeta, mu, g,gNA,gNB,muNA,muNB,NAmu,NAsigma)
k2<-K2_fast(zeta, mu, g,gNA,gNB,muNA,muNB,NAmu,NAsigma)
if(is.finite(k1) && is.finite(k2))
{
temp1<-zeta * q - k1
w<-sign(zeta) * (2 *temp1)^{1/2}
v<- zeta * (k2)^{1/2}
Z.test<-w + 1/w * log(v/w)
if(Z.test > 0){
pval<-pnorm(Z.test, lower.tail = FALSE,log.p=log.p)
} else {
pval= -pnorm(Z.test, lower.tail = TRUE,log.p=log.p)
}
} else {
if(log.p)
{
pval<- -Inf
} else {
pval<-0
}
}
return(pval)
}
Saddle_Prob_fast<-function(q, g,mu,gNA,gNB,muNA,muNB,Cutoff=2,alpha,output,nodes.fixed,nodes.init,log.p=FALSE)
{
m1<-sum(mu * g)
var1<-sum(mu * (1-mu) * g^2)
p1=NULL
p2=NULL
NAmu= m1-sum(gNB*muNB)
NAsigma=var1-sum(muNB*(1-muNB)*gNB^2)
if(output=="metaspline")
{
if(is.null(nodes.fixed)==TRUE)
{
nodes<-nodes.init
nodes<-getnodes_fast(nodes,mu,g,gNA,gNB,muNA,muNB,NAmu,NAsigma)$nodes
} else {
nodes<-unique(c(nodes.fixed,0))
nodes<-nodes[order(nodes)]
}
}
Score<- q-m1
qinv = -sign(q-m1) * abs(q-m1) + m1
# Noadj
pval.noadj<-pchisq((q - m1)^2/var1, lower.tail = FALSE, df=1,log.p=log.p)
Is.converge=TRUE
if(Cutoff=="BE"){
rho<-sum(((abs(g))^3)*mu*(1-mu)*(mu^2+(1-mu)^2))
B<-0.56*rho*var1^(-3/2)
p<-B+alpha/2
Cutoff=ifelse(p>=0.496,0.01,qnorm(p,lower.tail=F))
} else if(Cutoff < 10^-1){
Cutoff=10^-1
}
#
if(output=="metaspline")
{
splfun<-Get_Saddle_Spline_fast(nodes,mu,g,gNA,gNB,muNA,muNB,NAmu,NAsigma,m1)
}
if(abs(q - m1)/sqrt(var1) < Cutoff ){
pval=pval.noadj
} else {
out.uni1<-getroot_K1_fast(0, mu=mu, g=g, q=q,gNA=gNA,gNB=gNB,muNA=muNA,muNB=muNB,NAmu=NAmu,NAsigma=NAsigma)
out.uni2<-getroot_K1_fast(0, mu=mu, g=g, q=qinv,gNA=gNA,gNB=gNB,muNA=muNA,muNB=muNB,NAmu=NAmu,NAsigma=NAsigma)
if(out.uni1$Is.converge==TRUE && out.uni2$Is.converge==TRUE)
{
p1<-tryCatch(Get_Saddle_Prob_fast(out.uni1$root, mu, g, q,gNA,gNB,muNA,muNB,NAmu,NAsigma,log.p=log.p),error=function(e) {
if(log.p) return(pval.noadj-log(2))
else return(pval.noadj/2)})
p2<-tryCatch(Get_Saddle_Prob_fast(out.uni2$root, mu, g, qinv,gNA,gNB,muNA,muNB,NAmu,NAsigma,log.p=log.p),error=function(e) {
if(log.p) return(pval.noadj-log(2))
else return(pval.noadj/2)})
if(log.p)
{
pval = add_logp(p1,p2)
} else {
pval = abs(p1)+abs(p2)
}
Is.converge=TRUE
} else {
print("Error_Converge")
pval<-pval.noadj
Is.converge=FALSE
}
}
if(pval!=0 && pval.noadj/pval>10^3)
{
return(Saddle_Prob_fast(q, g,mu,gNA,gNB,muNA,muNB,Cutoff=Cutoff*2,alpha,output,nodes.fixed,nodes.init,log.p=log.p))
} else if(output=="metaspline")
{
return(list(p.value=pval, p.value.NA=pval.noadj, Is.converge=Is.converge, Score=Score,splfun=splfun,var=var1))
} else {
return(list(p.value=pval, p.value.NA=pval.noadj, Is.converge=Is.converge, Score=Score))
}
}
TestSPAfast<-function(G, obj.null, output, Cutoff=2,alpha,nodes.fixed,nodes.init,log.p=FALSE)
{
if(class(obj.null) != "SA_NULL"){
stop("obj.null should be a returned object from ScoreTest_wSaddleApprox_NULL_Model")
}
y = obj.null$y
mu = obj.null$mu
res = obj.null$res
n.g<-sum(G)
if(n.g/(2*length(G))>0.5)
{
G<-2-G
n.g<-sum(G)
}
NAset<-which(G==0)
G1<-G - obj.null$XXVX_inv %*% (obj.null$XV %*% G)
q<-sum(G1 * y)
g=G1
if(length(NAset)/length(G)<0.5)
{
out<-Saddle_Prob(q, mu=mu, g=g, Cutoff=Cutoff,alpha=alpha, output=output,nodes.fixed=nodes.fixed,nodes.init=nodes.init,log.p=log.p)
} else {
out<-Saddle_Prob_fast(q,g=g,mu=mu,gNA=g[NAset],gNB=g[-NAset],
muNA=mu[NAset],muNB=mu[-NAset], Cutoff=Cutoff,alpha=alpha,output=output,nodes.fixed=nodes.fixed,nodes.init=nodes.init,log.p=log.p)
}
return(out)
}
fast.logistf.fit <- function (x, y, weight = NULL, offset = NULL, firth = TRUE, col.fit = NULL,
init = NULL, control) {
n <- nrow(x)
k <- ncol(x)
if (is.null(init))
init = rep(0, k)
if (is.null(col.fit))
col.fit = 1:k
if (is.null(offset))
offset = rep(0, n)
if (is.null(weight))
weight = rep(1, n)
if (col.fit[1] == 0)
maxit <- 0
if (missing(control))
control <- fast.logistf.control()
maxit <- control$maxit
maxstep <- control$maxstep
maxhs <- control$maxhs
lconv <- control$lconv
gconv <- control$gconv
xconv <- control$xconv
beta <- init
iter <- 0
pi <- as.vector(1/(1 + exp(-x %*% beta - offset)))
evals <- 1
repeat {
beta.old <- beta
XW2 <- t(x * (weight * pi * (1-pi))^0.5)
myQR <- qr(t(XW2))
Q <- qr.Q(myQR)
h <- (Q*Q) %*% rep(1, ncol(Q))
if (firth)
U.star <- crossprod(x, weight * (y - pi) + h * (0.5 - pi))
else U.star <- crossprod(x, weight * (y - pi))
XX.covs <- matrix(0, k, k)
if (col.fit[1] != 0) {
XX.XW2 <- t(x[, col.fit, drop=FALSE] * (weight * pi * (1-pi))^0.5)
XX.Fisher <- crossprod(t(XX.XW2))
XX.covs[col.fit, col.fit] <- fast.invFisher(XX.Fisher) ###### HERE IS THE PROBLEM!!!
}
if(all(is.na(XX.covs)) == T) {
break
}
delta <- as.vector(XX.covs %*% U.star)
delta[is.na(delta)] <- 0
mx <- max(abs(delta))/maxstep
if (mx > 1)
delta <- delta/mx
evals <- evals + 1
if (maxit > 0) {
iter <- iter + 1
beta <- beta + delta
pi <- as.vector(1/(1 + exp(-x %*% beta - offset)))
}
if (iter == maxit | ((max(abs(delta)) <= xconv) & (all(abs(U.star[col.fit]) <
gconv))))
break
}
# Error catching (if chol(x) not positive definite)
if(all(is.na(XX.covs))==T) {
var <- XX.covs
list(beta = NA, var = var, pi = NA, hat.diag = NA,
iter = NA, evals = NA, conv = c(NA,
NA, NA))
} else {
var <- XX.covs
list(beta = beta, var = var, pi = pi, hat.diag = h,
iter = iter, evals = evals, conv = c(max(abs(U.star)),
max(abs(delta))))
}
}
fast.logistf.control <- function (maxit = 50, maxhs = 15, maxstep = 15, lconv = 1e-05,
gconv = 1e-05, xconv = 1e-05)
{
list(maxit = maxit, maxhs = maxhs, maxstep = maxstep, lconv = lconv,
gconv = gconv, xconv = xconv)
}
fast.logDet <- function (x) {
my.chol <- tryCatch(chol(x),error=function(e) {NA})
if (all(is.na(my.chol))==T) {
return(NA)
} else {
return (2 * sum(log(diag(my.chol))))
}
}
fast.invFisher <- function(x) {
my.chol <- tryCatch(chol(x),error=function(e) {NA})
if (all(is.na(my.chol))==T) {
return(NA)
} else {
return (chol2inv(my.chol))
}
#ifelse(is.na(my.chol), NA, chol2inv(my.chol))
}
getnodes<-function(init,mu,g)
{
nodes<-unique(c(init,0))
nodes<-nodes[order(nodes)]
flag=0
rep<-0
t<-c(2^seq(-2,10,0.5),-2^seq(-2,10,0.5),0)
t<-t[order(t)]
yt<-K1_adj(t,mu,g,0)
w<-rep(1,length(t))
w[which(t^2>1)]=1/(abs(t[which(t^2>1)]))^(1/3)
jump<-0.5
totres<-Inf
repeat{
y1<-K1_adj(nodes ,mu ,g,0)
y2<-K2(nodes ,mu ,g)
sfun<-splinefunH(nodes, y1 ,m=y2)
pred<-sfun(t,deriv=0)
res<-w*abs(pred-yt)
newnode=NULL
for(i in 1:length(nodes))
{
if(nodes[i]!=0)
{
if(i>1)
{
int1<-which(t<nodes[i] & t>nodes[i-1])
} else {
int1<-which(t<nodes[i])
}
if(i<length(nodes))
{
int2<-which(t<nodes[i+1] & t>nodes[i])
} else {
int2<-which(t>nodes[i])
}
if(length(int1)>0)
{
r1=res[int1]
t1=t[int1]
t1<-t1[min(which(r1==max(r1)))]
r1<-max(r1)
} else {
r1<-0
}
if(length(int2)>0)
{
r2=res[int2]
t2=t[int2]
t2<-t2[min(which(r2==max(r2)))]
r2<-max(r2)
} else {
r2<-0
}
if(r1==r2)
{
newnode<-c(newnode,nodes[i])
} else if(r1>r2)
{
jump1<- jump*(1-max(0.1,r2/r1))*abs(nodes[i])
newnode<-c(newnode,nodes[i]-jump1)
} else {
jump1<- jump*(1-max(0.1,r1/r2))*abs(nodes[i])
newnode<-c(newnode,nodes[i]+jump1)
}
}
}
newnode<-c(newnode,0)
if(rep>100 || (sum(res)>totres)) break
rep<-rep+1
totres<-sum(res)
nodes<-newnode[order(newnode)]
}
return(list(nodes=nodes,loss=totres))
}
getnodes_fast<-function(init,mu,g,gNA,gNB,muNA,muNB,NAmu,NAsigma)
{
nodes<-unique(c(init,0))
nodes<-nodes[order(nodes)]
flag=0
rep<-0
t<-c(2^seq(-2,10,0.5),-2^seq(-2,10,0.5),0)
t<-t[order(t)]
yt<-K1_adj_fast(t,mu, g, 0,gNA,gNB,muNA,muNB,NAmu,NAsigma)
w<-rep(1,length(t))
w[which(t^2>1)]=1/(abs(t[which(t^2>1)]))^(1/3)
jump<-0.5
totres<-Inf
repeat{
y1<-K1_adj(nodes ,mu ,g,0)
y2<-K2(nodes ,mu ,g)
sfun<-splinefunH(nodes, y1 ,m=y2)
pred<-sfun(t,deriv=0)
res<-w*abs(pred-yt)
newnode=NULL
for(i in 1:length(nodes))
{
if(nodes[i]!=0)
{
if(i>1)
{
int1<-which(t<nodes[i] & t>nodes[i-1])
} else {
int1<-which(t<nodes[i])
}
if(i<length(nodes))
{
int2<-which(t<nodes[i+1] & t>nodes[i])
} else {
int2<-which(t>nodes[i])
}
if(length(int1)>0)
{
r1=res[int1]
t1=t[int1]
t1<-t1[min(which(r1==max(r1)))]
r1<-max(r1)
} else {
r1<-0
}
if(length(int2)>0)
{
r2=res[int2]
t2=t[int2]
t2<-t2[min(which(r2==max(r2)))]
r2<-max(r2)
} else {
r2<-0
}
if(r1==r2)
{
newnode<-c(newnode,nodes[i])
} else if(r1>r2)
{
jump1<- jump*(1-max(0.1,r2/r1))*abs(nodes[i])
newnode<-c(newnode,nodes[i]-jump1)
} else {
jump1<- jump*(1-max(0.1,r1/r2))*abs(nodes[i])
newnode<-c(newnode,nodes[i]+jump1)
}
}
}
newnode<-c(newnode,0)
if(rep>100 || (sum(res)>totres)) break
rep<-rep+1
totres<-sum(res)
nodes<-newnode[order(newnode)]
}
return(list(nodes=nodes,loss=totres))
}
ScoreTest_SPA <-function(genos,pheno,cov,obj.null,method=c("fastSPA","SPA"),minmac=5,Cutoff=2,alpha=5*10^-8,missing.id=NA,beta.out=FALSE,beta.Cutoff=5*10^-7,log.p=FALSE)
{
method<-match.arg(method)
if(missing(obj.null))
{
if(missing(cov) || is.null(cov))
{
cov<-rep(1,length(pheno))
}
obj.null<-ScoreTest_wSaddleApprox_NULL_Model(as.matrix(pheno) ~as.matrix(cov))
}
cov<-obj.null$X1
pheno<-obj.null$y
genos<-as.matrix(genos)
if(ncol(genos)==1)
{
m<-1
genos<-t(genos)
} else {
m <- nrow(genos)
}
if(is.na(missing.id)==FALSE)
{
genos[which(genos==missing.id)]<-NA
}
p.value<-rep(NA,m)
p.value.NA<-rep(NA,m)
Is.converge<-rep(NA,m)
beta<-rep(NA,m)
SEbeta<-rep(NA,m)
for (i in 1:m)
{
try({
ina<-which(is.na(genos[i,]))
if(length(ina)>0)
{
genos[i,ina]<-mean(genos[i,],na.rm=TRUE)
}
MAC<-min(sum(genos[i,]),sum(2-genos[i,]))
if(MAC>=minmac)
{
if(method=="fastSPA")
{
re <- TestSPAfast(as.vector(genos[i,,drop=FALSE]), obj.null, Cutoff=Cutoff, alpha=alpha, output="P", log.p=log.p)
} else {
re <- TestSPA(as.vector(genos[i,,drop=FALSE]), obj.null,Cutoff=Cutoff, alpha=alpha, output="P",log.p=log.p)
}
p.value[i] <- re$p.value
p.value.NA[i] <- re$p.value.NA
Is.converge[i]<- re$Is.converge
if(beta.out==TRUE && p.value[i]<beta.Cutoff)
{
re.firth<-fast.logistf.fit(x=cbind(t(genos[i,,drop=FALSE]),cov),y=pheno,firth=TRUE)
beta[i]<-re.firth$beta[1]
SEbeta[i]<-sqrt(re.firth$var[1,1])
}
}
})
if(i%%1000==0) print(paste("Processed",i,"SNPs",sep=" "))
}
return(list(p.value=p.value,p.value.NA=p.value.NA,Is.converge=Is.converge,beta=beta,SEbeta=SEbeta))
}
ScoreTest_SPA_wMeta <-function(genos,pheno,cov,obj.null,method=c("fastSPA","SPA"),minmac=5,Cutoff=2,alpha=5*10^-8,missing.id=NA,beta.out=FALSE,beta.Cutoff=5*10^-7,output=c("P","metaZ","metaGC","metaspline"),nodes.fixed=NULL,nodes.init=c(-100,-10,-1,1,10,100))
{
method<-match.arg(method)
output<-match.arg(output)
if(missing(obj.null))
{
obj.null<-ScoreTest_wSaddleApprox_NULL_Model(as.matrix(pheno) ~as.matrix(cov))
}
cov<-obj.null$X1
pheno<-obj.null$y
if(is.null(nodes.fixed))
{
nodes.init<-unique(nodes.init[which(nodes.init!=0)])
} else {
nodes.init<-unique(nodes.fixed[which(nodes.fixed!=0)])
}
genos<-as.matrix(genos)
if(ncol(genos)==1)
{
m<-1
genos<-t(genos)
} else {
m <- nrow(genos)
}
if(is.na(missing.id)==FALSE)
{
genos[which(genos==missing.id)]<-NA
}
p.value<-rep(NA,m)
p.value.NA<-rep(NA,m)
GCmat<-matrix(NA,m,2)
Is.converge<-rep(NA,m)
beta<-rep(NA,m)
SEbeta<-rep(NA,m)
Score<-rep(NA,m)
#Skew<-rep(NA,m)
nodes<-matrix(NA,m,length(nodes.init))
splfuny1<-matrix(NA,m,length(nodes.init))
splfuny2<-matrix(NA,m,length(nodes.init))
Var<-rep(NA,m)
MAF<-rep(NA,m)
for (i in 1:m)
{
try({
ina<-which(is.na(genos[i,]))
if(length(ina)>0)
{
genos[i,ina]<-mean(genos[i,],na.rm=TRUE)
}
AC1<-sum(genos[i,])
AC2<-sum(2-genos[i,])
MAC<-min(AC1,AC2)
if(MAC>=minmac)
{
if(method=="fastSPA")
{
re <- TestSPAfast(as.vector(genos[i,,drop=FALSE]), obj.null, Cutoff=Cutoff, alpha=alpha,output=output, nodes.fixed=nodes.fixed,nodes.init=nodes.init)
} else {
re <- TestSPA(as.vector(genos[i,,drop=FALSE]), obj.null,Cutoff=Cutoff, alpha=alpha,output=output, nodes.fixed=nodes.fixed,nodes.init=nodes.init)
}
p.value[i] <- re$p.value
p.value.NA[i] <- re$p.value.NA
Is.converge[i]<- re$Is.converge
if(beta.out==TRUE && p.value[i]<beta.Cutoff)
{
re.firth<-fast.logistf.fit(x=cbind(t(genos[i,,drop=FALSE]),cov),y=pheno,firth=TRUE)
beta[i]<-re.firth$beta[1]
SEbeta[i]<-sqrt(re.firth$var[1,1])
}
if(output!="P")
{
GCmat[i,]<-c(length(which(round(genos[i,],0)==2)),length(which(round(genos[i,],0)==1)))
if(AC1>AC2) GCmat[i,]<-c(length(which(round(genos[i,],0)==0)),length(which(round(genos[i,],0)==1)))
Score[i]<-re$Score
MAF[i]<-MAC/(2*length(pheno))
if(output=="metaspline")
{
nodes[i,1:nrow(re$splfun)]<-re$splfun[,1]
splfuny1[i,1:nrow(re$splfun)]<-re$splfun[,2]
splfuny2[i,1:nrow(re$splfun)]<-re$splfun[,3]
Var[i]<-re$var
}
}
}
})
if(i%%1000==0) print(paste("Processed",i,"SNPs",sep=" "))
}
if(output=="P")
{
return(list(p.value=p.value,p.value.NA=p.value.NA,Is.converge=Is.converge,beta=beta,SEbeta=SEbeta))
} else if(output=="metaZ"){
return(list(p.value=sign(Score)*p.value,p.value.NA=sign(Score)*p.value.NA,Is.converge=Is.converge,beta=beta,SEbeta=SEbeta,MAF=MAF))
} else if(output=="metaGC"){
return(list(p.value=sign(Score)*p.value,p.value.NA=sign(Score)*p.value.NA,Is.converge=Is.converge,beta=beta,SEbeta=SEbeta,GCmat=GCmat))
} else {
return(list(p.value=sign(Score)*p.value,p.value.NA=sign(Score)*p.value.NA,Is.converge=Is.converge,beta=beta,SEbeta=SEbeta,GCmat=GCmat,spldata=cbind(Score,Var,nodes,splfuny1,splfuny2)))
}
}
pval_MAC1<-function(q,p,m1=m1,var1=var1,GenC,gcent,truep,Cutoff=2)
{
n0<-GenC[1]
n1<-GenC[2]
n2<-GenC[3]
re<-Saddle_Prob_MAC(q, p=p,m1=m1,var1=var1, n0=n0,n1=n1,n2=n2,g=gcent, Cutoff=Cutoff)
return(re$p.value-truep)
}
Korg_MAC<-function(t, p, n0,n1,n2,g)
{
n.t<-length(t)
out<-rep(0,n.t)
for(i in 1:n.t){
t1<-t[i]
out[i]<-n1*log(1 - p + p * exp(g[2]*t1))+n2*log(1 - p + p * exp(g[3]*t1))+
n0*log(1 - p + p * exp(g[1]*t1))
}
return(out)
}
K1_adj_MAC<-function(t, p, n0,n1,n2,g, q)
{
n.t<-length(t)
out<-rep(0,n.t)
for(i in 1:n.t){
t1<-t[i]
temp1<-n1*g[2]*p/((1 - p)* exp(-g[2]*t1) + p)+n2*g[3]*p/((1 - p)* exp(-g[3]*t1) + p)+
n0*g[1]*p/((1 - p)* exp(-g[1]*t1) + p)
out[i]<-temp1-q
}
return(out)
}
K2_MAC<-function(t, p,n0,n1,n2,g)
{
n.t<-length(t)
out<-rep(0,n.t)
for(i in 1:n.t){
t1<-t[i]
temp<-rep(NA,3)
temp[2]<-n1*(1-p) * g[2]^2* p * exp(-g[2]*t1)/(((1 - p)* exp(-g[2]*t1) + p)^2)
temp[3]<-n2*(1-p) * g[3]^2* p * exp(-g[3]*t1)/(((1 - p)* exp(-g[3]*t1) + p)^2)
temp[1]<-n0*(1-p) * g[1]^2* p * exp(-g[1]*t1)/(((1 - p)* exp(-g[1]*t1) + p)^2)
out[i]<-sum(temp,na.rm=TRUE)
}
return(out)
}
getroot_K1_MAC<-function(init,p,n0,n1,n2,g,q,m1,tol=.Machine$double.eps^0.25,maxiter=1000)
{
gtot<-c(n0,n1,n2)*g
g.pos<-sum(gtot[which(gtot>0)])
g.neg<- -g.pos
if(q>=g.pos || q<=g.neg)
{
return(list(root=Inf,n.iter=0,Is.converge=TRUE))
} else {
t<-init
K1_eval<-K1_adj_MAC(t,p,n0,n1,n2,g,q)
prevJump<- Inf
rep<-1
repeat
{
K2_eval<-K2_MAC(t,p,n0,n1,n2,g)
tnew<-t-K1_eval/K2_eval
if(is.na(tnew))
{
conv=FALSE
break
}
if(abs(tnew-t)<tol)
{
conv<-TRUE
break
}
if(rep==maxiter)
{
conv<-FALSE
break
}
newK1<-K1_adj_MAC(tnew,p,n0,n1,n2,g,q)
if(sign(K1_eval)!=sign(newK1))
{
if(abs(tnew-t)>prevJump-tol)
{
tnew<-t+sign(newK1-K1_eval)*prevJump/2
newK1<-K1_adj_MAC(tnew,p,n0,n1,n2,g,q)
prevJump<-prevJump/2
} else {
prevJump<-abs(tnew-t)
}
}
rep<-rep+1
t<-tnew
K1_eval<-newK1
}
return(list(root=t,n.iter=rep,Is.converge=conv))
}
}
Get_Saddle_Prob_MAC<-function(zeta, p,n0,n1,n2,g, q)
{
k1<-Korg_MAC(zeta, p,n0,n1,n2,g)
k2<-K2_MAC(zeta, p,n0,n1,n2,g)
if(is.finite(k1) && is.finite(k2))
{
temp1<-zeta * q - k1
w<-sign(zeta) * (2 *temp1)^{1/2}
v<- zeta * (k2)^{1/2}
Z.test<-w + 1/w * log(v/w)
if(Z.test > 0){
pval<-pnorm(Z.test, lower.tail = FALSE)
} else {
pval= -pnorm(Z.test, lower.tail = TRUE)
}
} else {
pval<-0
}
return(pval)
}
Saddle_Prob_MAC<-function(q, p, m1,var1,n0,n1,n2,g, Cutoff=2)
{
qinv = -sign(q-m1) * abs(q-m1) + m1
# Noadj
pval.noadj<-pchisq((q - m1)^2/var1, lower.tail = FALSE, df=1)
p1<-sign(q-m1)*pval.noadj/2
p2= -sign(q-m1)*pval.noadj/2
Is.converge=TRUE
if(Cutoff < 10^-2){
Cutoff=10^-2
}
#
if(abs(q - m1)/sqrt(var1) < Cutoff ){
pval=pval.noadj
} else {
out.uni1<-getroot_K1_MAC(0, p=p,n0=n0,n1=n1,n2=n2, g=g,q=q)
out.uni2<-getroot_K1_MAC(0, p=p,n0=n0,n1=n1,n2=n2, g=g,q=qinv)
if(out.uni1$Is.converge==TRUE && out.uni2$Is.converge==TRUE)
{
p1<-tryCatch(Get_Saddle_Prob_MAC(out.uni1$root, p,n0,n1,n2,g, q),error=function(e) {return(pval.noadj/2)})
p2<-tryCatch(Get_Saddle_Prob_MAC(out.uni2$root, p,n0,n1,n2,g, qinv),error=function(e) {return(pval.noadj/2)})
pval = abs(p1)+abs(p2)
Is.converge=TRUE
} else {
print("Error_Converge")
pval<-pval.noadj
Is.converge=FALSE
}
}
if(pval!=0 && pval.noadj/pval>10^3)
{
return(Saddle_Prob_MAC(q, p,m1,var1,n0,n1,n2,g, Cutoff=Cutoff*2))
} else {
return(list(p.value=pval, p.value.NA=pval.noadj, Is.converge=Is.converge, p1=p1, p2=p2))
}
}
SPAmeta<-function(pvalue.Z=NULL,MAF.Z=NULL,CCsize.Z=NULL,
pvalue.GC=NULL,GCmat=NULL,CCsize.GC=NULL,Cutoff.GC=2,
spldata=NULL,
Cutoff.meta=2)
{
#Check all the inputs are correct
info<-c(length(pvalue.Z),length(pvalue.GC))
if(info[1]>0)
{
if(info[2]==0 & is.null(spldata))
{
return(MetaSPA_Normal(pvalue=pvalue.Z,MAF=MAF.Z,CCsize=CCsize.Z))
} else {
if(is.null(MAF.Z)) stop("MAF is needed when hybridizing Z-score based meta-analysis with other methods")
if(length(which(is.na(MAF.Z)))>0 | length(MAF.Z)!=length(pvalue.Z)) stop("MAF is needed when hybridizing Z-score based meta-analysis with other methods")
re<-Hybrid_Normal(pvalue=pvalue.Z,MAF=MAF.Z,CCsize=CCsize.Z)
q_Normal<-re$q
Var_Normal<-re$Var
}
} else {
q_Normal<-0
Var_Normal<-0
}
if(info[2]>0)
{
GCmat<-matrix(GCmat,length(pvalue.GC),2)
re<-Hybrid_GC(pvalue=pvalue.GC,GCmat=GCmat,CCsize=CCsize.GC,Cutoff.GC=Cutoff.GC)
q_GC<-re$qjoint
mu_GC<-re$mujoint
g_GC<-re$gjoint
NAset<-re$NAset
} else {
q_GC<-0
mu_GC<-NULL
g_GC<-NULL
NAset<-NULL
}
if(!is.null(spldata))
{
if(is.vector(spldata)) spldata<-as.matrix(t(spldata))
if(is.matrix(spldata))
{
k<-ncol(spldata)
if(k%%3!=2) stop("spldata should have (3k+2) number of columns where the first two columns represent Score and Var, next 3k represent k nodes, k values of K1, and k values of K2")
if(sum(is.na(spldata))>0) stop("Cannot have NAs in spldata, if there are different number of nodes for different studies, input spldata in list format")
spldata<-split(spldata, seq(nrow(spldata)))
}
if(!is.list(spldata))
{
stop("spldata needs to be in one of the following formats: vector (only 1 study), matrix (all studies have same number of nodes), or list (different number of nodes in different studies)")
}
}
return(Hybrid_meta(q_Normal,Var_Normal,q_GC,mu_GC,g_GC,NAset,spldata,Cutoff.meta)$p.value)
}
MetaSPA_Normal <-function(pvalue,MAF=NULL,CCsize)
{
sgn<-sign(pvalue)
pvalue<-abs(pvalue)
if(is.vector(CCsize)) CCsize<-as.matrix(t(CCsize))
ncases<-CCsize[,1]
ncontrols<-CCsize[,2]
Z<-sgn*qnorm(pvalue/2,lower.tail=FALSE)
if(is.null(MAF))
{
Neff<-ncases*ncontrols/(ncases+ncontrols)
Zmeta<-sum(sqrt(Neff)*Z)/sqrt(sum(Neff))
pval<-pchisq(Zmeta^2,lower.tail=FALSE,df=1)
return(sign(Zmeta)*pval)
} else {
Neff<-2*ncases*ncontrols*MAF*(1-MAF)/(ncases+ncontrols)
Zmeta<-sum(sqrt(Neff)*Z)/sqrt(sum(Neff))
pval<-pchisq(Zmeta^2,lower.tail=FALSE,df=1)
return(sign(Zmeta)*pval)
}
}
Hybrid_Normal<-function(pvalue,MAF,CCsize)
{
sgn<-sign(pvalue)
pvalue<-abs(pvalue)
if(is.vector(CCsize)) CCsize<-as.matrix(t(CCsize))
ncases<-CCsize[,1]
ncontrols<-CCsize[,2]
if(!is.null(MAF) & (length(MAF)!=length(pvalue) | length(which(is.na(MAF)))>0)) MAF=NULL
Z<-sgn*qnorm(pvalue/2,lower.tail=FALSE)
Neff<-2*ncases*ncontrols*MAF*(1-MAF)/(ncases+ncontrols)
Zmeta<-sum(sqrt(Neff)*Z)
Var<-sum(Neff)
return(list(q=Zmeta,Var=Var))
}
Hybrid_GC<-function(pvalue,GCmat,CCsize,Cutoff.GC)
{
truep<-abs(pvalue)
sgn<-sign(pvalue)
if(is.vector(CCsize)) CCsize<-as.matrix(t(CCsize))
ncases<-CCsize[,1]
ncontrols<-CCsize[,2]
N2<-GCmat[,1]
N1<-GCmat[,2]
k<-length(ncases)
qadj<-rep(NA,k)
yjoint<-NULL
mujoint<-NULL
gjoint<-NULL
if(length(Cutoff.GC)==1) Cutoff.GC<-rep(Cutoff.GC,k)
NAset<-NULL
for(s in 1:k)
{
n<-ncases[s]+ncontrols[s]
mu<-rep(ncases[s]/n,n)
y<-c(rep(1,ncases[s]),rep(0,ncontrols[s]))
g<-c(rep(2,N2[s]),rep(1,N1[s]),rep(0,n-N2[s]-N1[s]))
g1<-g-mean(g)
m1<-mean(g)*ncases[s]
var1<-mu[1]*(1-mu[1])*sum(g1^2)
if(truep[s]>pnorm(Cutoff.GC[k],lower.tail=F)*2)
{
qadj[s]<-qnorm(truep[s]/2,lower.tail=F)*sgn[s]*sqrt(var1)
yjoint<-c(yjoint,y)
mujoint<-c(mujoint,mu)
NAset<-c(NAset,length(gjoint)+N2[s]+N1[s]+1:(n-N2[s]-N1[s]))
gjoint<-c(gjoint,g1)
} else {
gcent<-0:2-mean(g)
GenC=c(n-N2[s]-N1[s],N1[s],N2[s])
rootadj<-try(uniroot(pval_MAC1,sgn[s]*c(0,2*sum(y)),p=mu[1],m1=0,var1=var1,GenC=GenC,gcent=gcent,truep[s],Cutoff=Cutoff.GC[s])$root,silent=T)
if(class(rootadj) != "try-error")
{
qadj[s]<-rootadj
yjoint<-c(yjoint,y)
mujoint<-c(mujoint,mu)
NAset<-c(NAset,length(gjoint)+N2[s]+N1[s]+1:(n-N2[s]-N1[s]))
gjoint<-c(gjoint,g1)
}
}
}
qjoint<-sum(qadj,na.rm=T)
return(list(qjoint=qjoint,mujoint=mujoint,gjoint=gjoint,NAset=NAset))
}
Hybrid_meta<-function(q_Normal,Var_Normal,q_GC,mu_GC,g_GC,NAset,spldata,Cutoff)
{
gNA<-g_GC[NAset]
muNA<-mu_GC[NAset]
if(length(NAset)>0)
{
gNB<-g_GC[-NAset]
muNB<-mu_GC[-NAset]
} else {
gNB=NULL
muNB=NULL
}
var1<-sum(mu_GC * (1-mu_GC) * g_GC^2)
NAmu= -sum(gNB*muNB)
NAsigma=var1-sum(muNB*(1-muNB)*gNB^2)+Var_Normal
if(is.null(spldata))
{
q_spl=0
Var.spl=0
nodes=NULL
splfuny1=NULL
splfuny2=NULL
} else {
k<-length(spldata)
Score.spl<-rep(NA,k)
Var.spl<-rep(NA,k)
nodes<-list()
splfuny1<-list()
splfuny2<-list()
for(s in 1:k)
{
nk<-length(spldata[[s]])
if(nk%%3!=2) stop("rows/list elements of spldata must have (3k+2) elements, where k is the number of nodes")
nn<-(nk-2)/3
Score.spl[s]<-spldata[[s]][1]
Var.spl[s]<-spldata[[s]][2]
nodes[[s]]<-spldata[[s]][2+1:nn]
splfuny1[[s]]<-spldata[[s]][2+nn+1:nn]
splfuny2[[s]]<-spldata[[s]][2+2*nn+1:nn]
int<-findInterval(0,nodes[[s]])
nodes[[s]]<-c(nodes[[s]][1:int],0,nodes[[s]][(int+1):length(nodes[[s]])])
splfuny1[[s]]<-c(splfuny1[[s]][1:int],0,splfuny1[[s]][(int+1):length(splfuny1[[s]])])
splfuny2[[s]]<-c(splfuny2[[s]][1:int],Var.spl[s],splfuny2[[s]][(int+1):length(splfuny2[[s]])])
dup<-which(duplicated(nodes[[s]]))
nodes[[s]][dup]<-NA
splfuny1[[s]][dup]<-NA
splfuny2[[s]][dup]<-NA
}
q_spl=sum(Score.spl)
}
q<-q_Normal+q_GC+q_spl
qinv = -q
var2<-var1+sum(Var.spl)
# Noadj
pval.noadj<-pchisq(q^2/var2, lower.tail = FALSE, df=1)
Is.converge=TRUE
if(Cutoff < 10^-1){
Cutoff=10^-1
}
if(abs(q)/sqrt(var2) < Cutoff ){
pval=pval.noadj
} else {
out.uni1<-getroot_K1_meta(0, q=q,gNB=gNB,muNB=muNB,NAmu=NAmu,NAsigma=NAsigma,nodes=nodes,splfuny1=splfuny1,splfuny2=splfuny2)
out.uni2<-getroot_K1_meta(0, q=qinv,gNB=gNB,muNB=muNB,NAmu=NAmu,NAsigma=NAsigma,nodes=nodes,splfuny1=splfuny1,splfuny2=splfuny2)
if(out.uni1$Is.converge==TRUE && out.uni2$Is.converge==TRUE)
{
p1<-tryCatch(Get_Saddle_Prob_meta(out.uni1$root, q,gNB,muNB,NAmu,NAsigma,nodes,splfuny1,splfuny2),error=function(e) {return(pval.noadj/2)})
p2<-tryCatch(Get_Saddle_Prob_meta(out.uni2$root, qinv,gNB,muNB,NAmu,NAsigma,nodes,splfuny1,splfuny2),error=function(e) {return(pval.noadj/2)})
pval = abs(p1)+abs(p2)
Is.converge=TRUE
} else {
print("Error_Converge")
pval<-pval.noadj
Is.converge=FALSE
}
}
if(pval!=0 && pval.noadj/pval>10^3)
{
return(Hybrid_meta(q_Normal,Var_Normal,q_GC,mu_GC,g_GC,NAset,spldata,Cutoff=Cutoff*2))
} else {
return(list(p.value=sign(q)*pval, p.value.NA=sign(q)*pval.noadj, Is.converge=Is.converge))
}
}
Korg_meta<-function(t, gNB,muNB,NAmu,NAsigma,nodes,splfuny1,splfuny2)
{
n.t<-length(t)
out<-rep(0,n.t)
for(i in 1:n.t){
t1<-t[i]
if(is.null(gNB))
{
out[i]<-NAmu*t1+0.5*NAsigma*t1^2
} else {
temp<-log(1 - muNB + muNB * exp(gNB* t1))
out[i]<-sum(temp)+NAmu*t1+0.5*NAsigma*t1^2
}
if(!is.null(nodes))
{
out[i]<-out[i]+integrate(K1_adj_meta,0,t1,q=0,gNB=NULL,muNB=NULL,NAmu=0,NAsigma=0,nodes=nodes,splfuny1=splfuny1,splfuny2=splfuny2)$value
}
}
return(out)
}
K1_adj_meta<-function(t, q,gNB,muNB,NAmu,NAsigma,nodes,splfuny1,splfuny2)
{
n.t<-length(t)
out<-rep(0,n.t)
for(i in 1:n.t){
t1<-t[i]
temp3<-NAmu+NAsigma*t1
if(is.null(gNB))
{
out[i]<-temp3-q
} else {
temp1<-(1 - muNB)* exp(-gNB * t1) + muNB
temp2<-muNB *gNB
out[i]<-sum(temp2/temp1)+temp3-q
}
if(!is.null(nodes))
{
k<-length(nodes)
sval=0
for(s in 1:k)
{
sfun<-splinefunH(na.omit(nodes[[s]]) , na.omit(splfuny1[[s]]) ,m=na.omit(splfuny2[[s]]))
sval<-sval+sfun(t1,deriv=0)
}
out[i]<-out[i]+sval
}
}
return(out)
}
K2_meta<-function(t, gNB,muNB,NAmu,NAsigma,nodes,splfuny1,splfuny2)
{
n.t<-length(t)
out<-rep(0,n.t)
for(i in 1:n.t){
t1<-t[i]
if(is.null(gNB))
{
out[i]<-NAsigma
} else {
temp1<-((1 - muNB)* exp(-gNB * t1) + muNB)^2
temp2<-(1-muNB) * muNB * gNB^2 * exp(-gNB*t1)
out[i]<-sum(temp2/temp1,na.rm=TRUE)+NAsigma
}
if(!is.null(nodes))
{
k<-length(nodes)
var=0
for(s in 1:k)
{
sfun<-splinefunH(na.omit(nodes[[s]]) , na.omit(splfuny1[[s]]) ,m=na.omit(splfuny2[[s]]))
var<-var+sfun(t1,deriv=1)
}
out[i]<-out[i]+var
}
}
return(out)
}
getroot_K1_meta<-function(init,q,gNB,muNB,NAmu,NAsigma,nodes,splfuny1,splfuny2,tol=.Machine$double.eps^0.25,maxiter=1000)
{
t<-init
K1_eval<-K1_adj_meta(t,q,gNB,muNB,NAmu,NAsigma,nodes,splfuny1,splfuny2)
prevJump<- Inf
rep<-1
repeat
{
K2_eval<-K2_meta(t,gNB,muNB,NAmu,NAsigma,nodes,splfuny1,splfuny2)
tnew<-t-K1_eval/K2_eval
if(is.na(tnew))
{
conv=FALSE
break
}
if(abs(tnew-t)<tol)
{
conv<-TRUE
break
}
if(rep==maxiter)
{
conv<-FALSE
break
}
newK1<-K1_adj_meta(tnew,q,gNB,muNB,NAmu,NAsigma,nodes,splfuny1,splfuny2)
if(sign(K1_eval)!=sign(newK1))
{
if(abs(tnew-t)>prevJump-tol)
{
tnew<-t+sign(newK1-K1_eval)*prevJump/2
newK1<-K1_adj_meta(tnew,q,gNB,muNB,NAmu,NAsigma,nodes,splfuny1,splfuny2)
prevJump<-prevJump/2
} else {
prevJump<-abs(tnew-t)
}
}
rep<-rep+1
t<-tnew
K1_eval<-newK1
}
return(list(root=t,n.iter=rep,Is.converge=conv))
}
Get_Saddle_Prob_meta<-function(zeta, q,gNB,muNB,NAmu,NAsigma,nodes,splfuny1,splfuny2)
{
k1<-Korg_meta(zeta, gNB,muNB,NAmu,NAsigma,nodes,splfuny1,splfuny2)
k2<-K2_meta(zeta, gNB,muNB,NAmu,NAsigma,nodes,splfuny1,splfuny2)
if(is.finite(k1) && is.finite(k2))
{
temp1<-zeta * q - k1
w<-sign(zeta) * (2 *temp1)^{1/2}
v<- zeta * (k2)^{1/2}
Z.test<-w + 1/w * log(v/w)
if(Z.test > 0){
pval<-pnorm(Z.test, lower.tail = FALSE)
} else {
pval= -pnorm(Z.test, lower.tail = TRUE)
}
} else {
pval<-0
}
return(pval)
}
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.