Nothing
"emaxsimB" <-
function(nsim, genObj, prior, modType=4,binary=FALSE,seed=12357,
check=FALSE,nproc=parallel::detectCores(),
negEmax=FALSE,ed50contr=NULL, lambdacontr=NULL,testMods=NULL,
idmax=length(doselev),mcmc=mcmc.control(),customCode=NULL,
customParms=NULL,
description="")
{
#### assumes dose levels are sorted on input
#### placebo (comparator) is in the first position
#### primary test dose is in position idmax
if(! modType %in% c(3,4))stop("modType must be 3 or 4")
if(length(ed50contr)!=length(lambdacontr))stop('The number of ED50 and Lambda defining contrasts must be equal')
if(isTRUE( binary!=prior$binary ))stop('Binary specification in prior and model do not match')
if(inherits(prior,'emaxPrior'))localParm<-TRUE else localParm<-FALSE
if(exists('.Random.seed'))save.seed<-.Random.seed
save.rng<-RNGkind()[1]
on.exit( RNGkind(save.rng))
on.exit(if(exists('save.seed')).Random.seed<<-save.seed,add=TRUE)
### extract design parameters from genObj
n<-genObj$genP$n
doselev<-genObj$genP$doselev ## doselev should be sorted/unique
Ndose<-length(doselev)
if(!binary)ddf<-sum(n)-Ndose else ddf<-Inf
nfrac<-n/sum(n) ##allocate fractional obs for stability
nfrac2<-0.5*nfrac
dose<-genObj$genP$dose
if(localParm)estan<-selEstan('mrmodel') else estan<-selEstan('basemodel')
### set up emax-model contrasts for null hypothesis test
contMat<-NULL
if(missing(testMods)){
if(is.null(ed50contr)){
if(Ndose<=4){
ed50contr<-c((doselev[1]+doselev[2])/2,
(doselev[Ndose-1]+doselev[Ndose])/2)
lambdacontr<-rep(1,2)
}else{
ed50contr<-c((doselev[1]+doselev[2])/2,median(doselev),
(doselev[Ndose-1]+doselev[Ndose])/2)
lambdacontr<-rep(1,3)
}
}
parmscontr<-cbind(ed50contr,lambdacontr)
testMods<-Mods(sigEmax=parmscontr,doses=doselev, placEff = 0,
maxEff=1-2*(negEmax))
}
if(!binary) contMat <-optContr(testMods,w=n)
### simulation result holders
if(modType==3){nparm<-3
}else {nparm<-4}
### posterior intervals to compute and store
llev<-c(0.025,0.05,0.1)
ulev<-c(0.975,0.95,0.9)
### set up independent stream of random numbers for
### each simulation iteration.
RNGkind("L'Ecuyer-CMRG")
set.seed(seed)
rseed<-matrix(integer(nsim*7),ncol=7)
rseed[1,]<-as.integer(.Random.seed)
if(nsim>1){
for(i in 2:nsim){
rseed[i,]<-nextRNGStream(rseed[i-1,])
}
}
if(isTRUE(.Platform$OS.type=='unix') && missing(nproc))
stop('nproc must be specified on multi-user unix machines')
if(check)nproc<-1 else{
if(nproc>detectCores())stop("The number of processes requested exceeds the number of processors available.")
if(nproc>nsim){
warning(paste('The number of processors requested exceeded the number of simulations\n',
'This iikely a mistake. nproc set to 1'))
nproc<-1
}
}
if(mcmc$chains>1 && !check)stop('The number of chains should be 1 except when testing.')
### set up indices for computing consecutive blocks
### of simulations
nblock<-as.integer(trunc(nsim/nproc))
nleft<-nsim%%nproc
indmat<-matrix(integer(nproc*2),ncol=2)
indmat[1,1]<-1
indmat[1,2]<-nblock+1*(1<=nleft)
if(nproc>=2){
for (i in 2:nproc){
indmat[i,1]<-indmat[i-1,2]+1
indmat[i,2]<-indmat[i,1]+nblock-1 + 1*(i<=nleft)
}
}
inlist<-list(indmat=indmat,rseed=rseed,Ndose=Ndose,dose=dose,
ddf=ddf,doselev=doselev,nparm=nparm,modType=modType,
binary=binary,genObj=genObj,testMods=testMods,
contMat=contMat,negEmax=negEmax,
check=check,estan=estan,prior=prior,mcmc=mcmc,
customCode=customCode,customParms=customParms,
nfrac=nfrac,nfrac2=nfrac2,n=n,
ulev=ulev,llev=llev,localParm=localParm)
if(nproc==1){
simout<-simrepB(1,inlist)
if(check)return(simout)
simout<-list(simout)
}else{
cl<-makeCluster(nproc)
registerDoParallel(cl)
simout<-foreach(j=1:nproc, .packages='clinDR') %dopar%{
simrepB(j,inlist)
}
stopCluster(cl)
}
####################################
### assign output to matrices/vectors
predpop <- matrix(rep(NA, nsim * Ndose), ncol = Ndose)
fitpredv <- matrix(rep(NA, nsim * Ndose), ncol = Ndose)
fitdifv <- matrix(rep(NA, nsim * Ndose), ncol = Ndose)
sepredv <- matrix(rep(NA, nsim * Ndose), ncol = Ndose)
sedifv <- matrix(rep(NA, nsim * Ndose), ncol = Ndose)
mv <- matrix(rep(NA, nsim * Ndose), ncol = Ndose)
msSat<-rep(NA,nsim)
pVal<-rep(NA,nsim)
divergence<-rep(NA,nsim)
gofP<-rep(NA,nsim)
selContrast<-rep(NA,nsim)
colnames(predpop)<-doselev
colnames(fitpredv)<-doselev
colnames(fitdifv)<-doselev
colnames(sepredv)<-doselev
colnames(sedifv)<-doselev
colnames(mv)<-doselev
residSD<-numeric(nsim)
est<-matrix( rep(NA,nsim*nparm),ncol=nparm )
if(modType==3){
colnames(est)<-c("led50","emax","e0")
}else{
colnames(est)<-c("led50","lambda","emax","e0")
}
if(localParm) est<-cbind(est,difTarget=rep(NA,nsim))
np1<-ncol(est)
nlev<-length(llev)
estlb<-array(numeric(np1*nsim*nlev),dim=c(np1,nsim,nlev))
estub<-array(numeric(np1*nsim*nlev),dim=c(np1,nsim,nlev))
dimnames(estlb)<-list(colnames(est),1:nsim,llev)
dimnames(estub)<-list(colnames(est),1:nsim,ulev)
if(!binary){
sdv <- matrix(rep(NA, nsim * Ndose), ncol = Ndose)
}else sdv<-NULL
nd1<-Ndose-1
lb<-array(numeric(nd1*nsim*nlev),dim=c(nd1,nsim,nlev))
ub<-array(numeric(nd1*nsim*nlev),dim=c(nd1,nsim,nlev))
dimnames(lb)<-list(doselev[2:Ndose],1:nsim,llev)
dimnames(ub)<-list(doselev[2:Ndose],1:nsim,ulev)
pop<-NULL
popSD<-NULL
if(is.null(customCode))customOut<-NULL
else customOut<-vector("list",nsim)
for(j in 1:nproc){
predpop[indmat[j,1]:indmat[j,2],]<-simout[[j]]$predpop
fitpredv[indmat[j,1]:indmat[j,2],]<-simout[[j]]$fitpredv
fitdifv[indmat[j,1]:indmat[j,2],]<-simout[[j]]$fitdifv
sepredv[indmat[j,1]:indmat[j,2],]<-simout[[j]]$sepredv
sedifv[indmat[j,1]:indmat[j,2],]<-simout[[j]]$sedifv
mv[indmat[j,1]:indmat[j,2],]<-simout[[j]]$mv
sdv[indmat[j,1]:indmat[j,2],]<-simout[[j]]$sdv
msSat[indmat[j,1]:indmat[j,2]]<-simout[[j]]$msSat
pVal[indmat[j,1]:indmat[j,2]]<-simout[[j]]$pVal
divergence[indmat[j,1]:indmat[j,2]]<-simout[[j]]$divergence
gofP[indmat[j,1]:indmat[j,2]]<-simout[[j]]$gofP
selContrast[indmat[j,1]:indmat[j,2]]<-simout[[j]]$selContrast
residSD[indmat[j,1]:indmat[j,2]]<-simout[[j]]$residSD
est[indmat[j,1]:indmat[j,2],]<-simout[[j]]$est
estlb[,indmat[j,1]:indmat[j,2],]<-simout[[j]]$estlb
estub[,indmat[j,1]:indmat[j,2],]<-simout[[j]]$estub
lb[,indmat[j,1]:indmat[j,2],]<-simout[[j]]$lb
ub[,indmat[j,1]:indmat[j,2],]<-simout[[j]]$ub
customOut[indmat[j,1]:indmat[j,2]]<-simout[[j]]$customOut
pop<-rbind(pop,simout[[j]]$pop)
popSD<-c(popSD,simout[[j]]$popSD)
}
return(structure(list(description=description,localParm=localParm,
binary=binary,modType=modType,genObj=genObj,
pop=pop,popSD=popSD,mcmc=mcmc,prior=prior,
est=est,estlb=estlb,estub=estub,residSD=residSD,
pVal=pVal,selContrast=selContrast,
testMods=testMods,gofP=gofP,
negEmax=negEmax,
predpop=predpop,
mv = mv, sdv = sdv, msSat=msSat, fitpredv = fitpredv,
sepredv = sepredv, fitdifv = fitdifv, sedifv = sedifv,
lb=lb,ub=ub,divergence=divergence,
rseed=rseed, idmax=idmax,customOut=customOut
),class="emaxsimB") )
}
simrepB<-function(j,inlist)
{
indmat<-inlist$indmat
rseed<-inlist$rseed
Ndose<-inlist$Ndose
dose<-inlist$dose
ddf<-inlist$ddf
doselev<-inlist$doselev
nparm<-inlist$nparm
modType<-inlist$modType
binary<-inlist$binary
genObj<-inlist$genObj
testMods<-inlist$testMods
contMat<-inlist$contMat
negEmax<-inlist$negEmax
check<-inlist$check
estan<-inlist$estan
prior<-inlist$prior
mcmc<-inlist$mcmc
customCode<-inlist$customCode
customParms<-inlist$customParms
nfrac<-inlist$nfrac
nfrac2<-inlist$nfrac2
n<-inlist$n
llev<-inlist$llev
ulev<-inlist$ulev
localParm<-inlist$localParm
nrep<-indmat[j,2]-indmat[j,1]+1
### set up input dose variable
if(binary)din<-c(doselev,doselev) else din<-doselev
if(binary)sigsim<-1 ### default placeholder
predpop <- matrix(rep(NA, nrep * Ndose), ncol = Ndose)
fitpredv <- matrix(rep(NA, nrep * Ndose), ncol = Ndose)
fitdifv <- matrix(rep(NA, nrep * Ndose), ncol = Ndose)
sepredv <- matrix(rep(NA, nrep * Ndose), ncol = Ndose)
sedifv <- matrix(rep(NA, nrep * Ndose), ncol = Ndose)
mv <- matrix(rep(NA, nrep * Ndose), ncol = Ndose)
msSat<-rep(NA,nrep)
pVal<-rep(NA,nrep)
divergence<-rep(NA,nrep)
selContrast<-rep(NA,nrep)
colnames(predpop)<-doselev
colnames(fitpredv)<-doselev
colnames(fitdifv)<-doselev
colnames(sepredv)<-doselev
colnames(sedifv)<-doselev
colnames(mv)<-doselev
residSD<-numeric(nrep)
gofP<-numeric(nrep)
est<-matrix( rep(NA,nrep*nparm),ncol=nparm )
if(modType==3){
colnames(est)<-c("led50","emax","e0")
}else{
colnames(est)<-c("led50","lambda","emax","e0")
}
if(localParm) est<-cbind(est,difTarget=rep(NA,nrep))
np1<-ncol(est)
nlev<-length(llev)
estlb<-array(numeric(np1*nrep*nlev),dim=c(np1,nrep,nlev))
estub<-array(numeric(np1*nrep*nlev),dim=c(np1,nrep,nlev))
dimnames(estlb)<-list(colnames(est),1:nrep,llev)
dimnames(estub)<-list(colnames(est),1:nrep,ulev)
if(!binary){
sdv <- matrix(rep(NA, nrep * Ndose), ncol = Ndose)
}else sdv<-NULL
nd1<-Ndose-1
lb<-array(numeric(nd1*nrep*nlev),dim=c(nd1,nrep,nlev))
ub<-array(numeric(nd1*nrep*nlev),dim=c(nd1,nrep,nlev))
dimnames(lb)<-list(doselev[2:Ndose],1:nrep,llev)
dimnames(ub)<-list(doselev[2:Ndose],1:nrep,ulev)
pop<-NULL
popSD<-NULL
if(is.null(customCode))customOut<-NULL
else customOut<-vector("list",nrep)
if(!negEmax)trend<-'positive' else trend<-'negative'
for(i in 1:nrep) {
ioff<-i+indmat[j,1]-1
### generate simulated data set by calling
### the function in genObj
.Random.seed<<-rseed[ioff,]
### simulate design and data
gendat <- genObj$genFun(genObj$genP)
predpop[i,]<-gendat[['meanlev']]
pop<-rbind(pop,gendat[['parm']])
popSD<-c(popSD,gendat[['resSD']])
y<-gendat[['y']]
#### set up normal approximation to logit(phat) if
#### binary data
####
if(binary){
phat<-as.numeric(tapply(y,dose,sum))
phat<-(phat+nfrac2)/(n+nfrac) ###shrink to .5 for stability
yhat<-qlogis(phat)
V<-diag(1/(phat*(1-phat)*n))
contMat <-optContr(testMods,S=V)
}else{
yhat<-as.numeric(tapply(y,dose,mean))
ms<-summary(lm(y~factor(dose)))$sigma
msSat[i]<-ms^2
V<-diag(msSat[i]/n)
}
### compute p-value for preliminary test
holdC<-MCTtest(doselev, yhat, contMat = contMat,
S=V,df=ddf,type='general')$tStat
holdP<-attr(holdC,'pVal')
orderh<-order(holdC)
ncontr<-length(orderh)
pVal[i]<-holdP[orderh[ncontr]]
selContrast[i]<-orderh[ncontr]
### format data as counts if binary
if(binary){
cin<-tapply(y,dose,sum)
cin<- c(cin,n-cin)
yin<-c(rep(1,Ndose),rep(0,Ndose))
}else{
cin<-n
yin<-yhat
}
bfit<-fitEmaxB(yin,din,prior,modType,count=cin,binary=binary,
msSat=msSat[i],mcmc=mcmc,estan=estan,
diagnostics=check,nproc=1)
### return mcmc for preliminary checking
if(check){
if(binary){
return(list(estanfit=bfit$estanfit,parms=coef(bfit),pVal=pVal[i],dose=dose,y=y))
}else{
return(list(estanfit=bfit$estanfit,parms=coef(bfit),residSD=sigma(bfit),pVal=pVal[i],dose=dose,y=y))
}
}
### assign simulation output
sampler <- get_sampler_params(bfit$estanfit, inc_warmup = FALSE)
sampler<-sampler[[1]]
divergence[i] <- mean(sampler[,'divergent__'])
mv[i, ] <- tapply(y,dose,mean)
if(!binary){
sdv[i, ] <- sqrt(tapply(y,dose,var))
}
simout<-predict(bfit,doselev)
fitpredv[i,]<-simout$predMed
fitdifv[i,]<-simout$fitdifMed
sepredv[i, ] <-simout$se
sedifv[i, ] <-simout$sedif
### extract generated parameters
parms<-coef(bfit)
if(!binary){
sigsim<- sigma(bfit)
residSD[i]<-median(sigsim)
}
if(localParm){
parms1<- cbind(parms,coef(bfit,local=TRUE))
}else parms1<-parms
est[i,]<- apply(parms1,2,median)
simdif<-simout$simResp[,2:Ndose]-simout$simResp[,1]
for(k in 1:nlev){
qtile<-apply(simdif,2,quantile,c(llev[k],ulev[k]))
lb[,i,k]<-qtile[1,]
ub[,i,k]<-qtile[2,]
qtile<-apply(parms1,2,quantile,c(llev[k],ulev[k]))
estlb[,i,k]<-qtile[1,]
estub[,i,k]<-qtile[2,]
}
### compute gof test
gofP[i]<-bpchkMonoEmax(bfit,trend=trend)
### execute customized code
if(!is.null(customCode)){
if(!is.null(customParms)){
if(binary)customOut[[i]]<-customCode(parms,pVal[i],dose,y,
customParms=customParms )
else customOut[[i]]<-customCode(parms,sigsim,pVal[i],dose,y,
customParms=customParms )
}else{
if(binary)customOut[[i]]<-customCode(parms,pVal[i],dose,y)
else customOut[[i]]<-customCode(parms,sigsim,pVal[i],dose,y)
}
}
}
return(list(pop=pop,popSD=popSD,
est=est,estlb=estlb,estub=estub,residSD=residSD,
pVal=pVal,selContrast=selContrast,
gofP=gofP,predpop=predpop,
mv = mv, sdv = sdv, msSat=msSat, fitpredv = fitpredv,
fitdifv = fitdifv,sepredv = sepredv, sedifv = sedifv, lb=lb,ub=ub,
divergence=divergence,
customOut=customOut
))
}
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.