R/ancGLOB_pv_pts.R

ancGLOB_pv_pts <-
function(x1,x2,est=tmean,fr1=1,fr2=1,nboot=500,SEED=TRUE,iter=1000,nmin=12,MC=TRUE,alpha=.05,PRM=FALSE,pts=NULL,...){
#
#  Determine critical p-value when using the function ancGLOB and pts is specified.
#  Strategy: generage data from a normal distribution, NULL true
#  compute p-value, repeat
#  iter times (iter=1000 is default)
#
#  pts is used to indicate the covariate values where comparisons are to be made.
#  Example: pts=c(1,4,6) will compare regression lines at X=1, 4 and 6
# if pts is not specified, the function terminates with an error.
#
#
# returns:
# p.crit, the critical p-value for the specified alpha value
# if PRM=T, all p-values that were computed.
# ef.iter, the actual number of interations, which might differ from iter
# due to sample sizes where it makes no sense to compute a p-value
# based on the  generated data.
#
# Like ancGLOB_pv, only pts is specified and use data in x1 and x2
#
if(is.null(pts[1]))stop('pts is null, use ancGLOB_pv')
x1=elimna(x1)
x2=elimna(x2)
n1=length(x1)
n2=length(x2)

if(SEED)set.seed(45)
bvec=list()
np1=min(c(n1,n2))+1
nmax=max(c(n1,n2))
for(i in 1:iter){
bvec[[i]]=rmul(nmax,p=4)
if(n1!=n2)bvec[[i]][np1:nmax,1:2]=NA
bvec[[i]][1:n1,1]=x1
bvec[[i]][1:n2,3]=x2
}
prm=NA
if(MC){
library(parallel)
prm=mclapply(bvec,ancGLOB_sub4,fr1=fr1,fr2=fr2,est=est,SEED=SEED,nboot=nboot,pts=pts,...)
}
#if(!MC)prm=lapply(bvec,ancGLOB_sub4,fr1=fr1,fr2=fr2,est=est,SEED=SEED,nboot=nboot,pts=pts,...)
if(!MC){
for(ij in 1:length(bvec)){
bv=as.matrix(bvec[[ij]])
prm[ij]=ancGLOB_sub4(bv,fr1=fr1,fr2=fr2,est=est,SEED=SEED,nboot=nboot,pts=pts,nmin=nmin,...)
}}
prm=elimna(as.vector(matl(prm)))
ef.iter=length(prm)
p.crit=hd(prm,alpha)
prm=sort(elimna(prm))
if(!PRM)prm=NULL
list(p.crit=p.crit,prm=prm,ef.iter=ef.iter)
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.