Nothing
"summary.emaxsim" <-
function(object,testalpha=0.05,clev=0.9,seSim= FALSE,...)
{
binary<-object$binary
nsim<-length(object$fitType)
modType<-object$modType
doselev<-object$genObj$genP$doselev
Ndose<-length(doselev)
n<-object$genObj$genP$n
negEmax<-object$negEmax
dirEff<- 1-2*(negEmax)
pVal<-object$pVal
idmax<-object$idmax
fitdifv<-object$fitpredv-object$fitpredv[,1]
fitdifP<-object$predpop-object$predpop[,1]
noFit<- (modType==4 & apply(is.na(object$est4),1,any)) | (modType==3 & apply(is.na(object$est3),1,any))
propFit<-table(object$fitType)/nsim
## select best pairwise function (return index to dif)
if(!negEmax){f<-function(x){xmax<-max(x); which(x==xmax)[1]}
}else f<-function(x){xmin<-min(x); which(x==xmin)[1]}
cat(paste("\n",object$description,"\n\n"))
cat(paste("Number of simulations: ",nsim,'\n',sep=''))
if(seSim== TRUE){
cat(paste("Proportion failing to converge: ",round(mean(noFit),3),
"(",round(sqrt(mean(noFit)*(1-mean(noFit))/nsim),4),")",
"\n",sep=""))
if(object$switchMod){
cat(paste("Proportion converged but ED50>upper limit: ",round(mean(object$bigC),3),
"(",round(sqrt(mean(object$bigC)*(1-mean(object$bigC))/nsim),4),")",
"\n",sep=""))
cat(paste("Proportion converged but ED50<lower limit: ",round(mean(object$negC),3),
"(",round(sqrt(mean(object$negC)*(1-mean(object$negC))/nsim),4),")",
"\n",sep=""))
}
cat(paste("Proportions with modType:"))
print(round(propFit,2))
cat(paste("Simulation standard errors:"))
print(round(sqrt(propFit*(1-propFit)/nsim),2))
}else{
cat(paste("Proportion failing to converge: ",round(mean(noFit),3),
"\n",sep=""))
if(object$switchMod){
cat(paste("Proportion converged but ED50>upper limit: ",round(mean(object$bigC),3),
"\n",sep=""))
cat(paste("Proportion converged but ED50<lower limit: ",round(mean(object$negC),3),
"\n",sep=""))
}
cat(paste("Proportions with modType:"))
print(round(propFit,2))
}
### pairwise mean+sv ndose<-length(doselev)
mdifv <- object$mv - object$mv[, 1]
if(!binary){
vp<-scale((n[1] - 1) * object$sdv[, 1]^2 +
scale(object$sdv^2,center=FALSE,scale=1/(n-1)),center=FALSE,
scale=(n[1]+n-2))
semdifv<-sqrt(scale(vp,center=FALSE,scale=1/(1/n[1] + 1/n)))
}else{
### add 1/2 y/n if 0/1 rate for se
tol<-sqrt(.Machine$double.eps)
tmpm<-object$mv
for(i in 1:Ndose){
tmpm[object$mv[,i]<tol,i]<- 0.5/(n[i]+1)
tmpm[object$mv[,i]>1-tol,i]<-(n[i]+0.5)/(n[i]+1)
}
semdifv<-sqrt(tmpm[,1]*(1-tmpm[,1])/n[1] +
scale(tmpm*(1-tmpm),center=FALSE,scale=n))
}
### power
cat(paste("\nPower for 1-sided tests at level ", testalpha," :",sep=""))
pow<- mean(pVal<=testalpha)
mean.pow <- mean(dirEff*mdifv[,idmax]/semdifv[,idmax] >qnorm(1-testalpha))
if(seSim== TRUE){
cat(paste("\n Global null test based on MCP-mod: ",round(pow,3),
"(",round(sqrt(mean(pow)*(1-mean(pow))/nsim),4),")",
"\n",sep=""))
cat(paste(" Pairwise comparison (simple, no model) of high dose vs pbo: ",round(mean.pow,3),
"(",round(sqrt(mean(mean.pow)*(1-mean(mean.pow))/nsim),4),")",
"\n",sep=""))
}else{
cat(paste("\n Global null test based on MCP-mod: ",round(pow,3),
"\n",sep=""))
cat(paste(" Pairwise comparison (simple, no model) of high dose vs pbo: ",round(mean.pow,3),
"\n",sep=""))
}
### confidence interval coverage and bias
actlev <- apply(abs(fitdifv - fitdifP)/object$sedifv <=qnorm(clev+(1-clev)/2),2,mean)
mean.actlev<- apply(abs(mdifv - fitdifP)/semdifv <=qnorm(clev+(1-clev)/2),2,mean)
seldif<-apply(mdifv[,2:Ndose],1,f)
seldif<-cbind(seq_along(seldif),seldif+1) #skip pbo
bestm<-mdifv[seldif]
bestpop<-fitdifP[seldif]
sel.actlev <- mean(abs(bestm - bestpop)/semdifv[seldif]<=
qnorm(clev+(1-clev)/2))
sel.up.err<- mean((bestm - bestpop)/semdifv[seldif]>=
qnorm(clev+(1-clev)/2))
sel.low.err<- mean((bestm - bestpop)/semdifv[seldif]<=
-qnorm(clev+(1-clev)/2))
se.actlev<-sqrt( actlev*(1-actlev)/nsim )
se.mean.actlev<-sqrt( mean.actlev*(1-mean.actlev)/nsim )
cat(paste("\n\nCoverage probabilities for nominal ",
clev," intervals [Dose-PBO]:",sep=""))
names(actlev)<-doselev
names(mean.actlev)<-doselev
names(se.actlev)<-doselev
names(se.mean.actlev)<-doselev
cat(paste("\nDose response modeling:\n"))
print(round(actlev[-1],3))
if(seSim== TRUE){
cat(paste("Simulation standard errors:\n"))
print(round(se.actlev[-1],4))
}
cat(paste("Pairwise comparisons:\n"))
print(round(mean.actlev[-1],3))
if(seSim== TRUE){
cat(paste("Simulation standard errors:\n"))
print(round(se.mean.actlev[-1],4))
}
cat(paste("Most favorable pairwise comparison:\n"))
cat(paste(round(sel.actlev,3),' Intervals too low(',
round(sel.low.err,3),
') Intervals too low(',round(sel.up.err,3)),')\n',sep='')
bias<-apply(fitdifv-fitdifP,2,mean)
se.bias<-sqrt( apply(fitdifv,2,var)/nsim )
names(bias)<-doselev
names(se.bias)<-doselev
cat(paste("\n\nBias from dose response modeling [DOSE-PBO, EST-POP]:\n"))
print(round(bias[-1],2))
if(seSim== TRUE){
cat(paste("Simulation standard errors:\n"))
print(round(se.bias[-1],3))
}
sel.mbias<-mean(bestm-bestpop)
cat(paste("Bias in the most favorable pairwise comparison:\n"))
cat(paste(round(sel.mbias,2),'\n'))
### mean squared errors
sd.sedifv <- apply(object$sedifv, 2, mean)
names(sd.sedifv)<-doselev
sd.semdifv <- apply(semdifv, 2, mean)
names(sd.semdifv )<-doselev
mse.sedifv <- sqrt( apply((fitdifv-fitdifP)^2,2,mean) )
names(mse.sedifv)<-doselev
mse.pair<-sqrt( apply((mdifv-fitdifP)^2,2,mean) )
names(mse.pair)<-doselev
sel.mse.pair<-sqrt( mean((bestm-bestpop)^2) )
cat(paste("\nSquare Root Mean Squared Error:\n"))
cat(paste("Dose response modeling:","\n"))
print(round(mse.sedifv[-1],3))
cat(paste("\n","Pairwise comparisons:","\n"))
print(round(mse.pair[-1],3))
cat(paste("\nMost favorable pairwise comparison:","\n"))
cat(paste(round(sel.mse.pair,3)))
if(seSim== TRUE){
cat(paste("\nNote: (Standard errors) are simulation errors based ",
"on ",nsim," simulations","\n",sep=""))
cat(paste("Note: Distinguish (Standard errors) from ",
"reported non-linear estmation SE" ,"\n",sep=""))
}
return(invisible(list(propNoFit=mean(noFit),
propSwitchMod=mean(object$bigC),
propHighED50=mean(object$bigC),
propLowED50=mean(object$negC),
propFit=propFit,powMCPMOD=pow,
PowMean=mean.pow,covMod=actlev[-1],
covMean=mean.actlev[-1],covSelMean=sel.actlev,
errLowSelMean=sel.low.err,errUpSelMean=sel.up.err,
biasMod=bias[-1],biasSelMean=sel.mbias,
summarySEmod=sd.sedifv,summarySEmean=sd.semdifv,
mseMod=mse.sedifv[-1],mseMean=mse.pair[-1],
mseSelMean=sel.mse.pair)))
}
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.