#! /bin/sh
# This is a shell archive. Remove anything before this line, then unpack
# it by saving it into a file and typing "sh file". To overwrite existing
# files, type "sh file -c". You can also feed this as standard input via
# unshar, or by typing "sh 'README' <<'END_OF_FILE'
X THIS IS THE READ-ME FILE FOR THE
X NONLINEAR CALIBRATION DESIGN
X SHAREWARE
X
X version 5/23/95
X
X The purpose of the shareware is to allow testing through
Xsimulation of various possible designs for nonlinear calibration
Xproblems. A brief description of the nonlinear calibration
Xdesign problem is described below. When unpacked in the
XUNIX directory ~user/current, this share archive should
Xproduce the following directories and files:
X
X~user/current
X 3 files: 1 orientation [README]
X 1 source code [DDOCODE]
X 1 example code file [DDOEXAMP]
X~user/current/.Data
X 0 files.
X~user/current/.Data/.Help
X 21 files: 21 help files
X
X
X HOW TO BEGIN
X
XTo begin using these files, type
X Splus
XOnce in Splus, type
X source("DDOCODE")
XThis command should create 65 shareware-related files in
Xthe .Data directory: 64 functions, and 1 object [phrab.bestfit1].
XIt will also give a basic orientation to the system.
XTo repeat the basic orientation to the system type
X ddo
XA great deal of other help is contained in the help files. An index
Xto these help files can be obtained, once in Splus, by
Xtyping
X help(index.ddo)
X
X
X CONTINUING
X
X If you want to experiment with an example by working
Xthrough an analysis similar to the pH calibration design search
Xdescribed in section VII of the M.A. thesis referenced below,
Xyou can try copying code from the file DDOEXAMP into Splus.
XYou probably will not want to execute the entire DDOEXAMP file
Xat once because the large number of simulations done on
Xeach design would take a long time. If you are not interested
Xin working through this example, the DDOEXAMP file may be deleted.
X
X If you think you may want to implement your own
Xapplication you should examine the "application" help file
Xvia help(application)). It provides the primary instructions
Xfor creating your own application (assuming, of course,
Xyou are familiar with the nonlinear calibration design
Xproblem as described below). You may also wish to examine
Xthe comments contained at the top of most of the Splus
Xfunctions, and in the source file DDOCODE;
Xfinally, you may wish to examine the top of the
Xfile DDOCODE which contains various overall comments.
XThe file DDOCODE may be removed from ~user/current after
Xthe source("DDOCODE") command has been given, since all
Xroutines and objects which it supplies will now be contained
Xin the ~user/current/.Data directory.
X
X
X PORTING TO SPLUS FOR WINDOWS
X
X The shareware system is compatible with Splus for windows.
XThe porting of the help-files is left to the devices of the user;
XThe DDOCODE file may be directly transferred, and sourced into
Xthe windows _data directory just as it is sourced into the UNIX
X.Data directory.
X
X
X THE NONLINEAR CALIBRATION DESIGN PROBLEM
X
X Suppose we wish to measure aspect X of a
Xcollection of objects i={1,...n}. That is, we wish to measure
XX_i for i={1,...n}. But suppose our measurement machine
Xis such that we cannot measure X directly
Xand must instead measure Y_i for i={1,...n} and backcalculate
XX_i via a relationship Y_i = F_beta(X_i) where F is a nonlinear
Xfunction of a known family and beta is a vector of parameters.
XDetermining beta is the process of calibrating the machine,
Xand is done by measuring a collection of objects of known
XX values, determining the corresponding Y values, and then
Xchoosing beta such that F_beta best fits the known (X,Y) pairs.
XThe DESIGN problem which is addressed by the shareware is the
Xproblem of choosing the X values efficiently so that the
Xresulting fitted beta will provide good backcalculated X
Xvalues.
X
X REFERENCES/QUESTIONS
X
X For a complete discussion of the computer simulation
Xapproach to the nonlinear calibration problem which this
Xshareware implements, see the (copyrighted) Master's Thesis,
X"Nonlinear Calibration Design" by Doug Oman, University of
XCalifornia at Berkeley, Department of Biostatistics, 1995.
X
X This shareware was created by Doug Oman, who hereby gives
Xpermission to use the software for noncommercial purposes
Xand to freely distribute it. At some point I hope to publish
Xthe essentials of the software simulation approach in a
Xrefereed statistics journal; until that time, users with
Xquestions about the approach behind the software are encouraged
Xto obtain a copy of the Master's thesis referenced above.
XUsers are welcome to send comments or questions regarding
Xthe software to Doug Oman,
X
X oman@stat.berkeley.edu.
X
X -Doug Oman, Berkeley, California, May, 1995
END_OF_FILE
if test 4595 -ne `wc -c <'README'`; then
echo shar: \"'README'\" unpacked with wrong size!
fi
# end of 'README'
fi
if test -f 'DDOCODE' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'DDOCODE'\"
else
echo shar: Extracting \"'DDOCODE'\" \(134753 characters\)
sed "s/^X//" >'DDOCODE' <<'END_OF_FILE'
X##############################################################################
X#VERSION 5/23/95
X#This is the code for The experimental design task:
X#To initialize, type
X# > ddobegin("phcal")
X# or simply
X# > ddobegin()
X
X#Help files are archived separately. To obtain full
X#shar archive (for UNIX) send message
X# send ddo from S
X#to
X# statlib@lib.stat.cmu.edu
X#Send (software-specific) questions (although not requests for
X#general support) or comments to Doug Oman at oman@stat.berkeley.edu
X
X#This code should work for Splus for either UNIX or MS-WINDOWS
X# (has been used extensively on WINDOWS 6.0)
X
X########## Possible future changes (though not currently planned)
X# -allow ddomake(extend=T) to indicate extending
X# -create (in-code?) table of which subroutines use
X# which variables (e.g., sigma)
X# -in-code tree structures of which routines use which other routines.
X# -create seperate CAerep,CAemax, CAemin for repeat.measures
X
X
X############## LIST OF FUNCTIONS
X
X# PH-CALIBRATION SPECIFIC FUNCTIONS:
X# (REQUIRED)
X#ddo.phcal
X# function to initialize (or reset) phcal values
X#yofx.phcal
X# find fir as a function of ph, beta
X#xofy.phcal
X# find ph as a function of fir, beta
X#simulate.phcal
X# function to simulate a calibration experiment
X#
X# (OPTIONAL)
X#print.DDO.phcal
X# augmunts printDDO by printing out application-specific
X# default parameters.
X#printdetail.phcal
X# augments generic print.ddo function with phcal-specific
X# methods.
X#
X# (USED BY OTHER PHCAL-SPECIFIC; or AUTONOMOUS)
X#sigmaiidphcal
X# calculate the overall residual sigma assuming iid.
X# (used by discrepancy.phcal)
X#designphcal
X# initializes ph-calibration design tools
X#
X# DEFAULT FUNCTIONS WHICH CAN OFTEN SUBSTITUTE
X# FOR APPLICATION-SPECIFIC ROUTINES
X#discrepancyC.default
X# calculates calibration-level discrepancies
X#discrepancyD.default
X# calculates design-level discrepancies
X#gpar.showx.default
X# lays out graphing parameters for showx.default()
X#showx.default
X# plots a design and various auxilliary information
X#showxddo.default
X# a front end for showx when design to be shown is part of discrepancy
X# distribution object
X#
X#
X# GENERAL CALIBRATION FUNCTIONS:
X# HIGH LEVEL and CALLABLE BY USER
X#ddobegin
X# initialize a particular application
X#ddo
X# generic function to reset parameter values
X#intrinsic
X# calculates and plots intrinsic aberrance, asymptotic Ses,
X# and other quantities of interest which are intrinsic to
X# the entire measurement system.
X#ddomake
X# create an estimated distribution, a "discrepancy distribution object
X#roimake
X# create a range of interest as expressed using either menu name,
X# x vector, or y vector; weights optional
X#compareddo
X# compares several distribution objects [and called by plot(many args)].
X#plotddo
X# combines showddo() and plotdensity() [and called by plot(one arg)].
X#showddo
X# Plots a design and some summary statistics from a ddo.
X#focus
X# runs several analyses on one design: plotdensity, plotstrata, pairddo
X#plotdensity
X# plots discrepancy distribution from discrepancy distribution object
X#plotstrata
X# plot discrepancy curves, stratified by quantiles
X#pairddo
X# plots correlations of pairs of discrepancy distributions
X#mwtest
X# does a Mann-Whitney rank sum test comparison between
X# two discrepancy distribution objects
X#extend
X# extends an existing ddo to more trials, more discrepancy
X# functions, or a new range of interest
X#print.ddo [can be accessed by print() or by typing name of object]
X# function to print a discrepancy distribution object in
X# a readable output
X#plot.ddo [can be accessed by plot()]
X# calls either plotddo (if one ddo argument) or compareddo(two
X# or more ddo arguments)
X#make
X# takes a list of (named) designs and creates a group of ddo's.
X#
X# DESIGN MANAGEMENT (DDOFRAME) FUNCTIONS
X#
X# names.ddo
X# names.ddoframe
X# elements
X# dimnames.ddoframe
X# length.ddoframe
X# dim.ddoframe
X# names<-.ddoframe
X# c.ddo
X# c.ddoframe
X# [.ddo
X# [.ddoframe
X# [<-.ddoframe
X# print.ddoframe
X#
X# ACCESSED BY APPLICATIONS
X#print.DDO
X# prints out (universal) default parameters
X#
X# ACCESSED PRIMARILY BY OTHER GENERAL FUNCTIONS
X#xofydist
X# applies measurement system to a y vector, obtaining medians, etc.
X#lindim
X# computes the dimensions of a single character string with
X# embedded newline characters
X#chdim
X# computes the dimensions of a table of character strings in terms
X# of characters
X#tableplot
X# plots a table of character strings onto a specified area of plot
X# called only by compareddo, potentially of interest to applications
X#rankddo
X# ranks a list of ddos according to specified criteria
X#getqua
X# gets quantile estimates given discrepancy distribution (& solves a
X# windows/unix discrepancy)
X#getquaci
X# gets confidence intervals for quantiles of a discrepancy distribution
X#gpar.plotdensity
X# calculates preferred plotting limits for a discrepancy distribution
X# object
X#at
X# a function to create desired tag "(n=...)" from discrepancy dist. obj.
X#print.ddoorientation
X# a function to print out information orienting the user to the
X# calibration design set of routines
X#yofx
X#xofy
X#xofydist
X#discrepancyC
X#discrepancyD
X#gpar.showx
X#showx
X#showxddo
X# these are generic routines and call upon the appropriate
X# version of themselves (eg, yofx.phcal) as determined by
X# .DDO$application
X
X#
X#OTHER OBJECTS USED IN THIS CODE BESIDES THE CODE ITSELF:
X# ddo.phcal uses
X# phrab.bestfit1 for mu,sigma0,roi
X#
X#
X
X#OVERVIEW OF PLOTTING ROUTINES (from lowest to highest levels)
X#
X# SHOWING DESIGN:
X# gpar.showx - from several specific arguments (excluding xdesign),
X# get plot control parameters (no plot generated)
X# showx.default - from xdesign and several specific arguments with
X# defaults in .DDO, plot a design including auxillary
X# app-specific information (such as number of batches)
X# showxddo.default - from ddo, call showx, setting specific arguments
X# to values extracted from ddo
X# showddo - calls showxddo and also plots DDall table.
X#
X# SHOWING DENSITY:
X# gpar.plotdensity - from ddo, get plotting limits for density plots
X# (no plot generated)
X# plotdensity - from ddo, plot densities
X#
X# HYBRID (DESIGNS AND DENSITIES):
X# plotddo - from ddo, call showddo and then plotdensity
X# compareddo - from several ddo's, call showxddo and plotdensity
X# several times to do various comparisons
X# plot.ddo - determine whether one or several ddo's are given
X# as arguments, and call plotddo or compareddo accordingly
X#
X#
X#
X#
X##############################################################################
X#NONPORTABILITY IDIOSYNCRACY:
X#The system time functions have different names in UNIX and WINDOWS
Xoursys.time <- function(arg, opt = c("unix.time", "dos.time")){
X #Function to get time of execution from either unix or DOS
X os <- NA
X pbegin <- proc.time()
X for(nm in opt)
X if(exists(nm)) os <- nm
X if(!is.na(os)) {
X out <- list(do.call(os, list(arg)))
X names(out) <- os
X }
X else{
X do.call(os,list(arg))
X out <- list()
X }
X out$proc.time <- proc.time()-pbegin
X unlist(out)
X}
X##############################################################################
X#Initial data for experimental design.
Xphrab.bestfit1 <-list(
X coef=list(fixed=c(b1=0.9190219,b2=0.884284,b3=6.708686)),
X var.ran=0.0003479892,
X sigma=0.02919418
X )
X#####################################################################
Xddo.phcal <- function(
X mu = phrab.bestfit1$coef$fixed,
X n = 200,
X nforce = F,
X fam0 = "quasi-logistic",
X family = "quasi-logistic",
X batch = 2,
X N = 1,
X sigma0 = c(sqrt(phrab.bestfit1$sigma^2 + phrab.bestfit1$var.ran)),
X roi = "clinical",
X Xconceivable = c(5, 8.5),
X measurement.system = "truncation",
X CDname.store = c("CD1", "CD2", "CDinf","CA"),
X droperr = F,
X Xplotrange = c(5, 8), coscale = "allsame",
X Xdesignrange = c(5,8.5),
X CDname.all = discrepancyC(getCDname.all = T),
X roimenu = list(clinical = list(x = seq(7.4, 7.9, length = 5)),
X experimental = list(x = seq(6.9, 7.9, length = 9))),
X coscalemenu = list(allsame = list(c(CD1 = 1, CD2 = 1, CDinf = 1),
X c(CA = 1, CAmin = 1, CAmax = 1)),
X maxsmall = list(c(CD1 = 1, CD2 = 1, CDinf = 2/3),
X c(CA = 1, CAmin = 1, CAmax = 1))))
X{
X#function to initialize parameters for discrepancy software.
X#These parameters will be stored in object .DDO. Their meanings:
X# STATE OF NATURE:
X# sigma0: Default sigma used for design simulation
X# (simulate.phcal)
X# mu: Default mu used for design simulation (simulate.phcal)
X# SIMULATION CONTROL:
X# n: Default number of trials for random simulations
X# nforce: Default on whether number n is forced
X# fam0: Default generating family - function family for
X# generating the data.
X# family: Default analysis family - ie, function family for
X# fitting calibration curves. (xofy.phcal,..?)
X# yofx.phcal, ...
X# batch: Default number of (identical) batches of repeats of
X# design vector (ddomake)
X# N Default number of repeat measures taken of each eye.
X# EVALUATION CONTROL:
X#Xconceivable
X# Default Xconceivable range controlling aberrance cutoff points ()
X#measurement.system
X# Default "complete measurement system" to be used.
X# must be either (1) repeat.measures or (2) truncation
X#Xdesignrange
X# Default range for possible design points. Outside this,
X# software refuses to allow design points (ddomake)
X# roi Default range of interest for discrepancy function (discrepancy)
X#CDname.store
X# Default name of all discrepancy function names to store
X# automatically (ddomake)
X# droperr Default flag about whether to omit errors in calculation
X# of quantile-based Design-level Discrepancies (eg DDmedian)
X# GRAPHICS CONTROL:
X# coscale Default scaling groups (plotdensity)
X# Xplotrange Default x range for plotting (showx), designing (designphcal)
X# MENUS ETC.:
X#roimenu Named roi vectors (roimake)
X#CDname.all
X# Default name of all discrepancy function names (discrepancy)
X#
X#coscalemenu
X# Named coscale lists (gpar.plotdensity)
X list(sigma0 = sigma0,
X mu = mu,
X n = n,
X nforce = nforce,
X fam0 = fam0,
X family = family,
X batch = batch,
X N = N,
X roi = roi,
X Xconceivable = Xconceivable,
X measurement.system = measurement.system,
X CDname.store = CDname.store,
X droperr = droperr,
X Xplotrange = Xplotrange,
X Xdesignrange = Xdesignrange,
X coscale = coscale,
X CDname.all = CDname.all,
X roimenu = roimenu,
X coscalemenu = coscalemenu)
X}
X##############################################################################
Xyofx.phcal <- function(beta,ph,family){
X #program to calculate vector of fir values from vector of ph's
X #use b1 = fir(ph=b3) parameterization
X #We must be very careful about large ph values because of
X #a HORRENDOUS S FOR WINDOWS BUG
X #In the most recent version (8/94) of S FOR WINDOWS,
X #If you compute 10^100000 it does fine, returning "Inf"
X #But if you computer 10^100000.1, it interferes with the
X #OPERATING SYSTEM! All the plotting on the S window and
X #ANY OTHER OPEN WINDOW starts to turn multicolored and
X #illegible, and work is impossible, and the computer must
X #be rebooted. Do the STAT SCI company people in Washington
X #have good checking procedures for their code?
X b1 <- beta[1]
X b2 <- beta[2]
X b3 <- beta[3]
X fir <- switch(family,
X "quasi-logistic"=
X b1+b2*(1/(10^
X pmax(-100,pmin(100,b3-ph))
X +1)-.5),
X quadratic=(b1*ph+b2)*ph+b3,
X stop(cat("\nINVALID APPROXIMATING FAMILY (\"",
X family,"\")\n"))
X )
X fir
X }
X##############################################################################
Xxofy.phcal <- function(beta,fir,family){
X #function to calculate ph as a function of fir
X #inputs:
X # beta - a length three vector of parameters
X # fir - a scalar, matrix, or array of fir values
X #family - a character vector specifying function family
X #outputs:
X # ph - same dimension as fir, the ph that correspond.
X #use b1 = fir(ph=b3) parameterization
X #
X b1 <- beta[1]
X b2 <- beta[2]
X b3 <- beta[3]
X if(family=="quasi-logistic"){
X #use quasi-logistic family
X inlogpos <- (b2-2*b1+2*fir)/(b2+2*b1-2*fir)
X outdomain <- inlogpos < 1e-20
X inlogpos[outdomain] <- 1e-20
X ph <- b3 + log(inlogpos)/log(10)
X #when fir is out of domain of function,
X #send out appropriate na
X ph[outdomain] <- NA
X }
X else
X ph <- switch(family,
X quadratic={
X dis <- b2*b2-4*b1*b3+(4*b1)*fir
X dis <- array(sqrt(
X {d <- c(dis)
X d[d<0] <- NA
X d}), dim(as.array(dis)))
X #Select that branch of the parabola
X #for backcalculation which is closer
X #to the center of the range of interest.
X branch <- sign(mean(roimake(name=.DDO$roi)$x)+
X b2/(2*b1))
X -b2/(2*b1)+(branch/(2*abs(b1)))*dis
X }
X )
X ph
X }
X
X##############################################################################
Xsimulate.phcal <-
Xfunction(ph,batch=.DDO$batch,mu, family,fam0,sigma=sigmaiidphcal(),
X N=.DDO$N, sigmabatch=0,start=mu){
X #function to simulate a calibration experiment
X #inputs:
X # ph - design vector of ph values
X # batch - number of times the ph design vector is repeated
X # mu - the imagined true vector of betas relating fir to ph
X # sigma - the imagined true standard error of residual variance
X #sigmabatch - the imagined true standard error of batch
X # variation (assumed normally distributed)
X #family - the family used to analyze the data
X #fam0 - the family used to generate the data
X #outputs:
X # beta - a length three vector of fitted parameters
X # error - a flag, TRUE if there was failure of convergence
X # errinfo- a list of information relating to the error, if exists
X # arg - a list of the arguments to the function which are
X # unique to the application
X #
X #The usual model for generating simulation data is
X # fir = b1 + b2*(1/(10^(b3-ph)+1)-.5) + error
X # with error ~ iid as N(0,sigma^2)
X arg <- list(batch=batch,N=N,sigma=sigma,
X sigmabatch=sigmabatch,start=start)
X phuse <- rep(ph,batch)
X errbatch <- c(matrix(sigmabatch*rnorm(batch),nrow=length(ph),
X ncol=batch,byrow=T))
X errresid <- sigma*rnorm(length(phuse))
X #Generate the data using fam0
X fir <- yofx(mu,phuse,family=fam0) + errbatch + errresid
X #now analyze data:
X beta <- error <- F
X errinfo <- NULL
X if(family=="quasi-logistic"){
X #Fit using the logistic-like family.
X fit <- ms(~(fir-(
X beta[1] + beta[2]*(1/(10^(beta[3]-phuse)+1)-.5)
X ))^2,
X start=list(beta=start))
X #keep track of convergence:
X #See summary.ms, min.cvg.msg,
X #min.messages, and min.errno
X convflag <- fit$flags[1]
X error <- is.na(match(convflag,c(3,4,5,6,7)))
X #The above list of acceptable codes
X #was created by looking at
X #min.messages and min.errno
X if(error)
X errinfo <- list(message=min.cvg.msg(convflag))
X beta <- fit$par
X }
X else switch(family,
X quadratic={
X #Fit quadratic curves
X fit <- lm(fir~(phuse^2)+phuse+1,data=data.frame(fir,phuse))
X beta <- coef(fit)[c(2,3,1)]
X names(beta) <- c("b1","b2","b3")
X },
X stop(cat("\nINVALID APPROXIMATING FAMILY\n"))
X )
X
X return(beta=beta,error=error,errinfo=errinfo,arg=arg)
X }
X
X##############################################################################
Xprint.DDO.phcal <- function(DINIT=.DDO,menu=T){
X #[OPTIONAL FUNCTION: automatic test for existence done by print.DDO]
X # [in its absence, certain values simply
X # won't be printed out by ddo()]
X #function to print phcal-specific default values
X #These default values are:
X #STATE OF NATURE:
X # sigma0 sigma used for design simulation (simulate.phcal)
X #SIMULATION CONTROL:
X # batch number of identical times that design vector is repeated
X # N number of repeat measures taken of each object (eye)
X cat("STATE OF NATURE (phcal):\n")
X cat("sigma0 = ",DINIT$sigma0,"\n")
X cat("SIMULATION CONTROL (phcal):\n")
X cat("batch = ",DINIT$batch," (number of repeats of whole design)\n")
X cat("N = ",DINIT$N," (number of repeats taken of each eye)\n")
X if(menu)
X cat("FAMILIES (phcal):\n",
X "quasi-logistic: fir(ph) = b1 + b2*(1/(10^(b3-pH)+1)-.5)\n",
X "quadratic: fir(ph) = b1*pH^2 + b2*pH + b3\n")
X }
X##############################################################################
Xprintdetail.phcal <- function(ddo,which){
X if(which=="nature"){
X cat(" sigmaresid =",signif(ddo$evalcall$sigma,4))
X #cat(" sigmarabbit=",signif(ddo$evalcall$sigmarabbit),4)
X cat("\n")
X }
X if(which=="design"){
X cat(" ( x",ddo$evalcall$batch,"batches = ",ddo$evalcall$batch*
X length(ddo$evalcall$xdesign)," points)\n")
X }
X }
X##############################################################################
Xsigmaiidphcal <- function(sigmares=.DDO$sigma0,sigmaeye=0,sigmarabbit=0,
X N=.DDO$N){
X #This function is used for simulations where measurement
X #offsets from predicted value are assumed to be iid (that
X #is, independent of rabbit, eye, etc.)
X #This is function to calculate the operating (residual) sigma from
X #experimentally determined values and from the number
X #of (unsimulated) repeat measures of each eye
X #
X sqrt(sigmarabbit^2+sigmaeye^2+(sigmares^2)/N)
X }
X##############################################################################
Xdesignphcal <- function(roixy=roimake(name=.DDO$roi)){
X #function to initialize values used for creating of designs:
X rl <<- min(this <- roixy$x)
X rr <<- max(this)
X rm <<- mean(c(rl,rr))
X ol <<- .DDO$Xdesignrange[1]
X or <<- .DDO$Xdesignrange[2]
X om <<- .DDO$mu[3]
X ok <<- .DDO$mu[3]-log10(2+sqrt(3))
X #ok is the "knee" of the curve
X #ok is where the second derivative achieves its maximum, ie
X #the point of maximum curvature
X os <<- .DDO$mu[3]+log10(2+sqrt(3))
X #os is the "shoulder"
X names(os) <<- names(om) <<- names(ok) <<- NULL
X allp <<- sort(c(rl=rl,rr=rr,rm=rm,ol=ol,or=or,om=om,ok=ok,os=os))
X invisible()
X }
X##############################################################################
X#
X# DEFAULT APPLICATION-LEVEL ROUTINES
X# (SUFFICIENT FOR MOST OR ALL APPLICATIONS)
X#
XdiscrepancyC.default <-function(betahat, beta,
X roixy = roimake(mu = beta, fam = fam0),
X sigma = .DDO$sigma0,
X measurement.system = .DDO$measurement.system,
X Xconceivable = .DDO$Xconceivable,
X fam0 = .DDO$fam0,
X family = .DDO$family,
X getCDname.all = F)
X{
X#[DEFAULT ROUTINE, SUFFICIENT FOR MOST APPLICATIONS]
X#This function determines Calibration-Level Discrepancies
X#when given (at minimum) the true parameters and the fitted parameters.
X#This is the DEFAULT discrepancy routine which will be called
X#if particular applications do not have their own discrepancy
X#routine named discrepancy.[application]. For most applications
X#this default routine should be sufficient.
X#inputs:
X# betahat - imagined fitted curve
X# beta - imagined true curve
X# roixy - a complete range of interest object of name/w/x/y
X# for where to measure difference between true and
X# fitted curves
X# sigma - the residyal Y variance (used to calculate aberrance)
X#Xconceivable
X# - the conceivable X range (used for truncations)
X#measurement.system
X# - the "complete measurement system" to be used
X# [see paper]. Must be either
X# repeat.measures or
X# truncation
X# fam0 - the function family applied to x and beta to get true y
X# family - the function family applied to x and betahat to get y
X# getCDname.all
X# - if T, all the function returns is the names of the
X# discrepancy functions it produces.
X#outputs:
X#If getCDname.all=T, a vector of CD names; otherwise,
X#a vector of outcomes of various discrepancy functions
X#with each element named:
X# CD1 - the mean absolute deviation over roixy$y
X# CD2 - the root mean square deviation over roixy$y
X# CB3 - the maximum absolute deviation over roixy$y
X# CAmin - the probability of aberrance at min(roixy$y)
X# CAmax - the probability of aberrance at max(roiy)
X# CA - the maximum probability of aberrance over roixy$y
X#
X if(!getCDname.all) {
X #Get pointwise biases PB at each point in roixy$y
X # [roixy$y=[Ymin(betahat),Ymax(betahat) in paper]
X xm <- xofydist(y = roixy$y, betahat = betahat, family = family,
X sigma = .DDO$sigma,
X Xconceivable = Xconceivable,
X measurement.system = measurement.system)
X #Here are point-level bias discrepancies for points in roixy$y
X PB <- xm$Xmedian - roixy$x
X #Get calibration-level bias discrepancies
X CD1 <- sum(abs(PB) * roixy$w)
X CD2 <- sqrt(sum(PB^2 * roixy$w))
X CDinf <- max(abs(PB) * (roixy$w > 0)) #Get aberrance discrepancies
X Yroi <- range(roixy$y)
X pabt <- xm$PA[match(Yroi, roixy$y)]
X CAmin <- pabt[1]
X CAmax <- pabt[2]
X CA <- max(pabt)
X c(CD1 = CD1, CD2 = CD2, CDinf = CDinf,
X CA = CA, CAmin = CAmin, CAmax = CAmax)
X }
X else {
X#make CDname.all; include PDname.all for plotstrata
X out <- c("CD1", "CD2", "CDinf", "CA", "CAmin", "CAmax")
X attr(out, "PDname.all") <- c("PD", "PD", "PD", "PA", "PA", "PA")
X out
X }
X}
X##############################################################################
XdiscrepancyD.default <- function(CDall,use.default="DDmean",
X use.specified=c(CD1="DDmean",CD2="DDmean",CD3="DDmean",
X CA="DDinterq",CAmin="DDinterq",CAmax="DDinterq"),
X nerr=0){
X #function to calculate design-level discrepancies
X #from calibration-level discrepancies
X #if droperr=F then the "errors" when convergence was not
X #possible are included in the computation of quantile-based
X #discrepancies (all except mean, sd, max)
X n <- nrow(CDall)
X DDmean <- apply(CDall,2,mean)
X if(n>1){
X DDsd <- apply(CDall,2,function(a) sqrt(var(a)))
X DDmean.se <- DDsd/sqrt(n-1)
X }
X else
X DDsd <- DDmean.se <- rep(NA,ncol(CDall))
X DDmax <- apply(CDall,2,max)
X #now get quantile-based discrepancies, taking into account
X #the errors
X ntot <- n+nerr
X rsqn <- 1/sqrt(ntot-1)
X qget <- c(.25,.5,.75,.95)*ntot/n
X qgetuse <- qget<=1
X DDq <- matrix(NA,length(qget),ncol(CDall))
X if(any(qgetuse))
X DDq[qgetuse,] <- apply(CDall,2,quantile,probs=qget[qgetuse])
X DDmedian <- DDq[2,]
X DDq95 <- DDq[4,]
X DDinterq <- DDq[3,]-DDq[1,]
X DDmedian.se <- DDq95.se <- DDinterq.se <- rep(NA,ncol(DDq))
X #if n=1 we will leave these Ses as NA
X if(n>1)
X for(i in 1:ncol(DDq)){
X dro <- DDq[,i]
X if(all(is.na(dro))){
X #no quantiles at all are available. So clearly
X #se's are all unavailable
X densm <- densq <- c(NA,NA)
X }
X else{
X #need to avoid sending density NA in from= or to=.
X #at least one quantile is available. We will just
X #fill in its value for the densities which we can't
X #get anyway
X drona <- is.na(dro)
X dro[drona] <- (c(-1,1)+range(CDall[,i]))[c(1,1,2,2)][drona]
X densm <- density(CDall[,i],n=2,from=dro[2],to=dro[4])$y
X densq <- density(CDall[,i],n=2,from=dro[1],to=dro[3])$y
X #now replace NA where it should be
X densm[drona[c(2,4)]] <- NA
X densq[drona[c(1,3)]] <- NA
X DDmedian.se[i] <- rsqn*.5/densm[1]
X DDq95.se[i] <- rsqn*sqrt(.95*.05)/densm[2]
X DDinterq.se[i] <- rsqn*sqrt(3*sum(densq^2)-
X 2*prod(densq))*.25/prod(densq)
X }
X }
X
X #combine them all:
X DDall <- rbind(DDmean,DDmean.se,DDmedian,DDmedian.se,
X DDinterq,DDinterq.se,DDq95,DDq95.se,
X DDsd,DDmax)
X #Now pull out different summary measures
X #for different discrepancies:
X usen <- match(dncd <- dimnames(CDall)[[2]],names(use.specified))
X usnam <- use.specified[usen]
X usnam[is.na(usen)] <- use.default
X userow <- match(usnam,dimnames(DDall)[[1]])
X userowse <- match(paste(usnam,"se",sep="."),dimnames(DDall)[[1]])
X DDuse <- c(DDall)[(ind <- nrow(DDall)*
X seq(from=0,by=1,length=ncol(DDall)))+userow]
X DDuse.se <- c(DDall)[ind+userowse]
X DDallnew <- rbind(DDuse,DDuse.se,DDall)
X attr(DDallnew,"usename") <- usnam
X DDallnew
X }
X##############################################################################
Xgpar.showx.default <-function(mu, fam0, roixy, sigmaband, Xplotrange,
X TABdim, FAMdim, TITdim=c(1,0),TOPdim=c(0,0),SUBdim=c(0,0),
X TITcex=2.5, TOPcex=1.5, SUBcex=1.5, TITexp=TITcex/max(TOPcex,SUBcex),
X mex=1)
X{
X#function to plan overall layout of a plot of a design.
X#The inputs are various parameters; the output is a list of
X#various parameters which can serve as input to the routine
X#showx.default() which plots a design. This routine will
X#usually only be called by the showx routine.
X#NOTE: This routine assumes that the correct value of par(mfrow)
X#has been set prior to the call to this routine. It resets
X#(but ultimately leaves unchanged) other
X#parameters such as mar and mex, but it assumes that mfrow
X#is correct.
X#inputs:
X# mu - the imagined true calibration curve
X# fam0 - the generating family
X# roixy - a complete roi list
X# Xplotrange - the overall range of X to be plotted
X# sigmaband - the distances on either side of the curve
X# (vertically) to draw lines representing the residual
X# error.
X# TABdim - the dimension of the table in numbers of characters
X# FAMdim - the dimensions of the family information in numbers of
X# characters.
X#outputs are lists with various components:
X# OVA - specifies the overall plotting limits in units of X and Y
X# (all plotting limits are specified in these limits)
X# Components OVA$x, OVA$y.
X# TAB - specifying the plotting limits within which there will be
X# room to plot a (possibly large) table of values (used by
X# compareddo). Components TAB$x, TAB$y.
X# ROI - list giving:
X# X coords of ends of range of interest (ROI$X)
X# Y coords of top and bottom of each vertical line at
X# end of range of interest (ROI$Y, a 2x2 matrix
X# with columns representing lines)
X# Y coordinates of height at which to plot symbols
X# for roi points (ROI$SYM$ht) and the fractions
X# of X and Y space available for
X# plotting the symbols (ROI$SYM$frac)
X# Y coordinates of limits within which to plot the ROI
X# name (ROI$NAM)
X# FAM - specifies the plotting limits within which can be plotted
X# information regarding the generating family, the
X# analysis family, or possibly other specifications
X# of interest. Components FAM$x, FAM$y.
X#
X#The procedure used is as follows (and in the following order):
X#MAR Determine the margins needed for each figure as a function of
X# mfrow, and the desired title and top and bottom text.
X#OVA Determine overall plotting limits (in X, Y units, as are all limits).
X#TAB Determine preferred table quadrant for table based upon
X# location of middle of curve, whether increasing or decreasing.
X#ROI Determine roi trapezoid position based upon 1. Make
X# vertically opposite (below<->above) the table.
X# Add a small margin if necessary to overall.
X#FAM Determine space for family information, making it diagonally
X# opposite the table and making sure to avoid where the ROI
X# information is plotted.
X#
X#Set control constants:
X fracsym <- c(0,.05) #how far from edge are symbols
X fractext <- c(.05,.15) #how far from edge is ROI text
X cornerlookT <- seq(.2,.8,length=30) #fractions at which to try to place
X #table corner
X lengthlookF <- c(10,10,10) #number of possible X boundary points to examine
X #in each range for family information
X overlinepenaltyF <- c(1,.8,.7) #penalty to apply to cex size for lying over line
X minoverlineokF <- .5 #minimum character expansion to allow
X#
X mfrow <- par("mfrow")
X old.par <- c(par(mex=mex),list(mar=par("mar")))
X on.exit(par(old.par))
X #set up some values for later use
X xr <- range(roixy$x)
X #
X#
X# MAR
X#
X cextit <- TITcex
X cextop <- TOPcex
X cexsub <- SUBcex
X #iterate a couple of times on the process of
X #setting the margin --> determining the
X # character expansion --> resetting margin
X for(i in 1:3){
X mar <- c(SUBdim[1]*cexsub+2.1,2.1,TITdim[1]*cextit+
X TOPdim[1]*cextop,1)
X par(mar=mar)
X #now determine the number of characters' with in the
X #margin areas in "standardized" character units.
X #cinocex is these units; this variable will be used
X #also by TAB, FAM, etc.
X cinocex <- par("cin")/par("cex")
X avx <- par("fin")[1]/cinocex[1]
X cextit <- min(avx/max(.001,TITdim[2]),TITcex)
X cextop <- min(avx/max(.001,TOPdim[2]),cextit/TITexp,TOPcex)
X cexsub <- min(avx/max(.001,SUBdim[2]),cextit/TITexp,SUBcex)
X }
X MAR <- list(mar=mar,cex=c(cextit,cextop,cexsub),
X startline=c(TOPdim[1]*cextop+(TITdim[1]-1)*cextit,
X (TOPdim[1]-1)*cextop,2),mex=mex)
X
X #get final sizes of the box within which our figure is plotted
X boxpin <- par("pin")
X#
X# OVA
X#
X ycorner <- yofx(mu, Xplotrange, family = fam0)
X bandmax <- range(sigmaband, 0, - sigmaband)
X Yplotrange <- range(c(ycorner[1] + bandmax, ycorner[2] + bandmax))
X#
X# TAB
X#
X #vector elements are possible table edges
X xlook <- Xplotrange[1] + diff(Xplotrange) * sort(cornerlookT)
X nlook <- length(xlook)
X ylook <- yofx(mu, xlook, family = fam0)
X increasing <- ylook[length(ylook)] > ylook[1]
X inc <- -1 + 2 * increasing
X #now array the possiblities:
X xtryin <- rep(xlook, 2)
X above <- c(rep(F, nlook), rep(T, nlook))
X xtryout <- c(rep(Xplotrange[1 + increasing], nlook),
X rep(Xplotrange[2 - increasing], nlook))
X ytryin <- c(ylook + bandmax[1], ylook + bandmax[2])
X ytryout <- c(rep(Yplotrange[2 - increasing], nlook),
X rep(Yplotrange[1 + increasing], nlook))
X #evaluate the possiblities:
X xfracT <- abs(xtryout - xtryin)/diff(Xplotrange)
X yfracT <- abs(ytryout - ytryin)/diff(Yplotrange)
X xcexT <- boxpin[1]*xfracT/(cinocex[1]*TABdim[2])
X ycexT <- boxpin[2]*yfracT/(cinocex[2]*TABdim[1])
X critT <- xycexT <- pmin(xcexT,ycexT)
X bestiT <- match(max(critT),critT)[1]
X TAB <- list(x=sort(c(xtryin[bestiT], xtryout[bestiT])),
X y=sort(c(ytryin[bestiT], ytryout[bestiT])),
X cex=xycexT[bestiT])
X#
X# ROI
X#
X roiabv <- !above[bestiT]
X roiin <- yofx(mu, xr, family = fam0)
X roiout <- rep(Yplotrange[1 + roiabv], 2)
X roifrac <- abs(roiin + bandmax[1 + roiabv] - roiout)/diff(Yplotrange)
X #If less than max(fracmin) is available, increase
X #top or bottom margin to make it available
X Yadd <- 0
X if((h <- min(roifrac)) < max(fractext)) {
X Yadd <- (diff(Yplotrange)*(max(fractext)-h))/(1-max(fractext))
X Yadd <- Yadd * (-1 + 2 * roiabv)
X Yplotrange[1 + roiabv] <- roiout[1] <-
X roiout[2] <- Yplotrange[1 + roiabv] + Yadd
X }
X #hts is also used below in FAM
X hts <- Yplotrange[1 + roiabv] + (Yplotrange[2 - roiabv] -
X Yplotrange[1 + roiabv]) * c(fracsym, fractext
X )
X ROI <- list(x=xr,
X y=rbind(roiin,roiout),
X symbols=sort(hts[1:2]),
X name=sort(hts[3:4])
X )
X#
X# FAM
X#
X #Arrange with width fraction, height fraction,
X #and over-line indicator in three vectors which
X #are assembled into criterion.
X outhigh <- (increasing+roiabv) %% 2
X #outhigh is a flag=1 if the x-"outside" of where
X #we are working is the high value, ie Xplotrange[2]
X ohsign <- -1+2*outhigh
X xout <- Xplotrange[1+outhigh]
X xin <- Xplotrange[2-outhigh]
X xrout <- xr[1+outhigh]
X xrin <- xr[2-outhigh]
X yout <- hts[1]
X yin <- hts[4]
X #now generate the x points to examine
X xlookF <- unique(c(seq(from=xin,to=xrin,length=lengthlookF[1]),
X seq(from=xrin,to=xrout,length=lengthlookF[2]),
X seq(from=xrout,to=xout,length=lengthlookF[3])))
X ylookF <- yofx(mu,xlookF,family=fam0) + bandmax[1+roiabv]
X #now divide them up into regions relative to roi information area
X xinmost <- sign(xrin-xlookF)==ohsign
X xoutmost <- sign(xlookF-xrout)==ohsign
X xmid <- !(xinmost | xoutmost)
X sinmost <- sum(xinmost)
X soutmost <- sum(xoutmost)
X smid <- sum(xmid)
X #now for each region generate possible sets of edges
X xlookFi <- c(rep(xlookF[xinmost],2),
X rep(xlookF[xmid],2),
X xlookF[xoutmost])
X xlookFo <- c(rep(xrin,sinmost),rep(xout,sinmost),
X rep(xrout,smid),rep(xout,smid),
X rep(xout,soutmost))
X ylookFi <- c(rep(ylookF[xinmost],2),
X rep(ylookF[xmid],2),
X ylookF[xoutmost])
X ylookFo <- c(rep(yout,sinmost),rep(yin,sinmost),
X rep(yin,2*smid),
X rep(yout,soutmost))
X overline <- c(rep(0,sinmost),rep(2,sinmost),
X rep(0,smid),rep(1,smid),
X rep(0,soutmost))
X #now see the dimensions of the box that is available and
X #the character expansions (cex) that would be called for
X xfracF <- abs(xlookFo-xlookFi)/diff(Xplotrange)
X yfracF <- abs(ylookFo - ylookFi)*(sign(ylookFo-ylookFi)==
X sign(yout-yin))/diff(Yplotrange)
X xcexF <- xfracF*boxpin[1]/(FAMdim[2]*cinocex[1])
X ycexF <- yfracF*boxpin[2]/(FAMdim[1]*cinocex[2])
X xycexF <- pmin(xcexF,ycexF)
X#THE FOLLOWING WEIGHTING SYSTEM IS
X#ARBITRARY AND COULD BE CHANGED
X critF <- xycexF*overlinepenaltyF[1+overline]^(xycexF>minoverlineokF)
X bestiF <- match(max(critF),critF)[1]
X FAM <- list(x=sort(c(xlookFi[bestiF],xlookFo[bestiF])),
X y=sort(c(ylookFi[bestiF],ylookFo[bestiF])),
X cex=xycexF[bestiF])
X#
X#return results
X#
X return(MAR=MAR,OVA=list(x=Xplotrange,y=Yplotrange),
X TAB=TAB,ROI=ROI,FAM=FAM)
X}
X##############################################################################
Xshowx.default <- function(xdesign,
X batchshow = .DDO$batch,
X mu = .DDO$mu,
X roixy = roimake(name = .DDO$roi, mu = mu, fam0 = fam0),
X sigmaband = .DDO$sigma0,
X fam0 = .DDO$fam0,
X family=.DDO$family,
X measurement.system=.DDO$measurement.system,
X Xconceivable=.DDO$Xconceivable,
X main = "UNTITLED",
X Xplotrange = .DDO$Xplotrange,
X plotpoints = 100,
X TABinfo=matrix("",1,1),
X TABcol=1,
X graphsetpar = make.graphsetpar,
X blockframe = F,
X graphsetparreturn=F,
X ...)
X{
X#function to plot a picture of a univariate experimental design
X#input:
X# xdesign - vector of design points
X#batchshow - the number of times the design vector will be duplicated
X# (default =1)
X# mu - the imagined true calibration curve
X# sigma - the imagined true residual sd
X# roixy - the range of interest over which discrepancy calculated
X# fam0 - the generating family of functions
X# main - a name to give the design (Default: "DESIGN")
X#Xplotrange - the limits of ph for plotting
X#plotpoints - the number of points to plot for discretizing line
X#graphsetpar - a list giving several important calculated values
X# for where various objects should be plotted.
X#blockframe - whether or not to block this routine's tendency to
X# automatically generate a new frame on graphical device.
X# This routine will NOT reset par(mfrow) to prior value
X# because it wants to allow further information to be
X# plotted on top of the plot it generates.
X# ... - more arguments are allowed (for the sake of plotddo)
X# but will be ignored.
X#
X#output (if graphsetparreturn=T):
X# graphsetpar
X#side effect:
X#a plot is generated as a side effect.
X#
X if(!blockframe)
X par(mfrow = c(1, 1), ask = T)
X mfrow <- par("mfrow")
X parreset <- par("adj", "ask","mex","mar")
X on.exit(par(parreset))
X#
X#Now let us create a lot of the information which will later
X#be put in various places on the plot or the margins
X#
X #create the family info
X fnam <- family
X if(is.null(fnam))
X fnam <- "[default]"
X f0nam <- fam0
X if(is.null(f0nam))
X f0nam <- "[default]"
X FAMinfo <- matrix(c(
X paste("analysis family: ",fnam),
X paste("generating family: ",f0nam),
X paste("measurement system: ",measurement.system),
X paste("Xconceivable: [",Xconceivable[1],",",
X Xconceivable[2],"]",sep="")
X ),ncol=1)
X FAMdim <- chdim(FAMinfo)$chdim
X #
X #create the title info
X #
X TITinfo <- main
X TITdim <- lindim(main)
X #
X #create the top info, just below the title
X #
X TOPinfo <-
X nbatch <- 1
X if(is.null(batchshow))
X batchshowtxt <- ""
X else {
X nbatch <- batchshow
X batchshowtxt <- paste(" across", nbatch, "batches")
X }
X TOPinfo <- paste(nbatch * length(xdesign)," design points",
X batchshowtxt, sep = "")
X TOPdim <- lindim(TOPinfo)
X #
X #create subtitle info (below the plot)
X #(no info at present)
X #
X SUBinfo <- ""
X SUBdim <- lindim(SUBinfo)
X TABdim <- chdim(TABinfo)$chdim
X #get final sizes of the box within which our figure is plotted
Xif(missing(graphsetpar))
X graphsetpar <- gpar.showx(
X mu=mu,
X fam0=fam0,
X roixy=roixy,
X sigmaband=sigmaband,
X Xplotrange=range(c(Xplotrange,roixy$x,xdesign)),
X TABdim=TABdim, FAMdim=FAMdim, TITdim=TITdim,
X TOPdim=TOPdim, SUBdim=SUBdim)
X MAR <- graphsetpar$MAR
X OVA <- graphsetpar$OVA
X TAB <- graphsetpar$TAB
X ROI <- graphsetpar$ROI
X FAM <- graphsetpar$FAM
X#
X par(mex=MAR$mex, mar=MAR$mar)
X boxpin <- par("pin")
X xplot <- seq(from = OVA$x[1], to = OVA$x[2], length = plotpoints)
X yplot <- yofx(mu, xplot, family = fam0)
X offset <- matrix(c(0, sigmaband, - sigmaband),
X length(xplot), ncl <- (2 * length(sigmaband) + 1),
X byrow = T)
X allcurvelty <- c(1, rep(4, nsb <- 2 * length(sigmaband)))
X allcurvecol <- c(1, rep(3, nsb))
X allcurve <- array(yplot, dim(offset)) + offset
X matplot(xplot, allcurve, type = "l", lty = allcurvelty,
X col = allcurvecol, ylim = OVA$y)
X if(TITdim[1])
X mtext(TITinfo,side=3,line=MAR$startline[1],cex=MAR$cex[1])
X if(TOPdim[1])
X mtext(TOPinfo,side=3,line=MAR$startline[2],cex=MAR$cex[2])
X if(SUBdim[1])
X mtext(SUBinfo,side=1,line=MAR$startline[3],cex=MAR$cex[3])
X #plot subsequent text with center adjustment:
X par(adj = 0.5) #plot design points
X txnum <- table(xdesign)
X tx <- as.numeric(names(txnum))
X text(tx, yofx(mu, tx, family = fam0), labels = txnum, lty = 1, col = 2)
X #plot range of interest
X #vertical lines for boundary
X matlines(rbind(ROI$x,ROI$x),ROI$y, type = "l", lty = 2, col = 4)
X #symbols showing weights
X symsize <- sqrt(roixy$w)
X symsize <- symsize/max(symsize)
X #now get available space for symbols in inches
X avx <- (boxpin[1]*diff(ROI$x)/diff(OVA$x))/length(roixy$w)
X avy <- boxpin[2]*diff(ROI$symbols)/diff(OVA$y)
X lsi <- .3 * min(avx, avy)
X symbols(x = roixy$x, y = rep(mean(ROI$symbols), length(roixy$x)),
X circles = symsize, inches = lsi, lty = 1, col
X = 4, add = T) #now plot the text
X tableplot(matrix(paste("ROI=",roixy$name,sep=""),1,1),
X xlim=ROI$x,ylim=ROI$name,cex=2,
X adj=.5,col=4)
X #now plot the family information
X tableplot(FAMinfo,FAM$x,FAM$y,adj=.5,cex=FAM$cex,col = 2)
X #now plot the table if it is not missing
X if(!missing(TABinfo))
X tableplot(TABinfo,TAB$x,TAB$y,adj=0,cex=TAB$cex,
X col=TABcol)
X if(!graphsetparreturn)
X invisible()
X graphsetpar
X}
X##############################################################################
Xshowxddo.default <- function(dist, main = dist$evalcall$main,
X blockframe = F, graphsetparreturn = F, ...)
X{
X#function to show the design from a distribution object
X#It is a front end to showx.default but also plots n and DDall
X# which handles the universal parameters
X# mu,roi,main
X# which can also handle the following auxilliary argument:
X# sigma
X# sigmabatch
X# batch
X sigmaband <- dist$evalcall$sigma
X if((s2 <- dist$evalcall$sigmabatch) > 0)
X sigmaband <- c(sigmaband, s2)
X if(!graphsetparreturn)
X invisible()
X do.call("showx",
X c(list(sigmaband=sigmaband,main=main,
X blockframe=blockframe,
X graphsetparreturn=graphsetparreturn,...),
X dist$evalcall[c("xdesign","batch","mu","roixy",
X "measurement.system","Xconceivable","family","fam0")]))
X}
X##############################################################################
X#
X# GENERAL CALIBRATION FUNCTIONS:
X#
X
Xddobegin <- function(application="phcal"){
X #function to begin a calibration application.
X #It sets DDO equal to .DDO.[application name]
X #.DDO.[application name] should contain initial (default) values
X #and should have .DDO$application = [application name]
X #
X resetname <- paste("ddo",application,sep=".")
X if(exists(resetname)){.DDO <<- list(application=application)
X ddo(initialize=application)
X } else stop(
X cat("\nNo Such Application. Need object named",ddoname,"\n"))
X cat("\nNOW RUNNING APPLICATION ",application,";\n",
X "To get some default parameters type\n",
X " > ddo(\"parametername1\",\"parametername2\")\n",
X "To get all default parameters type\n",
X " > ddo()\n",
X "To reset default parameters type\n",
X " > ddo(parametername1=newvalue1,parametername2=newvalue2)\n\n")
X }
X##############################################################################
Xddo <- function(...,initialize=NULL,menu=F){
X #generic function to reset parameters in .DDO
X #If called with an "initialize=" argument, this will
X # initialize values for the application with the name given
X # in the argument.
X #If called with no arguments, this will print out the
X #names of all the default variables in a pleasing format.
X #
X #There must be a ddo.[application name] for each application
X #which outputs a list of the parameters needed for that
X #application.
X #
X resetname <- callargs <- NULL
X if(is.null(initialize)){
X if(!exists(".DDO"))
X stop(
X "\nMust initialize before resetting.\nType \"ddobegin()\"\n")
X argsnew <- list(...)
X #If there are no arguments, print out the default variables.
X if(length(argsnew)==0)
X print.DDO(menu=menu)
X #What follows is resetting proper.
X #Here we must fill in all missing arguments
X #with elements of .DDO
X apname <- .DDO$application
X resetname <- paste("ddo",apname,sep=".")
X argsnew <- list(...)
X nmargspossible <- names( amatch(
X get(resetname),expression(call(resetname)) ) )
X #pull out all unnamed arguments to return to user
X unnam <- if(is.null(names(argsnew)))
X rep(T,length(argsnew))
X else
X names(argsnew)==""
X retnams <- if(length(argsnew))
X nmargspossible[match(unlist(argsnew[unnam]),
X nmargspossible)]
X else
X nmargspossible
X argsnew <- argsnew[!unnam]
X #now prepare to reuse old arguments if no new ones given
X .DDOreuse <- is.na(match(nmargspossible,names(argsnew)))
X callargs <- c(.DDO[nmargspossible[.DDOreuse]],argsnew)
X }
X else {
X #Here we are initializing:
X apname <- initialize
X resetname <- paste("ddo",apname,sep=".")
X #Here, we could do checking to make sure all needed
X #objects exist
X callargs <- list(...)
X }
X #Preset a .DDO with certain minimal values
X #to make sure that reset functin works
X .DDO <<- list(classobject={clo<-apname;class(clo)<-apname;clo})
X #Call reset function
X .DDO <<- do.call(resetname,callargs)
X #Now add two more values to .DDO which are needed
X .DDO$application <<- apname
X .DDO$classobject <<- clo
X if(!is.null(initialize))
X return()
X if(!any(unnam))
X invisible()
X #now return the elements of .DDO named in retnams:
X if(length(retnams)==1)
X .DDO[[retnams]]
X else
X .DDO[retnams]
X }
X### Give this function a special class so that when its name is typed, information will be printed:
Xclass(ddo) <- "ddoorientation"
X##############################################################################
Xintrinsic <- function(
X xcon = list(absolute=c(-Inf,Inf),
X truncation=.DDO$Xconceivable),
X mu = .DDO$mu,
X roi = .DDO$roi, roixy = roimake(name = roi, mu = mu,
X fam0 = fam0),
X xplotlim = {
X r <- range(roixy$x)
X r <- r + (r - mean(r)) * 0.5
X xconr <- range(unlist(xcon))
X c(max(r[1], xconr[1]), min(r[2], xconr[2]))
X },
X xplot = sort(c(roixy$x, seq(from = xplotlim[1], to = xplotlim[2], length = 100))), epsilonpart = 0.001,
X sigma = .DDO$sigma0,
X fam0 = .DDO$fam0,
X doplot = T)
X{
X#Function to display the intrinsic aberrance of a measurement system
X#Will display:
X# -a plot over Xplotrange of the probability of aberrance
X# -the maximum over that value
X#inputs:
X# xcon - a (named) list of the ranges of conceivable X values
X# for which a plot is desired of intrinsic aberrance
X# mu - the imagined true calibration curve
X# roi - (optional) the name of the range of interested, if menu-selected
X# roixy - a complete roi object with name/w/x/y
X# xplotlim - the limits of the X variable to be plotted
X# xplot - a vector, the values of the X variable to be plotted
X# epsilonpart - used to calculate derivative of yofx
X# sigma - the error variance used for calculation of intrinsic aberrance
X# fam0 - the function family used for calculations of Y from X
X# doplot - if T (the default), plots are generated;
X# if F (used by plotstrata), vectors are returned
X#outputs (if doplot=F):
X# PA - intrinsic aberrance. Vector same length as xplot
X# Xse - asymptotic standard error of Xmeasured.
X# Vector same length as xplot.
X parmfrow <- list(mfrow = par("mfrow"))
X on.exit(par(parmfrow))
X par(mfrow = c(2, 1))
X y0 <- yofx(mu, xplot, family = fam0)
X#plot INTRINSIC ABERRANCE
X PA <- list()
X if(doplot){
X #create the subtitle, a Xconceivable list
X sbt <- lapply(xcon,function(xc)paste("= [",xc[1],",",
X xc[2],"]",sep=""))
X sbt <- paste(names(xcon),unlist(sbt))
X sbt <- paste(c("Xconceivable: ",sbt),collapse=" ")
X plot(x=range(xplot),xlim=range(xplot),ylim=c(0,1),
X xlab="X",ylab="PAintrinsic",
X main = "INTRINSIC ABERRANCE: PAintrinsic, CAintrinsic", type="n",
X sub = sbt)
X abline(v = range(roixy$x), lty = 1, col = 4)
X }
X for(xci in seq(along=xcon)){
X xc <- xcon[[xci]]
X ylim <- range(yofx(mu, xc, family = fam0))
X zmin <- (ylim[1] - y0)/sigma
X zmax <- (y0 - ylim[2])/sigma
X indroi <- (xplot >= min(roixy$x)) & (xplot <= max(roixy$x))
X pab <- pnorm(zmin) + pnorm(zmax)
X probmax <- max(pab[indroi])
X xpmax <- xplot[(probmax == pab) & indroi]
X xpnum <- length(xpmax)
X xpran <- range(xpmax)
X xpuse <- xpran[1:xpnum]
X if(doplot) {
X matlines(xplot, pab, type = "l",
X lty = c(3,1)[1+xci%%2], col = 2)
X abline(h = probmax, lty = 2, col = 2)
X text(xplotlim[1]+diff(xplotlim)*(xci-.5)/length(xcon),
X probmax + 0.04,
X paste("CA",nxc<-names(xcon)[xci]," = P[",nxc,
X " aberrance | X = ",
X signif(xpuse, 3), " ] = ",signif(
X probmax, 3), sep=""), lty = 1, col = 2)
X }
X PA <- c(PA,list(pab))
X }
X names(PA) <- names(xcon)
X#PDlinear of Xmeasured
X y0inc <- y0 + sigma * epsilonpart
X Xmeasinc <- xofy(mu, y0inc, fam0)
X dxmeasdy <- abs((Xmeasinc - xplot)/epsilonpart)
X dxdyran <- range(dxmeasdy[indroi])
X if(doplot) {
X matplot(xplot, abs(dxmeasdy), type = "l",
X main = " LINEAR APPROXIMATION OF SD(Xmeasured): PDlinear", lty =
X 1, col = 2, xlab = "X", ylab = "SD(Xmeasured)",
X ylim = range(c(0, dxdyran)), sub =
X paste(signif((dxdyran[2]/dxdyran[1] - 1) * 100, 2),
X "% variation (", signif(dxdyran[1],
X 4), " to ", signif(dxdyran[2], 4), ") over ROI = [",
X min(roixy$x),", ",max(roixy$x),"]",sep=""))
X abline(v = range(roixy$x), lty = 1, col = 4)
X }
X if(!doplot)
X return(PA = PA, Xse = abs(dxmeasdy))
X else
X invisible()
X}
X##############################################################################
Xddomake <-function(
X xdesign,
X mu = .DDO$mu,
X roixy = switch(mode(.DDO$roi),
X character=roimake(name = .DDO$roi, mu = mu, fam0 = fam0),
X list=roimake(nwxy=.DDO$roi),
X numeric=roimake(x=.DDO$roi,mu=mu,fam0=fam0),
X stop(cat("Improper specification of range of interest:\n",.DDO$roi,"\n"))),
X measurement.system = .DDO$measurement.system,
X Xconceivable = .DDO$Xconceivable,
X CDname.store = .DDO$CDname.store,
X droperr = .DDO$droperr,
X n = .DDO$n,
X nforce = .DDO$nforce,
X fam0 = .DDO$fam0,
X family = .DDO$family,
X extend.info = list(DDget=T),
X randomseed = .Random.seed,
X main = paste("UNTITLED (", date(), ")", sep = ""),
X tag = T,
X warn = T,
X ...)
X{
X#
X#Function to create an estimated distribution of discrepancy values.
X#The output will be a "discrepancy distribution object," a list with a
X#number of components described below. The object is generated
X#by simulation using routine simulation.[application name]
X#inputs:
X# DESIGN:
X# xdesign - the design used as a basis for the simulation
X#
X# EVALUATION:
X# roixy - range of interest specifications, a complete list
X# measurement.system
X# - the measurement system used for backcalculations
X# Xconceivable
X# - the conceivable X range, used for measurement.system=truncation
X# CDname.store - the names of the types of discrepancy functions to be
X# calculated and stored (default: .DDO$CDname.store)
X# droperr - if FALSE, design-level quantile-based discrepancies (such as
X# median CD) are calculated assuming that there is
X# infinite discrepancy for fittings that failed to converge
X# (ie, the errors which are stored in ddo$error).
X#
X# SIMULATION CONTROL (GENERIC):
X# mu - the imagined "true" parameter values; used for
X# evaluating discrepancies.
X# n - the number of trials to simulate
X# nforce - true forces n succesful simulations; false means
X# n trials, some possibly with nonconvergence.
X# fam0 - the family used to generate the data
X# family - the (approximating) family used to analyze the data
X# extend.info - this argument receives information from the extend()
X# function. It is a list which may contain
X# three arguments: $beta is beta values to
X# reuse from previous simulations; $nerr is the
X# number of errors from previous simulations;
X# DDget tells whether to evaluate DDall now.
X#randomseed - the random seed with which to begin the
X# simulation (default: current .Random.seed)
X#
X# SIMULATION CONTROL (APPLICATION-SPECIFIC):
X# ... - possibility of other arguments to be passed on to
X# simulation routine. This may include:
X# starting values
X# imagined "true" values of parameters
X# imagined "true" residual SD
X# imagined "true" other variance components
X# number of batches
X#
X# OUTPUT CONTROL:
X# main - a title used by later functions (eg, plotting, etc) when
X# outputting data computed from this distribution
X# tag - if TRUE, type date and "FINISHED" when finished executing
X# (default=T)
X# warn - if TRUE, will emit run-time warnings as they happen
X#
X#outputs (contained in list):
X# beta - the results of the simulation:
X# an unordered array of beta values at which the
X# discrepancy function was calculated
X# CDall - a matrix of Calibration level discrepancy values
X# (nrow = number of reps; ncol = # discrepancies)
X# DDall - a matrix of a few Design level discrepancies
X# and summary measures (ncol = # discrepancies)
X# evalcall - an evaluated list of the elements (almost
X# all of them) used in the call
X# or generated by default.
X# call - the call to the function
X# rseedend - the ending random seed
X#
X callform <- match.call(expand.dots = T)
X call0 <- list(xdesign = xdesign,
X mu = mu,
X roixy = roixy,
X measurement.system = measurement.system,
X Xconceivable = Xconceivable,
X CDname.store = CDname.store,
X droperr = droperr,
X n = n,
X nforce = nforce,
X fam0 = fam0,
X family = family,
X randomseed = eval(randomseed),
X main = main,
X ...)
X if(warn)
X options(warn = 1)
X#test: design points out of range?
X if(any(range(c(.DDO$Xdesignrange, unique(xdesign))) !=
X .DDO$Xdesignrange)) {
X stop(cat("\nDesign points out of Xdesignrange [",
X .DDO$Xdesignrange, "] range\n"))
X }
X#test: too few design points?
X if((this <- length(unique(xdesign))) < length(.DDO$mu)) {
X stop(cat("fewer design points (", this,
X ") than DOF (", length(.DDO$mu), ")\n"))
X }
X#now generate the simulated values
X if(is.null(extend.info$beta)) {
X#We are not reusing previous simulation values (ie,
X#we have not been called by expand()) - so we must make
X#new ones.
X betahatall <- numeric()
X .Random.seed <<- randomseed
X simname <- paste("simulate", .DDO$application, sep = ".")
X lastsim <- list()
X nerr <- 0
X j <- 0
X errmessage <- character()
X simcall <- call(simname, xdesign,
X mu=mu,fam0=fam0,family=family, ...)
X############################BEGIN TIMING OF SIMULATION
X simtime <- oursys.time(while(j < n) {
X lastsim <- eval(simcall)
X#record number of errors.
X nerr <- nerr + as.numeric(lastsim$error)
X if(!lastsim$error)
X betahatall <- rbind(betahatall, lastsim$beta)
X else
X errmessage <- c(errmessage, lastsim$errinfo$message)
X if(nforce & lastsim$error) {
X if(abs(floor(r <- (2 * nerr)/n) - r) < 1e-008)
X cat("\nForcing", n, "simulations; ", nerr, "errors,", j,
X "successes so far.\n")
X }
X else j <- j + 1
X }
X )
X############################END TIMING OF SIMULATION
X #record the arguments upon which this simulation was based
X call0 <- c(call0, lastsim$arg)
X }
X else {
X #we are reusing beta values, so there's
X #almost no simulation time
X simtime <- oursys.time({
X betahatall <- extend.info$beta
X nerr <- extend.info$nerr
X })
X#################
X }
X nstoreflag <- !is.na(match(.DDO$CDname.all, CDname.store))
X CDall <- numeric()
X######################BEGIN TIMING OF EVALUATION
X #if not simulations succeeded, we will return NA
X retNA <- F
X evaltime <- oursys.time({
X if(!is.null(nrow(betahatall))){
X for(npt in 1:nrow(betahatall)) {
X CD <- discrepancyC(betahatall[npt, ], mu,
X roixy = roixy)
X CDall <- rbind(CDall, CD[nstoreflag])
X }
X if(extend.info$DDget)
X DDall <- discrepancyD(CDall,nerr=nerr*(1-droperr))
X else
X DDall <- NULL
X }
X else{
X CDall <- DDall <- NULL
X retNA <- T
X }
X######################END TMING OF EVALUATION
X }
X )
X out <- list(beta = betahatall, CDall = CDall, DDall = DDall,
X rseedend = eval(.Random.seed),
X call = callform, evalcall = call0,
X error = {
X if(nerr)
X list(nerr = nerr, message = {
X if(length(errmessage)){
X ter <- table(errmessage)
X ter[sort.list(-ter)]
X }
X else NULL
X })
X else NULL
X },
X time = list(time.simulation = unlist(simtime),
X time.evaluation = unlist(evaltime)),
X application=.DDO$application)
X if(tag)
X cat("\n")
X if(retNA){
X ret <- NA
X attr(ret,"dump") <- out
X out <- ret
X if(tag)
X cat("NA PRODUCED (NO SUCCESSFUL SIMULATIONS)! ")
X else
X if(extend.info$DDget)
X warning("NA produced - no successful simulations")
X }
X class(out) <- "ddo"
X if(tag)
X cat("FINISHED ", main, "\n (simtime=",
X simtime, "evaltime=", evaltime, ")\n")
X out
X}
X##############################################################################
Xroimake <- function(mu=.DDO$mu,fam0=.DDO$fam0,...){
X #Function to make an object containing complete information
X #about RANGE OF INTEREST (roi), conditional upon mu and fam0.
X #Required input may be in several forms:
X # nwxy - as a complete object which will simply be re-output.
X # name - looks up arguments in .DDO$roimenu corresponding
X # to the name. These arguments are appended to
X # the argument list received by this function.
X # (And, since an x vector is sufficient to generate
X # the entire object, no other arguments need be given
X # if, as is usually the case, the menu item is of the
X # form name=list(x=c(...)) ).
X # x or y - if either x or y vector missing, it will be calculated
X # from the other using mu and fam0. At least one of
X # x or y must be supplied.
X #Optional input includes possible explicit specification of mu and fam0
X #to use for calculation of a missing x or y; and also relative weights:
X # w - relative weights. Is repeated to make vector the
X # same length as x or y. Default is w=1. Is normalized
X # to sum to 1 for output object.
X # mu - mu to use for calculation of missing x or y
X # (default is .DDO$mu)
X # fam0 - fam0 to use for calculation of missing x or y
X # (default is .DDO$fam0)
X #Outputs:
X # name - name argument, or else a constructed name of the
X # form "n points X in [xmin,xmax]".
X # x - roi expressed on X scale
X # y - roi expressed on Y scale
X # w - relative weights of different objects.
X arg <- list(...)
X if(!is.null(arg$nwxy))
X arg <- arg$nwxy
X else{
X #no complete object given. Check first for name.
X #if a name is given - pull down args from menu
X if(!is.null(arg$name))
X arg <- c(arg,.DDO$roimenu[[arg$name]])
X #now see which is missing (of x or y), if any, and replace
X if(is.null(arg$y)){
X if(is.null(arg$x))
X stop("need an x or y input specification")
X else
X arg$y <- yofx(mu,arg$x,fam0)
X }
X else
X if(is.null(arg$x))
X arg$x <- xofy(mu,arg$y,fam0)
X #make weights if needed
X if(is.null(arg$w))
X arg$w <- rep(1,length(arg$x))
X #normalize weights
X arg$w <- arg$w/sum(arg$w)
X #make name if needed
X if(is.null(arg$name))
X arg$name <- paste(length(arg$x)," points X in [",
X (r<-signif(range(arg$x),3))[1],
X ",",r[2],"]",sep="")
X }
X return(name=arg$name,w=arg$w,x=arg$x,y=arg$y)
X }
X##############################################################################
Xshowddo <- function(ddo,...){
X #function to show the design and basic information (n, DDall)
X #from a discrepancy distribution object (ddo).
X DDall <- ddo$DDall
X out <- rbind(c(at(ddo),dimnames(DDall)[[2]]),
X cbind(dimnames(DDall)[[1]],
X array(as.character(signif(DDall,3)),dim(DDall)) ))
X showxddo(ddo,TABinfo=out,...)
X invisible()
X }
X##############################################################################
Xcompareddo <- function(distlist, CDname = .DDO$CDname.all, CL=.95,
X mwspread = c(2, 100), analysis = "dem", coscale = .DDO$coscale,
X DDElist = c("DDmean","DDmedian","DDinterq","DDq95"),
X maxsize = 3, maxrow = 6, tag = date())
X{
X#function to compare several distributions
X#inputs:
X# distlist - list of distribution objects to be compared, or else a
X# class ddoframe object
X# mwspread - spread of comparison-points (see code) across which to
X# do Mann-Whitney statistics.
X# maxsize - largest # of frame rows and columns per page for D, d analyses
X# maxrow - largest number of rows per page for f, F analyses
X# CDname - the names of the discrepancy functions to be used in the
X# comparison.
X# coscale - scaling groups used by analysis f(see gpar.plotdensity()).
X# default: .DDO$coscale
X# analysis - Controls which analyses output.
X# D - show designs in original order
X# d - show designs ranked in order roughly according to # efficiency
X# e - plot efficiencies of DDuse
X# E - plot efficiencies of all DDs
X# f - plot density
X# m - calculate Mann-Whitney statistics
X# tag - what should be stamped in the lower right corner of
X# each plot
X if(!is.null(class(distlist)) && class(distlist)=="ddoframe")
X distlist <- eval(call("[.ddoframe",distlist,virtual=F),
X sys.parent())
X mystamp <- function(tag = tag) NULL
X #this function can be created modelled on stamp()
X#the number of possible analyses:
X anname <- c("D", "d", "f", "m","e","E")
X nanalysis <- length(anname)
X #anuse <- list()
X analysisc <- substring(analysis,1:nchar(analysis),1:nchar(analysis))
X anuse <- as.list(!is.na(match(anname,analysisc)))
X names(anuse) <- anname
X #for(i in 1:nanalysis)
X # anuse[[anname[i]]] <- !is.na(match(anname[i], analysisc)))
X ndist <- length(distlist)
X##
X##Check if we are trying to compare different ranges of interest
X##
X if(ndist < 2)
X stop(cat("Need at least two ddo objects to use compareddo;",
X "\nUse focus() to examine a single object.\n"))
X fcomp <- function(ec)
X list("roixy$x (ROI X values)"=ec$roixy$x,
X "roixy$w (ROI weights)"=ec$roixy$w,
X "roixy$name (ROI name)"=ec$roixy$name,
X "measurement.system"=ec$measurement.system,
X "Xconceivable"=ec$Xconceivable)
X ec1 <- fcomp(distlist[[1]]$evalcall)
X for(i in 2:ndist){
X eci <- fcomp(distlist[[i]]$evalcall)
X for(tp in seq(along=ec1))
X if(length(tp1 <- ec1[[tp]]) != length(tpi <- eci[[tp]])
X || any(tp1!=tpi))
X warning(paste("Incompatible",names(ec1)[tp],
X "between ddo",
X i,"and first ddo."))
X }
X #prepare for plotting
X oldpar <- par("ask","oma")
X on.exit(par(oldpar))
X par(oma=c(0,0,0,0),ask=T)
X#
X#now figure out what discrepancy functions we will plot:
X CDnameOR <- character()
X CDnameAND <- CDname
X for(ddo in distlist) {
X CDnameOR <- unique(c(CDnameOR, ddo$evalcall$CDname.store))
X CDnameAND <- CDnameAND[!is.na(match(CDnameAND, ddo$evalcall$CDname.store))]
X }
X CDnameOR <- CDnameOR[!is.na(match(CDnameOR,CDname))]
X#give them default names if they don't come with names
X maind <- character()
X for(i in 1:ndist) {
X nm <- distlist[[i]]$evalcall$main
X if(is.null(nm))
X nm <- paste("ANALYSIS", i)
X maind <- c(maind, nm)
X }
X numnam <- paste("[#", 1:ndist, "]", sep = "")
X maind <- paste(maind, numnam)
X#initialize framing parameters and save previous framing parameters:
X parmfrow <- list(mfrow=par("mfrow"))
X on.exit(par(parmfrow))
X sizef <- min(ceiling(sqrt(ndist)), maxsize)
X sizef2 <- sizef - (sizef * (sizef - 1) >= ndist)
X nrowf <- min(ndist, maxrow)
X####
X#### BEGIN ANALYSES:
X####
X#Analysis "D": Show them in original order:
X if(anuse$D) {
X par(mfrow = c(sizef, sizef2))
X for(i in 1:ndist)
X showxddo(distlist[[i]], main = maind[i], blockframe = T)
X mystamp()
X }
X#Now prepare for later analyses: Rank
X# Then all subsequent plots will be in the ranked order,
X# according to the variable sortall (below).
X#
X#We will use a ranking scheme from subroutine rankddo
X#CDname2 is that subset of the requested CDnames (in CDname)
X#which are usable for making the comparisons by virtue
X#of being calculated for all the discrepancies
X CDname2 <- CDname[!is.na(match(CDname, CDnameAND))]
X ddor <- rankddo(distlist, CDname = CDname2)
X#
X#Now create the ordering
X#
X sortall <- order(ddor$orank)
X#Analysis "d": plot the designs and tell ranks and summary DDuse.
X if(anuse$d) {
X par(mfrow = c(sizef, sizef2))
X for(ido in sortall) {
X dn <- distlist[[ido]]$evalcall
X atn <- at(distlist[[ido]])
X out1 <- c(atn, "overall", CDname2)
X out2 <- c("DDuse","",signif(distlist[[ido]]$DDall[
X "DDuse",CDname2],5))
X out3 <- c(numnam[ido], as.character(round(
X c(ddor$orank[ido],
X ddor$rank[ido, ]), 2)))
X outall <- cbind(out1, out2, out3)
X gsp <- showxddo(distlist[[ido]], main = maind[ido],
X TABinfo=outall,
X graphsetparreturn = T, blockframe = T)
X }
X mystamp()
X }
X#Analyses "e" and "E": Show efficiencies:
X# "e": DDuse only
X# "E": all DDs
X if(anuse$e | anuse$E){
X DDvec <- c("DDuse",DDElist)[(2-anuse$e):
X (1+anuse$E*length(DDElist))]
X sizee <- min(ceiling(sqrt(length(CDnameOR))),maxsize)
X sizee2 <- sizee - (sizee*(sizee-1)>=length(CDnameOR))
X z <- qnorm(.5*(1+CL))
X #now plot the appropriate DDs:
X for(DDnamu in DDvec){
X par(mfrow=c(sizee,sizee2))
X #set overall margins; later put in names
X par(oma=c(.1+ndist,0,1.4,0))
X #now collect in rows of DDm and DDse for each
X #distribution the appropriate DDmeasures and ses.
X DDm <- DDse <- numeric(0)
X for(i in 1:ndist){
X DDall <- distlist[[i]]$DDall
X CDnamu <- match(CDnameOR,
X distlist[[i]]$evalcall$CDname.store)
X DDrowu <- match(paste(DDnamu,c("",".se"),
X sep=""),dimnames(DDall)[[1]])
X DDm <- rbind(DDm,DDall[DDrowu[1],CDnamu])
X DDse <- rbind(DDse,DDall[DDrowu[2],CDnamu])
X }
X #now do actual plotting:
X for(j in seq(along=CDnameOR)){
X DDmdo <- DDm[,j]
X DDsedo <- DDse[,j]
X #DDL, DDU are vectors of up & low lims
X DDL <- DDmdo - z*DDsedo
X DDU <- DDmdo + z*DDsedo
X #Ddmfin, DDLUfin are vectors of indicators
X #of whether points to be plotted are finite
X DDmfin <- is.finite(DDmdo)
X DDLUfin <- is.finite(DDL)&is.finite(DDU)
X yra <- range(c(DDmdo[DDmfin],
X DDL[DDLUfin],DDU[DDLUfin]))
X if(all(is.na(yra)))
X yra <- c(0,1)
X if(diff(yra)==0)
X yra <- c(0,2*yra[2])
X yra <- range(c(0,yra))
X DDnamlab <- DDnamu
X if(DDnamlab=="DDuse")
X DDnamlab <- paste(attr(DDall,
X "usename")[
X match(CDnameOR[j],
X dimnames(DDall)[[
X 2]])],"(=DDuse)")
X #plot the central points
X nds <- 1:ndist
X matplot(nds,
X DDmdo,pch=15,
X col=2,type="p",ylab="X units",
X xlab=paste("Design Number ([#1]-[#",ndist,"])",sep=""),
X xlim=c(.5,ndist+.5),ylim=yra)
X title(main=paste(CDnameOR[j],"\n[summarized by ",
X DDnamlab,"]",sep=""))
X #plot the upper and lower Conf. limit lines
X matlines(rbind(c(nds,nds)-
X .25,c(nds,nds)+.25),rbind(
X c(DDL,DDU),c(DDL,DDU)),type="l",lty=1,col=4)
X #horizontal line at zero:
X abline(h=0,col=3)
X #plot percentage offsets
X DDmdobest <- min(DDmdo[DDmfin])
X DDmpct <- round(100*(DDmdo/DDmdobest-1),0)
X text(nds+.25,DDmdo,paste(as.character(DDmpct),
X c("","%")[2-is.na(match(DDmdo,DDmdobest))
X ],sep=""), col=2,
X cex=min(1,.20/par("cxy")[1]),adj=0)
X #plot line connecting upper, lower
X if(any(DDLUfin))
X matlines(rbind(nds,nds)[,DDLUfin],
X rbind(DDL,DDU)[,DDLUfin],type="l",
X col=4,lty=1)
X #for nonfinite values plot "?" marks:
X if(any(!DDmfin))
X matpoints(nds[!DDmfin],rep(mean(yra),
X sum(!DDmfin)),type="p",
X pch="?",col=3)
X if(any(!DDLUfin))
X matpoints(c(rbind(nds,nds)[,!DDLUfin]),rep(
X yra+c(1,-1)*par("cxy")[2],
X sum(!DDLUfin)),
X pch="?",col=3)
X #end looping on type of CD
X }
X #end looping on type of DD by writing names
X # in overall margins
X mtext(paste("Percent Loss in Efficiency and",round(CL,2),
X "CIs for DESIGN-level DISCREPANCIES:",
X DDnamu),
X side=3,line=.1,outer=T,cex=1.25)
X for(h in 1:ndist)
X mtext(paste(maind[h]," ",at(distlist[[h]])),
X side=1,line=-.9+h,outer=T)
X }
X #end efficiency analysis
X par(oma=c(0,0,0,0))
X }
X#Analysis "m": Do some Mann-Whitney tests:
X#options are determined by mwspread, which controls the range
X#down from the best (pointwise) test to compare with others
X#(default: mwspread = c(2,100))
X#Under the model points~N(m,s^2) we fit s. Then from the design
X# with the best(smallest) number of points we make all pairwise
X# comparisons where better is <= best+s*mwspread[1],
X# second-best is <= best+s*mwspread[2] points.
X if(anuse$m) {
X ctpoint <- min(ddor$opoints) + mwspread *
X sqrt(var(ddor$opoints)/length(ddor$opoints))
X rrnd <- round(ddor$orank, 2)
X numlike <- sum(ddor$opoints <= ctpoint[1])
X numdo <- sum(ddor$opoints <= ctpoint[2])
X for(i in 1:numlike) {
X ido <- sortall[i]
X dn <- distlist[[ido]]$evalcall
X #first display the design
X out1 <- paste(rrnd[ido], "ranked design", numnam[ido])
X out2 <- CDnameAND
X #make all the same length for prettier output.
X for(i2 in 1:numdo) {
X if(i2 == i) {
X out1n <- if(numdo > 1) " " else
X "other designs can't compete"
X out2n <- rep(" ", length(CDnameAND))
X }
X else {
X mwtst <- mwtest(distlist[[ido]],
X distlist[[sortall[i2]]])
X out1n <- paste("P[ ", rrnd[ido], numnam[ido],
X at(distlist[[ido]]), "<", rrnd[
X sortall[i2]], numnam[sortall[i2]],
X at(distlist[[sortall[i2]]]), "] =",
X sep = "")
X out2n <- as.character(round(mwtst$mw, 4))
X sig <- c(" ", "* ", "**")[1 +
X (mwtst$pval < 0.05) +
X (mwtst$pval < 0.01)]
X out2n <- paste(out2n, sig, sep = "")
X }
X out1 <- c(out1, out1n)
X out2 <- rbind(out2, out2n)
X }
X outall <- cbind(out1, out2)
X showxddo(distlist[[ido]], main = maind[ido],
X TABinfo = outall,
X graphsetparreturn = T)
X mystamp()
X }
X }
X
X#Analysis "f": Show densities
X#gpar.plotdensity() returns a list of one or two vectors, each of same
X#length as the CDname argument to gpar.plotdensity (which we will
X#set equal to CDnameOR so as to plot every discrepancy that
X#appears even once in any object). The CDlim vector is the
X#plot limits for the discrepancy (CD) axis, and the
X#denlim vector, returned if argument maxden=T, gives the
X#plot limits for the density axis.
X if(anuse$f) {
X denlimuse <- rep(0, length(CDnameOR))
X names(denlimuse) <- CDnameOR
X CDlimuse <- denlimuse
X for(i in 1:ndist) {
X plimits <- gpar.plotdensity(distlist[[i]], maxden = T,
X CDname = CDnameOR, coscale = coscale)
X CDlimuse <- pmax(CDlimuse, plimits$CDlim)
X denlimuse <- pmax(denlimuse, plimits$denlim)
X }
X par(mfrow = c(nrowf, length(CDnameOR)))
X for(i in sortall) {
X plotdensity(distlist[[i]], tlim =
X list(CDlim = CDlimuse, denlim = denlimuse),
X blockframe = T, CDname = CDnameOR)
X mtext(maind[i],side=4,line=0)
X }
X mystamp()
X }
X invisible()
X}
X##############################################################################
Xfocus <- function(ddo,...){
X #front end to run three functions for looking in depth
X #into a single ddo:
X showddo(ddo,...)
X plotdensity(ddo,...)
X pairddo(ddo,...)
X plotstrata(ddo,...)
X invisible()
X }
X##############################################################################
Xplotdensity <- function(ddo,
X CDname = ddo$evalcall$CDname.store,
X m = 20,
X qua = c(0.5, 0.95),
X lty = c(4, 3),
X maxdfrow = 6,
X coscale = .DDO$coscale,
X excl = 10,
X tlim = F,
X blockframe = F,
X gparextra=list(),
X ...)
X{
X#function to plot density of discrepancy distribution
X#Inputs
X# ddo - a discrepancy distribution object (generated by ddomake)
X# CDname - a vector of the names of the type of discrepancy
X# function to be plotted. Default: all
X# m - the number of equally spaced points to be calculated
X# for density function if using S function density()
X# qua - the quantiles to be plotted and printed on distribution
X# default: c(.5,.95)
X# numeric vector - use this vector as quantiles
X# lty - the line-types with which to plot the quantiles
X#maxdfrow- controls the maximum number of frame rows in a page of plots
X# coscale- the coscaling groups (see gpar.plotdensity() for details)
X# default: .DDO$coscale. If tlim is given then
X# coscale need not be supplied (ie, will not be used).
X# excl - if 10 (or any pos integer), will exclude all values of
X# discrepancy more than 10*median(discrepancy)
X# default: 10
X# tlim - list of one or two vectors: plotting limits
X# for discrepancy, and perhaps for density.
X# If not suppplied, discrepancy plotting limits set to
X# default: max(3*median discrepancy, max nonexcluded discrepancy)
X# plus .5*median discrepancy
X# as described in gpar.plotdensity()
X#blockframe - unless true, the program will automatically reset par(mfrow)
X# to appropriate dimensions to plot all on one page
X#gparextra - used to pass mfg parameters from plotddo. It is a list of extra
X# graphical parameters to be appended to each plot.
X#outputs: none
X#side effects:
X# generates 1 or 2 plots for each discrepancy function in CDname
X#
X#first test CDname:
X if(any(this <- is.na(match(CDname, .DDO$CDname.all)))) stop({
X cat("\nError: Unknown discrepancy function name(s): ")
X cat(CDname[this], "\nterminating\n")
X }
X )
X atn <- at(ddo)
X ndis <- dim(ddo$CDall)[1]
X n <- nrow(ddo$CDall)
X if(!blockframe) {
X parmfrow <- list(mfrow = par("mfrow"))
X on.exit(parmfrow)
X par(mfrow = c(min(maxdfrow, length(CDname)), 1))
X }
X#now get plotting limits
X if(is.logical(tlim))
X plimits <- gpar.plotdensity(ddo, excl = excl,
X m = m, coscale = coscale,
X CDname = CDname[!is.na(match(
X CDname, ddo$evalcall$CDname.store))])
X else plimits <- tlim
X for(CDnameuse in CDname) {
X #get the extra graphical parameters for this plot:
X if(length(gparextra)){
X do.call("par",c(list(new=T),gparextra[[1]]))
X gparextra <- gparextra[-1]
X }
X #get the plotting limits for this discrepancy type:
X xprange <- c(0, plimits$CDlim[CDnameuse])
X yprange <- c(0, plimits$denlim[CDnameuse])
X if(is.na(match(CDnameuse, ddo$evalcall$CDname.store))) {
X plot(0, 0, xlim = xprange, ylim = c(0, 1), type = "n")
X text(mean(xprange), 0.5, "Not calculated",
X lty = 1, col = 4)
X }
X else {
X CD <- sort(ddo$CDall[, CDnameuse])
X qqout <- getqua(qua, CD, n)
X qoutci <- getquaci(qua, CD, n)
X #Drop excluded values:
X maxuseCD <- excl * median(CD)
X CD <- CD[CD <= maxuseCD] #Now plot the density
X #use the S density() function:
X if(n==1){
X plot(0,xlim=xprange,ylim=yprange,type="n")
X text(mean(xprange),mean(yprange),labels=
X "Can't calculate density for n=1")
X }
X else{
X #case n>1, we can plot densities
X matplot((this <- density(CD, n = m))$x,
X this$y,
X type = "l", lty = 1, col = 4, xlim = xprange,
X ylim = yprange)
X maint <- paste(CDnameuse, " density ",
X paste(atn,";",sep="")," ",attr(
X ddo$DDall, "usename")[match(CDnameuse,
X dimnames(ddo$CDall)[[2]])], "[",
X CDnameuse, "]=", signif(ddo$DDall["DDuse",
X CDnameuse], 4), sep = "")
X mtext(maint,side=3,line=0,cex=pmin(par("cex"),
X diff(par("usr")[1:2])/par("cxy")[1]/nchar(maint)))
X mtext(CDnameuse,side=1,line=0)
X #now put markers at same quantile levels as for dist
X scaleh <- max(this$y)
X#### plot the quantiles:
X for(q in seq(along = qua)) {
X ht <- yprange[1] + diff(yprange) * (0.9 + 0.1 *
X (.5- q%%pmin(length(qua),10)) )
X htphalf <- ht + diff(yprange) * 0.035
X if((!(is.na(qqout[q]))) && (qqout[q] <= maxuseCD)) {
X abline(v = qqout[q], lty = lty[q], col = 1)
X text(qqout[q], htphalf, paste("q(", qua[q], ")=",
X signif(qqout[q], 3), sep = ""),
X cex=pmin(par("cex"),
X diff(yprange)*.07/par("cxy")[2]), lty = 1, col = 1)
X }
X matlines(c(qoutci$qmin[q], qoutci$qmax[q]),
X ht2 <- c(ht, ht), type = "l", lty
X = 1, col = 1)
X if(!qoutci$qokmin[q])
X matlines(c(qoutci$qmin[q], xprange[1]), ht2,
X type = "l", lty = 2, col = 2)
X if(!qoutci$qokmax[q])
X matlines(c(qoutci$qmax[q], xprange[2]),
X ht2, type = "l", lty = 2, col = 2)
X }
X#######end of quantile plotting
X #end of plotting case n>1
X }
X#Done with "else" branch when data is present
X }
X }
X#We're done with all discrepancy names - exit
X invisible()
X}
X##############################################################################
Xplotstrata <-function(ddo, CDname = ddo$evalcall$CDname.store,
X PDname = attr(.DDO$CDname.all, "PDname.all")[match(CDname, .DDO$
X CDname.all)],
X qua = c(0.05, 0.1, 0.25, 0.5, 0.75, 0.9, 0.95),
X xplotlim = {
X r <- range(ddo$evalcall$roixy$x)
X r <- r + (r - mean(r)) * 0.5
X xconr <- range(unlist(xcon))
X c(max(r[1], xconr[1]), min(r[2], xconr[2]))
X },
X xplot = sort(c(ddo$evalcall$roixy$x,
X seq(from = xplotlim[1], to = xplotlim[2], length = 100))),
X xcon = list(absolute=round(c(-1e10,1e10)*max(abs(
X ddo$evalcall$roixy$x)),1),
X truncation=.DDO$Xconceivable),
X ...)
X{
X#function to plot the predicted curve residuals, stratified by
X# discrepancy function quantiles
X#inputs:
X# ddo - the discrepancy distribution object
X# CDname - the names of the discrepancy functions to be output
X# PDname - the names of the pointwise discrepancy to be plotted.
X# PD indicates discrepancy (median bias of Xmeasured)
X# PA indicates aberrance
X# qua - the quantile levels serving as cutpoints between strata
X# xplotlim - the max and min of the x values to be plotted
X# xplot - the vector of X values to be plotted
X# xcon - the range of conceivable X values
X# ... - extra arguments ignored
X#output: none
X#side effect:
X# For each discrepancy in CDname,
X#separates calibration curves into strata, then for each strata plots
X#X against the median of X measured (solid lines). Also on graph are plotted
X#bands showing the asymptotic standard error of X measured (broken line).
X lpt <- length(xplot)
X ec <- ddo$evalcall
X#get the asymptotic standard error of X measured:
X intrin <- intrinsic(xcon = xcon, mu = ec$mu, roixy = ec$roixy,
X xplotlim = xplotlim, xplot = xplot,
X sigma = ec$sigma, fam0 = ec$fam0, doplot = F)
X parreset <- par("ask", "mfrow","oma")
X on.exit(par(parreset))
X par(oma=c(1.1,0,2.1,0),ask=T)
X#
X#Now begin looping on type of discrepancy function:
X#
X for(di in seq(along = CDname)) {
X CDnameuse <- CDname[di]
X PDnameuse <- PDname[di]
X Dall <- ddo$CDall[, CDnameuse]
X maxD <- max(Dall)
X #Decide upon the vertical limits we will force on every plot:
X #look at the top few
X #branch by type of pointwise aberrance
X yra <- switch(PDnameuse,
X PD = c(-1, 1) * max(ddo$DDall["DDmax", "CDinf"],
X intrin$Xse),
X PA = c(0, 1),
X stop(paste("Improper pointwise discrepancy name (",
X PDnameuse, ")\n", sep = "")))
X #qneed is the quantiles levels wanted
X# (so it is in units 0 to 1)
X#qa is the estimated values of discrepancey function producing them
X# (so it is in discrepancy units)
X qneed <- sort(qua)
X qa <- quantile(Dall, qneed)
X#use only those quantiles that have values, ie,
X# are not too large or too small; and use corresponding
X# cumulative probability levels:
X qado <- qa[!is.na(qa)]
X s1do <- qneed[!is.na(qa)]
X#the number of inter-quantile strata including 0, infinity
X nmplot <- length(qado) + 1
X#the number of squares per side to get enough graphs:
X szplot <- ceiling(sqrt(nmplot))
X szplot2 <- szplot - (szplot * (szplot - 1) >= nmplot)
X #make flag to ask for overall header and footer:
X oflag <- T
X par(mfrow = c(szplot, szplot2))
X #
X #now loop on strata:
X #
X for(str in 1:nmplot) {
X Dmin <- c(0, qado)[str] #the lower limit of discrepancy
X Dmax <- c(qado, maxD)[str] #the upper limit of discrepancy
X Qmin <- c(0, s1do)[str] #the lower limit of cum prob
X Qmax <- c(s1do, 1)[str] #the upper limit of cum prob
X#determine which betas to pull out for this stratum;
X stwhich <- (Dmin < Dall) & (Dall <= Dmax)
X if(!any(stwhich)) {
X plot(x = 1, y = 1, type = "n", xlim = xplotlim,
X ylim = yra, xlab="X", ylab=PDnameuse,
X main = paste(signif(Dmin, 3), "<", CDnameuse,
X "< ", signif(Dmax, 3)), sub = paste(Qmin,
X "< quantile <=", Qmax))
X }
X else {
X#We have calibrations to plot. First make plot
X stbeta <- ddo$beta[stwhich, , drop = F]
X intrinplot <- switch(PDnameuse,
X PD = cbind("-SDasymptotic"=- intrin$Xse,
X "+SDasymptotic"= intrin$Xse),
X PA = sapply(intrin$PA,I))
X #make overall footer if needed
X matplot(xplot, intrinplot, type = "l", lty = 2, col = 2,
X xlim = xplotlim, ylim= yra,xlab = "X",ylab = PDnameuse,
X main = paste(signif(Dmin, 3), "<",
X CDnameuse, "<= ", signif(Dmax, 3), "\n(",
X nrow(stbeta), " calibrations)", sep= ""),
X sub = paste(Qmin, "< quantile <=", Qmax))
X if(oflag){
X #plot overall titles if needed
X oflag <- F
X mtext(ddo$evalcall$main,3,0,outer=T,
X cex=2*par("mex"))
X PDnameintrin <- switch(PDnameuse,
X PD="SD(Xmeasured)",
X PA="pointwise aberrance PA")
X mtext(paste("Dotted lines are intrinsic ",
X PDnameintrin," (",paste(
X dimnames(intrinplot)[[2]],collapse=", "),
X ")",sep=""),1,0,outer=T)
X }
X #augument the plot with the x-axis (==original curve):
X abhoriz <- switch(PDnameuse,
X PD = 0,
X PA = c(0, 1))
X abline(h = abhoriz, lty = 1, col = 1)
X abline(v = range(ec$roixy$x), lty = 3, col = 4)
X#now dothe individual "curves":
X for(ci in 1:nrow(stbeta)) {
X xm <- xofydist(y = yofx(ec$mu, xplot, family = ec$fam0),
X betahat = stbeta[ci,], family = ec$family,
X sigma = ec$sigma,
X Xconceivable = ec$Xconceivable,
X measurement.system = ec$measurement.system)
X xD <- switch(PDnameuse,
X PD = xm$Xmedian - xplot,
X PA = xm$PA)
X matlines(xplot, xD, type = "l", lty = 1, col = 4)
X }
X#finish with the "else" for not zero curves:
X }
X }
X#finish looping on type of discrepancy:
X }
X#finish function:
X invisible()
X}
X##############################################################################
Xpairddo <- function(ddo, x = ddo$CDall, label0 = dimnames(x)[[2]],
X labels = paste(ddo$evalcall$main,"\n",label0),
X CAintrinsic = max(
X do.call("xofydist", c(list(
X y = range(ddo$evalcall$roixy$y),
X betahat = ddo$evalcall$mu,
X family = ddo$evalcall$fam0),
X ddo$evalcall[c("sigma", "Xconceivable",
X "measurement.system")]))$PA),
X panel = points,
X invert = TRUE,
X ...) {
X#function for plotting pairs of discrepancies
X#created by modifying code for pairs.default in Splus
X#input:
X# ddo - a discrepancy distribution object
X# x - a matrix of discrepancy values (usually obtained from ddo).
X# label0 - the names of the discrepancies for use in matching with
X# name of intrinsic aberrance (usually obtained from ddo)
X# labels - the plotting names of the discrepancies (usually obtained from
X# ddo).
X# CAintrinsic - the intrinsic aberrance. It is calculated from
X# the ddo usually.
X# other arguments
X# - not investigated when constructing this routine.
X#output:
X# none
X#side effect
X# Creates plots of pairs of discrepancies much like in pairs()
X#except that a dotted line at CA=CAintrinsic is plotted.
X doaxis <- function(which, dolabel = T)
X axis(which, outer = T, line = -0.5, labels = dolabel)
X setup <- function(x, y, ...)
X .S(plot(range(x[!is.na(x)]), range(y[!is.na(y)]), type = "n", axes = F,
X ...), "plot")
X x <- as.matrix(x)
X if(is.character(panel))
X panel <- get(panel, mode = "function")
X n <- ncol(x)
X oldpar <- par("oma", "mar", "cex", "tck", "mfg", "mgp", "mex", "mfrow")
X oldcex <- par("cex")
X CEX <- oldcex * max(7.7/(2 * n + 3), 0.6)
X par(mfrow = c(n, n), mgp = c(2, 0.8, 0),
X oma = rep(3, 4), mar = rep(0.5, 4), tck = -0.03/n)
X on.exit({
X par(oldpar)
X }
X )
X par(cex = CEX)
X if(length(labels) < n)
X labels <- paste(deparse(substitute(x)), "[,", 1:n, "]", sep =
X "")
X if(par("pty") == "s") {
X dif <- diff(par("fin"))/2
X if(dif > 0)
X par(omi = c(dif * n, 0, dif * n, 0) + par("omi"))
X else par(omi = c(0, ( - dif) * n, 0, ( - dif) * n) + par("omi")
X )
X }
X order <- if(invert) 1:n else n:1
X for(i in order) {
X for(j in 1:n) {
X setup(as.vector(x[, j]), as.vector(x[, i]), ...)
X box()
X if(i == 1)
X doaxis(3, j %% 2 == 0)
X if(i == n)
X doaxis(1, j %% 2 == 1)
X if(j == 1)
X doaxis(2, i %% 2 == 0)
X if(j == n)
X doaxis(4, i %% 2 == 1)
X if(i != j) {
X panel(as.vector(x[, j]), as.vector(x[, i]), ...
X )
X #The following if and abline statements were
X #introduced to pairddo:
X if(label0[i] == "CA")
X abline(h = CAintrinsic, lty = 2, col = 2)
X if(label0[j] == "CA")
X abline(v = CAintrinsic, lty = 2, col = 2)
X }
X else {
X par(usr = c(0, 1, 0, 1))
X text(0.5, 0.5, labels[i], cex = 1.5 * CEX)
X }
X }
X }
X invisible()
X}
X##############################################################################
Xmwtest <- function(dist1,dist2,CDname=.DDO$CDname.store,droperr=F){
X #Function to do Mann-Whitney tests for pair of discrepancy distribution objs.
X #inputs:
X #dist1,dist2 - the two discrepancy distributions
X # dname - the name of the discrepancy functions to compare
X # droperr - if True, ignore the fact that some values did not
X # converge and thus may be considered to have had
X # infinite discrepancy
X #output: list with entries
X # mw - vector of length=length(CDname) with each entry being
X # P[X= 2)
X do.call("compareddo", c(distlist = list(ddolist), args))
X invisible()
X}
X##############################################################################
Xmake <- function(prefix= .,makelist= .,n=.DDO$n,override=F,
X where=1,...){
X #Function to make a collection of ddos from a list of design specifications
X #inputs:
X # prefix - the prefix for the name of the designs which are created
X # makelist - a list of design specifications. Each design should be in form
X # title = xdesign or
X # title = list(xdesign, argname1=arg1,...argnamek=argk)
X # where args are arguments for ddomake
X # n - the number of simulations to perform
X # override - if FALSE, function refuses to overwrite preexisting objects
X # where - the data base in which to write the objects when created.
X #To add:
X #check for enough space in memory, if possible
X invisible()
X leave <- F
X if(missing(prefix) && !exists("prefix",where=where)){
X cat("prefix missing; set prefix for created ddo's with",
X "\n\tprefix <- \"desired_prefix\"",
X "\nor enter as argument with",
X "\n\tmake(prefix=\"desired_prefix\")\n")
X leave <- T
X }
X else
X if(missing(prefix))
X prefix <- get("prefix",where=where)
X if(missing(makelist) && !exists("makelist",where=where)){
X cat("makelist missing; specify design list with",
X "\n\tmakelist <- list(name1=design1,name2=design2,...)",
X "\n\or enter as argument with",
X "\n\tmake(makelist=desired_makelist)\n")
X leave <- T
X }
X else
X if(missing(makelist))
X makelist <- get("makelist",where=where)
X if(leave)
X return()
X nd <- length(makelist)
X nmall <<- c(paste(prefix,1:nd,sep=""))
X #errorcheck prior existence of objects:
X if(!override){
X obj <- objects(where=where)
X nmallp <- c(prefix,nmall)
X already <- nmallp[!is.na(match(nmallp,obj))]
X if(length(already)){
X cat("Objects with prefix",prefix,"already exist:\n",
X already,"\nUse make(override=T) to override\n")
X return()
X }
X }
X cat(paste("Making ddo class objects ",nmall[1]," through ",
X nmall[length(nmall)]," with n=",n,
X ";\n Accessable through class ddoframe (virtual) object ",
X prefix,"[]:\n",sep=""))
X lnewnam <- list()
X for(di in 1:nd){
X nmi <- paste(prefix,di,sep="")
X xd <- darg <- makelist[[di]]
X args <- list(main=names(makelist)[di],n=n,...)
X assign(nmi,NA,where=where)
X if(!is.numeric(xd)){
X if(!is.list(darg)){
X cat(paste("Error: element ",di,
X " of design list must be numeric or",
X " list.\n",nmi," given value NA\n",
X sep=""))
X next
X }
X xd <- darg[[1]]
X args <- c(args,darg[-1])
X }
X if(!is.numeric(xd) || !is.vector(xd)){
X cat(paste("Error: element ",di,
X " of design list must specify a design",
X "\n\twhich is a numeric vector\n",nmi,
X " give value NA\n",sep=""))
X next
X }
X ri <- do.call("ddomake",c(args,list(xdesign=xd,tag=F)))
X assign(nmi,ri,where=where)
X cat(nmi,"made on",date(),"\n")
X lnewnam <- c(lnewnam,nmi)
X }
X fmake <- do.call("c.ddoframe",lnewnam)
X assign(prefix, fmake,where=where)
X cat(prefix,"made on",unclass(fmake)$maketime,"\n")
X }
X##############################################################################
X#Design canditate management functions:
X
X# names.ddo
X# names.ddoframe
X# elements
X# dimnames.ddoframe
X# length.ddoframe
X# dim.ddoframe
X# names<-.ddoframe
X# c.ddo
X# c.ddoframe
X# [.ddo
X# [.ddoframe
X# [<-.ddoframe
X# print.ddoframe
X##############################################################################
Xnames.ddo <- function(dd){
X #names all ddos should include:
X alln <- c("beta","CDall","DDall","rseedend","call",
X "roix","roiy","roiw","roiname")
X #names particular to the application:
X apn <- names(dd$evalcall)
X apn <- apn[apn!="roixy"]
X #names if error occurred:
X ern <- names(dd$error)
X #names of time for creation:
X tin <- names(dd$time)
X sort(c(apn,alln,tin,ern))
X }
X##############################################################################
Xnames.ddoframe <- function(ddf)
X unclass(ddf)$row.names
X##############################################################################
Xelements <- function(ddf){
X if(is.null(class) || class(ddf)!="ddoframe")
X stop("requires class ddoframe object")
X unclass(ddf)$col.names
X }
X##############################################################################
Xdimnames.ddoframe <- function(ddf)
X list(unclass(ddf)$row.names,unclass(ddf)$col.names)
X##############################################################################
Xlength.ddoframe <- function(ddf)
X length(unclass(ddf)$row.names)
X##############################################################################
Xdim.ddoframe <- function(ddf) c(length(unclass(ddf)$row.names),length(unclass(ddf)$col.names))
X##############################################################################
X"names<-.ddoframe" <- function(ddf,nam){
X if(length(ddf)!=length(nam))
X stop(paste("New names must be proper length (",
X length(ddf),") but instead are length ",
X length(nam),sep=""))
X out <- unclass(ddf)
X out$row.names <- nam
X class(out) <- "ddoframe"
X return(out)
X }
X##############################################################################
Xc.ddo <- function(...){
X #redirect to c.ddoframe
X sc <- sys.call()
X sc[[1]] <- as.name("c.ddoframe")
X eval(sc,sys.parent())
X }
X##############################################################################
Xc.ddoframe <- function(...,application=.DDO$application){
X #function to combine ddo and ddoframe class objects
X # into a ddoframe class object
X #get the call names of the arguments; these will
X #be needed for construction of any ddo
X caln <- match.call()
X mode(caln) <- "list"
X caln <- caln[-1]
X if(!is.null(names(caln)))
X caln <- caln[names(caln)!="application"]
X #now go through objects one by one:
X out <- list(row.names=character(),col.names=character())
X for(ddi in seq(along=caln)){
X if(is.character(caln[[ddi]]) && length(caln[[ddi]])==1)
X caln[[ddi]] <- as.name(caln[[ddi]])
X dd <- eval(caln[[ddi]],sys.parent())
X #transform dd to be an unclassed virtual ddoframe
X ddfvu <- switch(
X if(is.null(class(dd))) "NONE"
X else class(dd),
X ddo={ if(!is.name(caln[[ddi]]))
X stop("class ddo arguments must be in form of name")
X list(row.names=as.character(caln[[ddi]]),
X col.names=names(dd),
X application=dd$application)},
X ddoframe=unclass(dd),
X stop("All arguments must be of class ddo or ddoframe")
X )
X if(ddfvu$application != application)
X stop(paste("All objects must be for application",
X application,"but object",ddi,
X "is for application",ddfvu$application))
X out$row.names <- c(out$row.names,ddfvu$row.names)
X out$col.names <- c(out$col.names,ddfvu$col.names)
X }
X if(any(duplicated(out$row.names)))
X stop(paste("Duplicated design names:\n",
X out$row.names[duplicated(out$row.names)]))
X out$col.names <- unique(out$col.names)
X out$maketime <- date()
X out$application <- application
X class(out) <- "ddoframe"
X return(out)
X }
X##############################################################################
X"[.ddo" <- function(dd,elements,drop=T){
X #Function to extract elements from a ddo.
X out <- list()
X for(el in elements){
X oe <- dd$evalcall[[el]]
X if(!is.na(match(el,
X c("beta","CDall","DDall","rseedend","call"))))
X oe <- dd[[el]]
X if(!is.na(mt <- match(el,c("roix","roiy","roiw","roiname"))))
X oe <- dd$evalcall$roixy[[c("x","y","w","name")[mt]]]
X if(!is.na(match(el,c("nerr","message")))){
X oe <- dd$error[[el]]
X if(is.null(oe))
X oe <- list(nerr=0,message=numeric()
X )[[el]]
X }
X if(!is.na(match(el,c("time.simulation","time.evaluation"))))
X oe <- dd$time[[el]]
X out <- c(out,list(oe))
X }
X if(drop & length(elements)==1)
X return(out[[1]])
X names(out) <- elements
X return(out)
X }
X##############################################################################
X"[.ddoframe" <- function(ddf,ddowhich=ddowhich.default,
X elements=character(),simplify=simplify.default,
X virtual=missing.elements,drop=drop.default){
X #Function to extract elements from a ddoframe class object.
X #Returns either (I)
X # (i) a ("virtual") object of class ddoframe which takes
X # very little space,
X # (ii) a (null class) list of class ddo objects
X # (iii) a (null class) list of lists of collections of ddo # elements
X # ddf - a ddoframe
X # ddowhich - which ddo objects to extract from the ddoframe.
X # Can be either logical, numeric, or character.
X # Default is to extract all of them
X # elements - which elements to extract from each ddo object.
X # Default is to extract ENTIRE ddo object
X # without simplifying at all
X # (and leaving as class ddo)
X # drop - if TRUE and there is only one ddo object, return
X # the ddo rather than list(ddo)
X # simplify - if TRUE and elements are being extracted, collapse
X # elements of different ddo objects together.
X # For example output will include a vector of
X # values for n.
X # virtual - if TRUE returns a class ddoframe rather than
X # a simple list. If elements are specified, then
X # it would not be possible to return a ddoframe
X # object, so virtual is ignored.
X # ddowhich.default - this is inserted as a convenience to keep the
X # default definition in one place. The problem is
X # that there may be two ways to insert the default
X # for ddowhich:
X # if ddoframe[] then ddowhich assumes default value
X # if ddoframe[,] then ddowhich is mode missing
X #outputs:
X # If no elements are specified, return a list of ddo objects
X # (with names corresponding to the objects)
X # (if drop = T and there is only one object specified,
X # then returns that object itself, not list)
X # If elements are specified, returns a list of length equal to
X # the number of elements specified. Each of these
X # list members will in turn be a collection of the
X # corresponding elements from across all ddos.
X # If simplify=F this collection will be in the form
X # if a list of length (number of ddo objects);
X # If simplify=T this collection will be in simplified
X # form if possible (according to function sapply()),
X # ie, all scalars will be collapsed from a list of
X # scalars to a vector, etc.
X ###WARNING: Horrendous (aren't they all?) bug found in the
X ###internal code in function "[": It mismatches arguments
X ###unless the argument drop is made to be the last argument!!!
X ###Even so, further mismatching seems to take place; rather than
X ###substituting in default values for missing arguments, the evaluator
X ###will often instead give a value which is present in the first
X ###entry in sys.frame(). Thus I have hand-substituted all default
X ###values for backup.
X drop.default <- T
X if(missing(drop))
X drop <- drop.default
X simplify.default <- drop
X if(missing(simplify))
X simplify <- simplify.default
X ddfu <- unclass(ddf)
X #get the names of the ddo objects we will access
X ddowhich.default <- rep(T,length(ddfu$row.names))
X if(missing(ddowhich))
X ddowhich <- ddowhich.default
X ddon <- switch(mode(ddowhich),
X numeric=ddfu$row.names[ddowhich],
X logical=ddfu$row.names[ddowhich],
X character={
X nom <- is.na(match(ddowhich,ddfu$row.names))
X if(any(nom))
X warning(paste(c(
X "reference to nonexistent ddo in ddoframe:",
X ddowhich[nom]),collapse=" "))
X ddowhich[!nom]}
X )
X #get names of elements we will access, and a flag for zero
X if(missing(elements))
X elements <- character()
X missing.elements <- length(elements)==0
X #if no elements are specified, and object is virtual,
X #simply use the c.ddo
X if(missing(virtual))
X virtual <- missing.elements
X if(virtual){
X if(missing.elements){
X doc <- c(list(as.name("c.ddoframe")),ddon)
X mode(doc) <- "call"
X return(eval(doc,sys.parent()))
X }
X warning("Can't specify virtual object when extracting elements")
X }
X #object is not virtual, so assemble it into nvobj
X nvobj <- list()
X for(dn in ddon)
X nvobj <- c(nvobj,list(get(as.name(dn))))
X
X #if no elements specified,
X #simply send back the assembled (nonvirtual) object,
X #dropping if necessary
X if(missing.elements){
X if(drop & length(nvobj)==1)
X return(nvobj[[1]])
X names(nvobj) <- ddon
X return(nvobj)
X }
X #elements are extracted. Go through elements one by one
X out <- list()
X for(enm in elements){
X if(!simplify | enm=="call"){
X #We need not simplify. So just use lapply.
X out <- c(out,list(lapply(nvobj,"[",elements=enm)))
X next
X }
X #We must simplify.
X #How we do it depends upon the element
X if(enm=="message"){
X #treat message as a special case
X outelem <- lapply(nvobj,"[",elements=enm)
X if(simplify){
X outelem <- unlist(outelem)
X outelem <- rev(sort(tapply(outelem,
X names(outelem),sum)))
X }
X out <- c(out,list(outelem))
X next
X }
X #for all variables besides message and call
X #use the power of sapply().
X outelem <- sapply(nvobj,"[",elements=enm)
X if(!is.list(outelem)){
X #must reshape - get an exemplar in exemp
X #We want to preserve:
X # the dimensions and
X # the dimension names, if they exist.
X # the miscellaneous attributes of the object
X # dim(exemp) preserves dimensions
X # exnm will contain dimnames
X # exat will contain attributes
X ###DIMENSIONS
X exemp <- as.array(do.call("[.ddo",list(nvobj[[1]],
X elements=enm)))
X #In Windows, the as.array
X #Command leaves in a names
X #Attribute, wrongly
X attr(exemp,"names") <- NULL
X ###DIMENSION NAMES
X exnm <- dimnames(exemp)
X #Create names if needed
X if(is.null(exnm))
X exnm <- lapply(as.list(dim(exemp)),function(a)1:a)
X #Now COMBINE together into new array of combined
X #dimensions (this is object outelem).
X outelem <- array(outelem,
X dim=c(dim(exemp),length(ddon)),
X dimnames=c(exnm,list(ddon)) )
X #preserve attributes
X ###ATTRIBUTES
X exat <- attributes(exemp)
X exat <- exat[is.na(match( names(exat),
X names( attributes(outelem)) ))]
X #exat now contains all attributes of
X #the exemplar which were lost in the
X #sapply function
X attributes(outelem) <- c(attributes(outelem),exat)
X ###Final step: drop dimensions of one, and turn into scalar or
X #vector rather than array if possible
X #(by dropping dim altogether)
X dn1 <- dim(outelem)!=1
X ao <- attributes(outelem)
X if(any(dn1)){
X ao$dim <- ao$dim[dn1]
X ao$dimnames <- ao$dimnames[dn1]
X if(length(ao$dim)==1){
X ao$names <- ao$dimnames[[1]]
X ao$dim <- ao$dimnames <- NULL
X }
X }
X else
X ao$dim <- ao$dimnames <- NULL
X attributes(outelem) <- ao
X }
X #store the result
X out <- c(out,list(outelem))
X #Finish looping on element (variable enm)
X }
X if(drop & length(elements)==1)
X return(out[[1]])
X names(out) <- elements
X return(out)
X }
X##############################################################################
X"[<-.ddoframe" <- function(ddf,ddowhich,value){
X #function for replacing elements in a ddoframe
X ddfu <- unclass(ddf)
X oldnames <- ddf$row.names
X repnames <- switch(mode(ddowhich),
X numeric=oldnames[ddowhich],
X character=ddowhich,
X logical=oldnames[ddowhich])
X err <- "replacement value must be name of ddo or evaluate to class ddo object"
X valuenames <- switch(class(value),
X ddo={mc <- match.call()[["value"]]
X if(!is.name(mc))
X stop(err)
X as.character(mc)},
X ddoframe=unclass(value)$row.names,
X stop(err)
X )
X if(length(repnames)!=length(valuenames))
X stop(paste("replacement ddoframe is length",length(valuenames),
X "but should be",length(repnames)))
X oldnames[match(repnames,oldnames)] <- valuenames
X do.call("c.ddoframe",as.list(oldnames))
X }
X##############################################################################
Xprint.ddoframe <- function(ddf,width=.Options$width,sep=" "){
X #function to print ddoframe object
X flw <- function(v,w,sep){
X ncl <- floor(w/(wn<-max(nchar(v)+nchar(sep))))
X if(ncl==0)
X return(paste(v,collapse="\n"))
X v <- paste(sep,format(v),sep="")
X v <- c(v,rep("",(-length(v))%%ncl))
X v <- matrix(v,nrow=ncl)
X v <- rbind(v,rep("\n",ncol(v)))
X paste(c(v),collapse="")
X }
X cat(paste("Application=",unclass(ddf)$application,": ",sep=""),
X length(ddf),"ddo objects,",dim(ddf)[2],"elements\n")
X cat(paste("Objects:\n",flw(names(ddf),width,sep),sep=""))
X cat(paste("Elements:\n",flw(dimnames(ddf)[[2]],width,sep),sep=""))
X cat("Time created:",unclass(ddf)$maketime,"\n")
X invisible()
X }
X##############################################################################
Xplotddo <- function(...)
X{
X#function to plot a design and its discrepancy distributions
X#can supply it with arguments used by showxddo or plotdensity
X argls <- list(...) #The default is that we don't block reframing
X if(is.null(argls$blockframe))
X nowblock <- F
X else nowblock <- argls$blockframe
X argls$blockframe <- T
X showddopar <- plotdenpar <- list()
X if(!nowblock) {
X argspd <- amatch(plotdensity, call("plotdensity"))
X if(is.null(argls$CDname))
X argls$CDname <- eval(argspd$CDname[[2]],
X local = list(ddo = argls[[1]]))
X nden <- length(argls$CDname)
X if(nden<=5)
X sz <- sz2 <- 3
X else{
X sz <- 1+floor(nden/2)
X sz2 <- nden-sz
X }
X #overall dimensions will be sz rows by sz2 cols;
X #The design plot will take up as much as possible.
X sdes <- sz-(nden>0)
X sdes2 <- sz2-(nden>sz)
X showddopar <- list(omd=c(0,sdes2/sz2,1-sdes/sz,1,),mfrow=c(1,1))
X plotdenpar <- list()
X for(i in 1:sz2)
X plotdenpar <- c(plotdenpar,list(list(mfg=c(sz,i,sz,sz2))))
X for(i in seq(sz-1,1,by=-1))
X plotdenpar <- c(plotdenpar,list(list(mfg=c(i,sz2,sz,sz2))))
X plotdenpar <- plotdenpar[1:nden]
X }
X if(!nowblock){
X parlist <- par("omd","mfrow","mar")
X on.exit(parlist)
X do.call("par",showddopar)
X }
X do.call("showddo", argls)
X if(!nowblock)
X par(omd=c(0,1,0,1), mar=c(2.1,2.1,2.1,1),new=T)
X do.call("plotdensity", c(argls,list(gparextra=plotdenpar)))
X invisible()
X}
X##############################################################################
Xprint.DDO <- function(DINIT=.DDO,menu=F){
X #Function to print the default variables (contained in .DDO)
X #input:
X # DINIT contains the initial default list
X # menu flag indicating whether menus should be printed
X #Calls upon application methods (eg, print.DDO.[application])
X #To print application-specific default variables.
X #Variables to be printed:
X #
X #STATE OF NATURE:
X # mu imagined "true" calibration
X #SIMULATION CONTROL:
X # n number of repetitions
X # fam0 family of functions for data generation
X # family (approximating) family of functions
X # for data analysis
X #EVALUATION CONTROL:
X # Xconceivable Xconceivable range controlling
X # aberrance cutoff points
X # measurement.system "Complete measurement system" to be used.
X # Must be either (1) repeat measurements
X # or (2) truncation.
X # roi Range of Interest
X # CDname.store Names of Calibration-Level Discrepancy functions
X # which are stored.
X # droperr Whether nonconvergent simulations are retained
X # in calculation of quantile-based
X # Design-level discrepancies (eg, DDmedian)
X #GRAPHICS CONTROL:
X # coscale Groups of discrepancies to be scaled together.
X # Xplotrange X range for plotting (showx, etc.).
X #MENUS, ETC.:
X # roimenu Named sufficient (contains either x or y)
X # roi vectors accessable from roi.
X # coscalemenu Named coscaling lists accessable from coscale.
X # CDname.all Names of all legitimate Calibration-Level
X # Discrepancy functions.
X cat("\n\nDEFAULT PARAMETERS FOR APPLICATION ",DINIT$application,":\n\n")
X cat("STATE OF NATURE (UNIVERSAL):\n")
X cat("mu = ",DINIT$mu,"(imagined true parameters)\n")
X cat("SIMULATION CONTROL (UNIVERSAL):\n")
X cat("n = ",DINIT$n," (number of repetitions) nforce =",
X DINIT$nforce,"(force even if nonconvergence?)\n")
X if(is.null(DINIT$fam0)&is.null(DINIT$family))
X cat("fam0=family=NULL (generating and analysis families both default)\n")
X else{
X fam0 <- family <- "NULL [default]"
X if(!is.null(DINIT$fam0)) fam0 <- DINIT$fam0
X if(!is.null(DINIT$family)) family <- DINIT$family
X cat("fam0 =",fam0," (generating family)\nfamily =",
X family," (approximating family)\n")
X }
X cat("EVALUATON CONTROL (UNIVERSAL):\n")
X cat("Xconceivable = [",DINIT$Xconceivable[1],",",DINIT$Xconceivable[2],
X "] (conceivable X range; outside is aberrant)\n")
X cat("measurement.system = ",DINIT$measurement.system,
X " (response to aberrant measurements)\n")
X cat("droperr = ",DINIT$droperr," (omit errors from DD's?\n")
X cat("roi = ")
X menucat <- function(nm,nmlist,lookup=is.character(nm)){
X #Function which looks up named items in menu and prints name
X if(!lookup){
X mout <- nm
X cat("(unnamed)")
X }
X else{
X mout <- nmlist[[nm]]
X cat(nm)
X if(is.null(mout))
X cat(" (not found in menu!)")
X }
X mout
X }
X ra <- menucat(DINIT$roi,DINIT$roimenu)
X cat(" (Range of Interest)\n")
X this <- rbind(x=ra$x,y=ra$y,w=ra$w)
X for(i in 1:nrow(this))
X cat(" ",format(this)[i,],dimnames(this)[[1]][i],"\n")
X cat("CDname.store (Calibration-Level Discrepancies Stored):\n ",
X DINIT$CDname.store,"\n")
X cat("GRAPHICS CONTROL (UNIVERSAL):\n")
X cat("coscale (groups of CDs scaled together on plots) = ")
X coscale <- menucat(DINIT$coscale,DINIT$coscalemenu)
X cat("\n")
X print(coscale)
X cat("Xplotrange = ",DINIT$Xplotrange,"\n")
X if(menu){
X cat("MENUS, ETC. (UNIVERSAL):\n")
X cat("roimenu (menu of possible ranges of interest) =\n")
X print(DINIT$roimenu)
X cat("coscalemenu (menu of possible coscaling group choices) =\n")
X print(DINIT$coscalemenu)
X cat("CDname.all (all recognized Calibration Discrepancy Names) =\n",
X DINIT$CDname.all)
X }
X #now do the application-specific variables
X apname <- paste("print.DDO",DINIT$application,sep=".")
X if(exists(apname)){
X cat("\nAPPLICATION-SPECIFIC (",DINIT$application,") DEFAULT VALUES:\n")
X do.call(apname,list(DINIT=DINIT,menu=menu))
X }
X else
X cat("")
X cat("\n")
X invisible()
X }
X##############################################################################
Xxofydist.default <- function(y, betahat, family, sigma, Xconceivable,
X measurement.system)
X{
X#function to describe what happens when a measurement system is applied
X#to noisy measurements around a vector of Y values. Assuming that we
X#repeatedly measure these same Y values under Y-noise specified by sigma,
X#what is the distribution of Xmeasured that will be returned by the
X#"complete" measurement system?
X#inputs:
X# y - vector of "true" Y values
X# betahat - how Xmeasured is backcalculated from Ymeasured
X# family - family for backcalculation
X# xcon - conceivable X range (used by measurement.system=truncation)
X# sigma - the noise level (SD of Ymeasured | Y)
X# measurement.system -
X# - The "complete" measurement system used [see paper].
X# must be either
X# truncation or repeat.measures
X#outputs:
X# Xmedian - vector of medians of Xmeasured
X# PA - vector of expected number of aberrant measurements
X# produced in the course of measurement
X Yconceivable <- yofx(betahat, Xconceivable, family = family)
X Ymedian <- switch(measurement.system,
X repeat.measures = roixy$y + sigma * qnorm(0.5 * (pnorm((Yconceivable[1] - y)/sigma) + pnorm((
X Yconceivable[2] - y)/sigma))),
X truncation = pmax(pmin(y, Yconceivable[2]), Yconceivable[1]))
X Xmedian <- xofy(betahat, Ymedian, family = family)
X Ycon <- range(Yconceivable)
X pab <- pnorm((Ycon[1] - y)/sigma) + pnorm((y - Ycon[2])/sigma)
X enumab <- switch(measurement.system,
X repeat.measures = pab/(1 - pab),
X truncation = pab)
X return(Xmedian, PA = enumab)
X}
X##############################################################################
Xlindim <- function(info){
X #function to return dimensions of a line of text
X #which may be seperated by many newline characters
X infolin <- infolen <- 0
X if(!(is.null(info) || length(info)==0 || info=="" )){
X infolin <- 1+sum(infosep <- "\n"==substring(info,
X this<-1:nchar(info), this))
X this <- c(0,seq(along=infosep)[infosep],1+length(infosep))
X infolen <- max(diff(this))-1
X }
X c(infolin,infolen)
X }
X##############################################################################
Xchdim <- function(tabl,nsep=2,rch=1.25){
X nc <- ncol(tabl)
X nr <- nrow(tabl)
X nch <- array(nchar(c(tabl)),dim(tabl))
X #the number of character-columns needed in each column:
X nchcol <- nsep+apply(nch,2,max)
X nchcol[nc] <- nchcol[nc] - nsep
X coltot <- sum(nchcol)
X rowtot <- ((nr-1)*rch+1)*(nr>0)
X return(chdim=c(rowtot,coltot),
X nchcol=nchcol,
X pcol=c(0,cumsum(nchcol[-nc]))
X )
X }
X##############################################################################
Xtableplot <- function(tabl,xlim,ylim,cex=par("cex"),adj=0,tabadj=c(.5,.5),
X borderleave=c(1,.5),rch=1.25,col=1){
X #function to plot in table-format a character matrix
X # - will have left-justified margins
X #inputs:
X # tabl - a matrix of character strings (varying lengths ok)
X # xlim - the x limits in the current plot (will be left-justified)
X # ylim - the y limits in the current plot (will be top-justified)
X # cex - the maximum cex to allow
X # adj - adjustment applied to text within each column (default=0=left)
X # tabadj - adjustment (x,y) applied to the entire table, if smaller than
X # available space
X # borderleave
X # - the minimum margins (x,y) between the table and the margins
X # of the plotting area. The table may be slightly shrunk in
X # order to leave the necessary border
X # rch - the number of character-heights to move to get to new row
X #
X #prepare to properly justify the text
X adjhold <- par("adj")
X on.exit(par(adj=adjhold))
X par(adj=adj)
X basecxy <- par("cxy")/par("cex")
X #get new cex:
X td <- chdim(tabl,rch=rch)
X checkborder <- T
X cexcol <- (xlim[2]-xlim[1])/(td$chdim[2]*basecxy[1])
X cexrow <- (ylim[2]-ylim[1])/(td$chdim[1]*basecxy[2])
X cexnew <- min(cex,cexcol,cexrow)
X cxynew <- basecxy*cexnew
X #now find what borders must be added
X xextra <- diff(xlim)-td$chdim[2]*cxynew[1]
X yextra <- diff(ylim)-td$chdim[1]*cxynew[2]
X xlimuse <- xlim+xextra*c(tabadj[1],tabadj[1]-1)
X ylimuse <- ylim+yextra*c(1-tabadj[2],-tabadj[2])
X usr <- par("usr")
X xchmarglack <- borderleave[1] - c(1,-1)*(xlimuse-usr[1:2])/cxynew[1]
X ychmarglack <- borderleave[2] - c(1,-1)*(ylimuse-usr[3:4])/cxynew[2]
X xchadd <- xchmarglack*(xchmarglack>0)
X ychadd <- ychmarglack*(ychmarglack>0)
X #recompute the character size
Xif(sum(xchadd)>0){
X xlim <- xlim+xchadd*cxynew[1]*c(1,-1)
X cexcol <- cexcol*td$chdim[2]/(td$chdim[2]+sum(xchadd))
X }
X if(sum(ychadd)>0){
X ylim <- ylim+ychadd*cxynew[2]*c(1,-1)
X cexrow <- cexrow*td$chdim[1]/(td$chdim[1]+sum(ychadd))
X }
X cexnew <- min(cex,cexcol,cexrow)
X cxynew <- basecxy*cexnew
X xextra <- diff(xlim)-td$chdim[2]*cxynew[1]
X yextra <- diff(ylim)-td$chdim[1]*cxynew[2]
X #get the coordinate matrices:
X #We don't subtract from 1 because
X #it appears that when characters are plotted with text(),
X #their centers are approximately -.1667*par("cxy")[2] below
X #the y coordinate at which they are plotted, and thus their
X #(including a margin) tops are .3333*par("cxy")[2] above
X #the y coordinate.
X txrow <- ylim[2] + ((1-.3333)-row(array(0,dim(tabl))))*rch*cxynew[2] -
X tabadj[2]*yextra
X colpos <- td$pcol+adj*td$nchcol
X txcol <- xlim[1] +matrix(colpos,nrow=nrow(tabl),ncol=ncol(tabl),
X byrow=T)*cxynew[1] + tabadj[1]*xextra
X text(c(txcol),c(txrow),labels=c(tabl),cex=cexnew,col=col)
X invisible()
X}
X##############################################################################
Xrankddo <- function(ddolist,CDname,rtype="poolrank"){
X #function to rank discrepancy distributions
X # according to specified ranking procedure
X #inputs:
X # ddolist - the list of discrepancy distribution objects
X # CDname - the name of the discrepancy functions
X # which will be ranked
X # rtype - the method of ranking:
X # poolrank: a Mann-Whitney variant
X # median: a "point scheme" for medians
X #outputs:
X # points - length(ddolist) by length(CDname) matrix of (deprecatory)
X # points attained by each ddo at each discrepancy fnctn
X # ranks - same dimensions as points; The ranks within each CD fnctn
X # of each element of points
X # opoints - vector size length(ddolist). Mean across CD fnctn type
X # of points. O stands for overall.
X # oranks - same size as opoints. Ranks of elements of opoints.
X #
X ndist <- length(ddolist)
X ndf <- length(CDname)
X allpoints <- allranks <- numeric()
X for ( CDnameuse in CDname){
X #
X #We will do one element of CDname at a time, then combine them.
X #
X lst <- list()
X points <- numeric(ndist)
X ninf <- ntot <- numeric(0)
X #first unpack the CDnameuse distributions from the ddos
X for(ddo in ddolist) {
X lst <- c(lst,list(ddo$CDall[,CDnameuse]))
X ni <- ddo$error$nerr
X if (is.null(ni)) ni <- 0
X ninf <- c(ninf,ni)
X ntot <- c(ntot,ni+dim(ddo$CDall)[1])
X }
X #
X #now branch by ranking type
X # and get points vector (one value each ddo)
X #
X pointsout <- switch(rtype,
X poolrank={
X ###
X ###rank in pooled Mann-Whitney system
X ###
X vals <- indx <- numeric()
X for(i in 1:ndist){
X indx <- c(indx,rep(i,length(lst[[i]])))
X vals <- c(vals,lst[[i]])
X }
X rvals <- rank(vals)
X infrnk <- length(vals)+1+.5*sum(ninf)
X for(i in 1:ndist) points[i] <- (sum(rvals[indx==i])+
X infrnk*ninf[i])/ntot[i]
X points
X },
X stop("poolrank is only ranking system implemented")
X )
X #######
X #######We have points for requested analysis type. Now get ranks.
X #######
X dfrank <- rank(pointsout)
X allpoints <- cbind(allpoints,pointsout)
X allranks <- cbind(allranks,dfrank)
X }
X allpoints <- matrix(allpoints,nrow=ndist,ncol=ndf)
X allranks <- matrix(allranks,nrow=ndist,ncol=ndf)
X dimnames(allpoints) <- dimnames(allranks) <- list(NULL,CDname)
X list(points=allpoints,ranks=allranks,
X opoints=apply(allpoints,1,mean),oranks=apply(allranks,1,mean))
X }
X##############################################################################
Xgetqua <- function(qua,dis,ndis=length(dis)){
X #function to get quantile estimates
X # of discrepancy distribution
X #inputs:
X # qua - a vector of desired quantile values
X # dis - unordered vector of discrepancy values
X # ndis - the number of discrepancy values calculated. Usually
X # is equal to length of dis (unless fittings out of
X # range)
X #output:
X # estimated qua quantiles, same order as original qua vector.
X #
X##################
X#note terrible bug in Splus for Windows: When two values in the
X#quantile argument are the same, the computer falls into an
X#infinite loop and must be rebooted.
X qv <- pmax(0,pmin(1,qua*ndis/length(dis)))
X qvu <- sort(unique(qv))
X res <- quantile(dis,qvu)
X res[match(qv,qvu)]
X }
X
X##############################################################################
Xgetquaci <- function(qua,dis,ndis=length(dis),CL=.95){
X #function for calculating confidence intervals for quantiles
X # inputs:
X # qua -vector of quantiles needing confidence intervals
X # CL -the confidence level for the confidence intervals
X # (default = .95)
X # dis - unordered vector of discrepancy values
X # ndis - the total number of discrepancy values calculated
X # (usually equal to length(dis), but sometimes some
X # fittings were out of range)
X #outputs:
X # qmin - vector of minimum bounds
X # qmax - vector of maximum bounds
X # (these bounds together comprise a 2*alpha level CI
X # or each alone constitutes an alpha-level CB)
X # qokmin - Logical vector. True if requested quantile for qmin
X # within [0,1] and thus was not truncated.
X # qokmax - Same as qokmin but for qmax.
X zlim <- qnorm(.5+.5*CL)
X #get the quantile difference away from estimated quantile
X qdiff <- sqrt(qua*(1-qua)/length(dis))*zlim
X qmin <- getqua(qua-qdiff,dis,ndis)
X qmax <- getqua(qua+qdiff,dis,ndis)
X qokmin <- 0<=(qua-qdiff)
X qokmax <- 1>=(qua+qdiff)
X return(qmin=qmin,qmax=qmax,qokmin=qokmin,qokmax=qokmax)
X }
X##############################################################################
Xgpar.plotdensity <- function(ddo,CDname=.DDO$CDname.store,excl=10,m=20,
X coscale=.DDO$coscale,
X coscalemenu=.DDO$coscalemenu,
X maxden=T){
X #Function to calculate the preferred plotting limits for
X # a discrepancy distribution object.
X #gpar.plotdensity() returns a vector containing, for each
X #coscaling group:
X # always:
X # one element for the X (discrepancy) plotting limit;
X # one element for the Y (density) plotting limit.
X #This function is used primarily by plotdensity; it is also used
X # by compare
X #inputs:
X # tdo - the discrepancy distribution object for which to calculate
X # the plotting limits
X # CDname - the discrepancy functions to be used
X # excl - the (plotdensity) excl value for which we are preparing
X # m - the (plotdensity) n value for which we are preparing
X # (how many points for density estimation)
X # maxden - T says we should calculate the maximum of the density
X # F says only calculate the maximum of the discrepancy axis
X # coscale - scaling standardization groups (default: .DDO$coscale).
X #outputs:
X # CDlim - a vector same length as CDname giving
X # plotting limits for respective variables.
X # denlim - (only produced if maxden=T) a vector same
X # length as CDname giving plotting limits for
X # vertical axes (densities) of respective
X # variables.
X #
X ###########>>>>
X #This function is best understood in the context of the
X #overall scaling system used in the programs. The algorithm
X #by which scales are chosen for compareddo() plots is as follows:
X #
X # 1. Obtain specifications (from appropriate source)
X #of which plot types should be scaled together, that is, what
X #are the coscaling groups. Furthermore, obtain the relative
X #scales which are to be forced (usually same size for all
X #plots, but sometimes different, as in example below). Example:
X # coscale=list(c(CD1=1,CD2=1,CDinf=2/3),
X # c(CA=1)
X # )
X #specifies that there should be two coscaling groups.
X #CD1, CD2, and CDinf are all to be scaled together, and CDinf
X #is to be forced to be plotted with units only 2/3 as large
X #as the other two. CA will also be plotted, but its scale
X #is allowed to "float" free of the others and expand or
X #contract as needed so as to fit all the data comfortably
X #into the available space in the graph.
X # 2. For each ddo determine the preferred
X #plotting coefficients for each individual coplotting group.
X #This involves for every variable (every "discrepancy"):
X # a. Eliminating any extremely large outliers
X # b. Finding the maximum of remaining values, but
X # going slightly larger for comfortable visibility.
X # c. Finding the minimum plotting coefficient which will
X # accomodate all variables within the group.
X #The plotting coefficient will be a number which when DIVIDED
X #by the coscaling vector (example above) for the coscaling group
X #will produce the plotting limits in the units of the individual
X #variables (however, density must be MULTIPLIED by coscaling
X #factor to preserve areas across plots).
X # 3. For each coplotting group, go through the
X #various ddo's and determine the coscaling coefficient
X #for each group which will accomodate all the ddo's
X #(this is the maximum of the coefficients needed by the
X #individual objects).
X # 4. Plot all objects and all variables according to
X #the above determined limits.
X #
X #Implementation:
X # 1. Argument coscale fed to compareddo() and plotdensity()
X #specifies the coscaling groups. A default exists in
X #.DDO$coscale.
X # 2. Function gpar.plotdensity() does this, returning plot limit
X #coefficients (coordinated within scaling groups)
X #for each variable within one individual ddo.
X # 3. Function compareddo() repeatedly calls gpar.plotdensity()
X #and keeps a running record of the maximum needed
X #coefficients for each variable.
X # 4. All actual plotting is done by function plotdensity(),
X #which is called by compareddo() when compareddo() plots densities.
X #
X couse <- coscalemenu[[coscale]]
X
X maxd <- rep(0,length(CDname))
X names(maxd) <- CDname
X CDlimuseall <- maxd
X for (dfnuse in CDname){
X #First find out which other CDs are in same
X #Coscaling group, and find the coscaling factor
X #For this particular CD. If this CD does not
X #belong to any coscaling group in the list in
X #coscale, it will be turned into its own coscaling
X #group with implied coscaling factor 1.
X csg <- 1
X names(csg) <- dfnuse
X #(the above is the default)
X for(cg in couse){
X if(!is.na(match(dfnuse,names(cg))))
X csg <- cg
X }
X #now turn csg into indices coordinated with CDname
X csgflag <- !is.na(match(CDname,names(csg)))
X csgfactor <- csg[CDname[csgflag]]
X
X #return to doing the actual computations.
X CD <- sort(ddo$CDall[,dfnuse])
X medianCD <- median(CD)
X maxuseCD <- excl*medianCD
X CD <- CD[CD<=maxuseCD]
X CDlimuse <- .5*medianCD+max(CD,3*medianCD)
X #Adjust all limits in the same coplotting group
X CDlimuseall[csgflag] <- pmax(CDlimuseall[csgflag],
X CDlimuse/csgfactor)
X if (maxden){
X if(length(CD)>1)
X md <- max(density(CD,n=m)$y)
X else
X md <- 0
X maxd[csgflag] <- pmax(maxd[csgflag],md*csgfactor)
X }
X #finish looping on dfnuse:
X }
X if (maxden)
X return(CDlim=CDlimuseall,denlim=maxd)
X else
X return(CDlim=CDlimuseall)
X }
X##############################################################################
Xat <- function(ddo){
X #function to create "(n=..)" format
X #input is a discrepancy distribution object
X nerr <- ddo$error$nerr
X ec <- ddo$evalcall
X if (is.null(nerr)) nerr <- 0
X aerr <- if(nerr==0) "" else paste(
X c("-","F+")[1+ec$nforce],nerr,sep="")
X q <- if (ec$n==nrow(ddo$beta)+nerr*(!ec$nforce)) "" else "?"
X paste("(n=",ec$n,aerr,q,")",sep="")
X }
X##############################################################################
Xprint.ddoorientation <- function(obj,width=pmin(80,options("width")[[1]]-1)){
X orie <- list(
X "[INTRODUCTION TO CALIBRATION DESIGN SOFTWARE - type \"ddo \"]",
X paste(" Welcome to the calibration experiment",
X "design package of functions",
X "for the Splus language. These functions, listed via help(index.ddo),",
X "are for creating and manipulating",
X "what we call \"discrepancy distribution objects\" or \"ddo\"s for short",
X "[see help(object.ddo)]. Each",
X "of these objects contains the results of repeated",
X "simulations of a calibration experiment. All simulations recorded in",
X "one ddo are the result of using one specified design with",
X "particular set of X values [see help(ddomake)]. To compare the",
X "performance of various different designs, the package allows you to assemble",
X "together several ddos into a type of object called a \"discrepancy distribution",
X "object frame\" or \"ddoframe\" for short [see help(ddoframe)].",
X collapse=" "),paste(
X " There are many specifications besides the design points which must be",
X "given in order to generate and analyze simulations.",
X "Rather than asking you respecify these parameters each time",
X "a simulation is performed",
X "the software insists that they be given default values.",
X "In order to set or observe these default parameters values you can",
X "use the function ddo() which operates",
X "similarly to the S function par() [see help(ddo)].",collapse=" "),paste(
X " To initialize background parameters the first time for a particular",
X "application, type ddobegin(\"application.name\"). The package comes equipped",
X "with a sample application which enables the user to evaluate designs for",
X "calibrating a machine for measuring pH level using fluorescene intensity ratio",
X "(fir). To begin this application or to get more information about it, type",
X "respectively ddobegin(\"phcal\") or help(phcal). For information about how to",
X "set up for your own new application, type help(application).",collapse=" "),paste(
X " For the theoretical background of the software,",
X "see \"Nonlinear Calibration Design\",",
X "an M.A. thesis by Doug Oman in Biostatistics at the University of California",
X "at Berkeley, 1995.", collapse=" "))
X for(str in orie){
X stmore <- T
X a<-1:nchar(str)
X stbeg <- 0
X st1 <- c(a[substring(str,a,a)==" "],nchar(str))
X while(stmore){
X if(length(st1)>1)
X stn <- (b<-st1[st1<=stbeg+width])[length(b)]
X else stn <- stbeg+width
X cat(paste(substring(str,stbeg+1,stn),"\n"))
X stmore <- stn'DDOEXAMP' <<'END_OF_FILE'
X#APPENDIX I from M.A. Thesis by D.Oman:
X#COMPUTER CODE FOR pH CALIBRATION ANALYSIS
X
X#Note: This file would take a very long time to execute as
X#a source() file. Before experiment with it, you may want
X#to change the values of nstart and njackup assigned below.
X
X
X################################################################
X################################################################
X# This file of Splus code contains the commands used
X# to generate FIGURES 3, 5, 6, 7 and the example in
X# Section VII, "EXAMPLE: CHOOSING A DESIGN FOR pH CALIBRATION"
X# in the Master's Thesis, "NONLINEAR CALIBRATION DESIGN,"
X# by Doug Oman, University of California at Berkeley, 1995.
X#
X# Many selections of the narrative text from the thesis are
X# interspersed as comments in the code file so that the
X# correspondence between the code and the narrative is easy
X# to reconstruct. The results of the simulations, however,
X# will not be identical because of differences in the random
X# seed at the time of the original execution.
X################################################################
X#FIG 3: INTRINSIC ABERRANCE AND LINEARLY APPROXIMATED POINTWISE DISCREPANCY
X#FIG 5: DISPLAY OF DESIGN AND DESIGN-LEVEL SUMMARIES
X#FIG 6: INFERENTIAL RESULTS IDENTIFYING TWO BEST CANDIDATES
X#FIG 7: RELATIVE EFFICIENCIES OF (5 DESIGNS) $times$ (3 OR 2 BATCHES)
X################################################################
X#S code for initializing setup
X#
X #Initialize various values specifying how we will do the simulation
Xnstart <- 200 #number for initial simulations
Xnjackup <- 1000 #number for more extensive simulations
X
XXdesignrange <- c(6.2,8.4) #used in function designphcal() (see below)
XXconceivable <- c(6.6,8.2)
X
X #make sure the code knows we want to construct and evaluate
X #according to functions named "phcal"
Xddobegin("phcal")
X #now create particular values specified in current application
Xroix <- seq(from=6.8,to=8.0,length=15)
Xroiw <- 1+(roix<=7.3)
Xddo(roimenu=c(.DDO$roimenu,list(fall94=list(x=roix,w=roiw))))
Xddo(roi="fall94")
X
X #now specify the plotting range and also specify that we
X #will force there to be the number of COMPLETED simulations
X #that we specify; failure of convergence is always counted
X #as a "failure", and because nforce=T will need to be redone
X #and will not count toward the n repetitions
Xddo(Xplotrange=c(5.5,8.3),nforce=T)
X
X #specify the range outside of which we truncate
Xddo(Xconceivable = Xconceivable)
X
X #Using specially created (for this pHcal application) function designphcal(),
X #create values such as rl, om, etc. to be used in specifying designs.
Xdesignphcal()
X#note om(6.71) ~ rl(6.8) and os(7.28) ~ rm(7.4)
Xol <- Xdesignrange[1]
Xor <- Xdesignrange[2]
X
X
X
X################################################## BEGINNING OF NARRATIVE
X
X #before we do anything, look at intrinsic aberrance
X ###################### FIG 3
X # "The pH calibration machine had a change of 430% betwen the
X # largest and smallest SDs for $PD sub linear (X)$ over
X # ROI=[6.8, 8.0] [see FIG 3, which was generated by a
X # single software command in Splus]."
Xintrinsic() #produces FIG 3
X
X #Specify designs of interest
X #Our strategy will be to try designs of 7 points
X #at 3 identical batches to do a preliminary
X #screening.
X #later, for the good designs we will experiment with only 2 batches
Xddo(batch=3)
X
X
X
X
X # "The first step in answering the above questions was to
X # generate some candidates for good designs to be
X # subjected to evaluation by the software. This step was perhaps the most
X # challenging. We generated a wide array of candidate designs
X # by drawing upon an amalgam of heuristics from both traditional
X # clinical practice and statistical theory
X # ....
X # This process resulted in twelve candidate designs of seven points
X # each which were to be repeated either two or three times, depending
X # upon how many batches were wanted; we hoped by using this diverse
X # collection of starting points that we would not miss finding the
X # particular local minimum of the design-level discrepancy function
X # $DD( delta )$ which contained the global minimum. These calibration
X # designs were entered into the Splus calibration design software
X # environment."
X
Xdall <- list(
X"NARROW,WIDE~{3(RL),1(RM),3(RR)}" = c(rep(rl,3),rm,rep(rr,3)),
X"NAR,WIDE,TIEa~{1(OL),3(RL),3(RR)}" = c(ol,rep(rl,3),rep(rr,3)),
X"NAR,WIDE,TIEb~{3(RL),3(RR),1(OR)}" = c(rep(rl,3),rep(rr,3),or),
X"NARROW,SPREAD~{7(RL:RR)}" = seq(rl,rr,length=7),
X"NARROW,CLUSTER~{3(RL),2(RM),2(RR)}" = c(rep(rl,3),rep(rm,2),rep(rr,2)),
X"WIDE,SPREADa~{7(OL:OR)}" = seq(ol,or,length=7),
X"WIDE,SPREADb~{7:(OL:RR)}" = seq(ol,rr,length=7),
X"WIDE,SPREADc~{7(RL:OR)}" = seq(rl,or,length=7),
X"WIDE,CLUSTERa~{2(OL),3(RM),2(OR)}" = c(rep(ol,2),rep(rm,3),rep(or,2)),
X"WIDE,CLUSTERb~{2(OL),3(RM),2(RR)}" = c(rep(ol,2),rep(rm,3),rep(rr,2)),
X"WIDE,CLUSTERa~{2(RL),3(RM),2(OR)}" = c(rep(rl,2),rep(rm,3),rep(or,2)),
X"WIDE,CLUSTERb~{2(RL),3(RM),2(RR)}" = c(rep(rl,2),rep(rm,3),rep(rr,2)),
X )
Xdtrad <- 4 #number of the traditional design
X
X # "We used the calibration design software to simulate at n=200
X # repeats each of the twelve candidate designs."
Xddo(n=nstart)
Xmake("f",dall)
X
X ###################### FIG 6
X # "...two of the designs seemed superior, with all other designs being
X # significantly worse (p<.05) than the two superior designs in one or more
X # of the design-level discrepancies $DD( delta )$. See [FIG 6] for the
X # command output displaying the results of the comparison...
Xcompareddo(f) # <- FIG 6
X #For figure 6 alone: compareddo(f,anal="m")
X
Xbestf <- c(6,2,11,8)
X #To look at the best among all:
X#compareddo(f[bestf])
X #To examine all of our default values:
X#ddo()
X
X # "Of the remaining ten designs, two seemed best, so we selected these
X # four designs and resimulated at n=1000 repeats."
Xddo(n = njackup)
Xmake("fju",dall[bestf])
X #Using the software's data management commands,
X #we could also specify this with make("fju",dall[bestf],n=njackup)
X #
X #To find out if the same relationships hold between the design efficiencies
X #under this larger number of repetitions, use the easy manipulation
X #of ddoframe objects to compare them side-by-side:
X
X # "Essentially the same inferential results held, although one of the two best designs
X # ('Wide/Spread/a') was found to be preferable to the other according to one of
X # the design-level measures ($DD sub 1$) at p<.05)
Xcompareddo(c(fju,f[bestf])[c(1,5,2,6,3,7,4,8)])
X
X # "We wondered if slight variations of the positions of the points
X # would make a difference, so we generated some variations of the best
X # designs (we 'tinkered') by moving some of the design points in ways
X # we thought could improve the design."
X
Xdref <- c(dall[c(6,2)],list(
X"tink:lrNWTa~{1(OL),3(RL:7.1),3(OR)}" = c(ol,seq(rl,7.1,length=3),rep(or,3)),
X"tink:l-NWTa~{1(OL),3(RL:7.1),3(RR)}" = c(ol,seq(rl,7.1,length=3),rep(rr,3)),
X"tink:-rNWTa~{1(OL),3(RL),3(OR)}" = c(ol,rep(rl,3),rep(or,3)),
X ))
X
X
X
Xmake("fref",dref)
X #The plot command, when given a ddoframe object, operates
X #equivalently to compareddo()
X # "We simulated these five designs (re-re-simulating the two originals,
X # and simulating the three new variants) at 1000 repeats
X # and found that on this simulation the best design was one of the
X # newly generated designs, entitled 'tink:lrNWTa~{1(OL),3(RL:7.1),3(OR)}'
X # (see FIG 5 for its solo display).
Xplot(fref)
X ###################### FIG 5
X #FIG 5: DISPLAY OF DESIGN AND DESIGN-LEVEL SUMMARIES
X # [out of order]
X # "One command displays the design and numerically summarizes the
X # design-level summaries of the simulated distributions [FIG 5]."
Xshowddo(fref3) ## <- FIG 5
X
X #could look at the scatterplots and strata of the discrepancies
X#focus(fref1)
X #focus() calls in turn showddo(), plotdensity(), plotstrata(),
X # and pairddo()
X
X#It seems that the most interesting five designs are
X# the four best from the original list and new design
X# we generated by tinkering. For these designs, let us try
X# some further refining and see what happens.
X
X
X#What happens to efficiency if we reduce from 3 batches to 2 batches?
X #create a 1000-simulation of the traditional design,
X #using command to create single design simulation
Xftrad <- ddomake(dall[[dtrad]],n=njackup,main=names(dall)[[dtrad]])
X #pull out the most interesting old fittings
Xfthree <- c(ftrad,fref[1:3])
X #pullout the corresponding designs
Xdtwo <- c(dall[dtrad],dref[1:3])
X #an alternative way to do this is by pulling out titles
X #and designs using extraction commands
Xdtwoalso <- fthree[,"xdesign",simplify=F,drop=T]
Xnames(dtwoalso) <- fthree[,"main"]
Xcat("\ndtwo = dtwoalso:\n")
Xprint(dtwo)
Xcat("\nall.equal(dtwo,dtwoalso):",all.equal(dtwo,dtwoalso),"\n\n")
X #refit and look for the best of the two-batch designs
Xddo(batch=2)
Xmake("ftwo",dtwo)
Xplot(ftwo)
X #now compare the efficiencies of the three-batch and
X #two-batch designs:
Xplot(c(fthree,ftwo),analysis="Ee") ### <- FIG 7
X #For figure alone: plot(c(fthree,ftwo),analysis="e")
X ###################### FIG 7
X #FIG 7: RELATIVE EFFICIENCIES OF (4 DESIGNS) $times$ (3 OR 2 BATCHES)
X # "To determine whether the efficiency gained from the irregular design
X # was worth the energy, we examined the relative efficiencies of the
X # designs (see [FIG 7], designs 1-4, i.e., the four left-hand points in
X # each plot).
X # ...
X # Designs 5 through 8 of [FIG 7] show the same designs as discussed
X # above resimulated at two batches for n=1000 repeats....efficiency was
X # reduced by slighly less than 20% in each design discrepancy measure..."
X
Xefdif <- (ftwo[,"DDall"][1,,] / fthree[,"DDall"][1,,] -1 )*100
Xdimnames(efdif)[[2]] <- c("Traditional","Best of 12 (chosen)",
X "Second best of 12","Very best (newly generated)")
Xprint(rbind(t(round(efdif,1)),"mean (all)"=apply(efdif,1,mean),
X "mean (all but Traditional)"=apply(efdif[,-1],1,mean)))
END_OF_FILE
if test 10042 -ne `wc -c <'DDOEXAMP'`; then
echo shar: \"'DDOEXAMP'\" unpacked with wrong size!
fi
# end of 'DDOEXAMP'
fi
if test ! -d '.Data' ; then
echo shar: Creating directory \"'.Data'\"
mkdir '.Data'
fi
if test ! -d '.Data/.Help' ; then
echo shar: Creating directory \"'.Data/.Help'\"
mkdir '.Data/.Help'
fi
if test -f '.Data/.Help/application' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'.Data/.Help/application'\"
else
echo shar: Extracting \"'.Data/.Help/application'\" \(7300 characters\)
sed "s/^X//" >'.Data/.Help/application' <<'END_OF_FILE'
X.BG D
X.FN application
X.TL
XHow to implement an application (such as "phcal") using the software.
X.PP
XImplementation of an application requires the creation of several functions,
Xall with suffixes ".apname" where "apname" is the name of the application.
XWe shall illustrate the method of implementation through the phcal (pH Calibration)
Xexample, which is contained in the software. The following routines must be
Xcreated:
X
X.ce
XREQUIRED ROUTINES
X.IP ddo.apname
XFunction to initialize or reset default parameter values.
XThis function should have one argument for each default parameter
Xto be specified. For example, ddo.phcal has 18 arguments and its first line
Xbegins "ddo.phcal <- function(mu = phrab.bestfit1$coef$fixed, n = 200,
Xnforce = F,...". This assigns a default value to mu (which is drawn from
Xthe object phrab.bestfit1, a fitting of the original pH calibration data); sets the
Xdefault number of repeat simulations for each object to n=200; etc. The body
Xof the function should then assemble these arguments into a single list item, which is
Xoutput (see ddo.phcal). Each ddo.apname must define the following parameters:
X
X.ce
XRequired arguments to ddo.apname:
X.LP
Xmu, n, nforce, fam0, family, Xconceivable, measurement.system,
XCDname.store, droperr, CDname.all,
XXdesignrange, Xplotrange, coscale, coscalemenu, roi, roimenu.
XMost of these parameters (mu, n, nforce, fam0, family,
XXconceivable, measurement.system, CDname.store, droperr) are default parameters used
Xby ddomake(). The others have the following meanings: CDname.all is a vector
Xof all calibration-discrepancy names calculable by discrepancyC.apname, and should have
Xan attribute "PDname.all" which is a character vector with elements either "PD" or "PA" which
Xspecifies (for plotstrata()) whether the CD is a bias-like discrepancy (PD) or an aberrance (PA).
XXdesignrange is a
Xdefault range for possible design points. Outside of Xdesignrange,
Xthe ddomake() routine refuses to allow design points; Xplotrange gives the
Xdefault plotting range used by plotting routines (intrinsic, showx.default, gpar.showx.default);
Xcoscale gives specifications for how the axes of the various calibration-level densities
Xshould be scaled together (used by compareddo, plotdensity, gpar.plotdensity), and coscalemenu
Xis a menu of options for coscale. The argument roi specifies (but is not itself) the default range of interest
Xused by ddomake(). The roi parameter is run through the roimake() function by ddomake: if of
Xmode character, roi specifies a menu entry in roimenu; if of mode list, roi must be equivalent to an output
Xof roimake; if of mode numeric, roi is interpreted as a vector of X points to be used (equally weighted)
Xas the range of interest.
X
X.ce
XOptional arguments to ddo.apname (specific to phcal):
X.LP
XSeveral parameters in the phcal example are optional and need not be used in every application, but
Xare illustrative of types of extra default parameters of interest.
XThe argument sigma0 is the default residual variance used by simulation.phcal (see discussion below).
XThe argument N is the number of repeat measures done at each design point (the effect is simply to
Xreduce the error variance used in the simulation by a factor of N). The argument batch is the default
Xnumber of repeat batches to make of the design vector xdesign passed from ddomake().
X
X.IP yofx.apname
XFunction to calculate Y from X.
XMust take three arguments (beta, X, family). Arguments can have any names but must
Xbe in that order. The meanings of the arguments: (1) beta is a length three vector of parameters.
X(2) X is a scalar, matrix, or array of X values. (3)
Xfamily is a character vector specifying function family. The output is of the same
Xdimension as X, and is the Y values corresponding to the X value according
Xto beta and family.
X
X.IP xofy.apname
XFunction to backcalculate X from Y.
XMust take three arguments (beta, Y, family). Arguments can have any names but must
Xbe in that order. The meanings of the arguments: (1) beta is a length three vector of parameters.
X(2) Y is a scalar, matrix, or array of Y values. (3)
Xfamily is a character vector specifying function family. The output is
Xof the same dimension as Y, and is the backcalculated values of X
Xaccording to beta and family.
XIf Y is aberrant, an NA should be output.
X
X.IP simulate.apname
XFunction to simulate (and fit) a (single) calibration experiment
XRequires at least four argument in general (xdesign, mu, fam0, family),
Xbut individual applications such as phcal may include other arguments.
XThe four arguments should the first four arguments and in the above order.
XAny other arguments may be specified and should either have defaults or
Xthe user should be prepared to specify them each time ddomake() or make()
Xis called. A very important type of other variable is the specification of the
Xvariance components to be used in the simulation. These are not named and
Xpassed through ddomake() by name because it is possible that some
Xapplications may want to use models involving random effects. The output
Xshould be a list with four elements: beta, a vector of fitted parameters;
Xerror, TRUE if there was a failure of convergence; errinfo, NULL if convergence
Xachieved but if nonconvergence, a list of information about the nonconvergence;
Xarg, a list of the evaluated parameters used by simulate.apname which are
Xapplication-specific (must be passed to it through the ... argument in ddomake()).
X.PP
XExample: simulate.phcal contains the other arguments sigma, N sigmabatch,
Xand start. Each of these is given a default value in the first line of simulate.phcal,
Xbut their default values come from a variety of sources.
XThe sigma default comes from a separate function sigmaiidphcal() which
Xin turn draws upon several sources, including ddo("sigma0");
XN comes from ddo("N"); sigmabatch
Xis set to 0 as a default because the possibility of batch variation was never
Ximplemented in phcal; and start (starting values for fitting beta) is set to mu,
Xthe true parameter values already passed to simulate.phcal as another argument.
X.ce
XOPTIONAL ROUTINES:
X.PP
XAll of the following routines are optional and may be omitted. For applications similar to phcal,
Xusers may want to write only the first two routines which update the printing methods to
Xinclude application-specific parameters. These two methods have been written for phcal; For
Xthe other four, phcal relies upon default methods.
X.IP printDDO.apname
XPrint application-specific default parameters (when ddo()) is invoked).
X.IP printdetail.apname
XPrint application-specific parameters used for constructing a particular
Xclass ddo object (when print(object.ddo) is invoked).
X.IP discrepancyC.apname
XCalculate calibration-level discrepancies for constructing a class ddo object.
X.IP discrepancyD.apname
XCalculate design-level discrepancies for constructing a class ddo object
X.IP gpar.showx.apname
XGet default graphing parameters such as margins, etc.
X(used by routine showx.apname) for plotting a design.
X.IP showx.apname
XPlot a design and various auxiliary information.
X.IP showxddo.apname
XPlot the design associated with a class ddo object.
X.KW sysdata
X.WR
X
X
END_OF_FILE
if test 7300 -ne `wc -c <'.Data/.Help/application'`; then
echo shar: \"'.Data/.Help/application'\" unpacked with wrong size!
fi
# end of '.Data/.Help/application'
fi
if test -f '.Data/.Help/compareddo' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'.Data/.Help/compareddo'\"
else
echo shar: Extracting \"'.Data/.Help/compareddo'\" \(3078 characters\)
sed "s/^X//" >'.Data/.Help/compareddo' <<'END_OF_FILE'
X.BG
X.FN compareddo
X.TL
XFunction to compare several class ddo objects.
X.DN
XThe compareddo() function runs several analyses which are specified
Xin the "analysis" argument (or its default) upon the class ddo objects
Xwhich are specified in its first argument. These analyses are each
Xgiven one letter names:
X.IP (D)
Xshow designs in original order;
X.IP (d)
Xshow designs ranked in order roughly according to efficiency;
X.IP (e)
Xplot efficiencies of DDuse, the preferred design level discrepancy
Xcorresponding to each Calibration-Level discrepancy (e.g., Interquartile
Xfor Calibration-level aberrance);
X.IP (E)
Xplot efficiencies of all DDs;
X.IP (f)
Xplot densities;
X.IP (m)
Xcalculate Mann-Whitney statistics.
X.CS
Xcompareddo(distlist, CDname=.DDO$CDname.all, CL=0.95,
X mwspread=c(2, 100), analysis="dem", coscale=.DDO$coscale,
X DDElist=c("DDmean", "DDmedian", "DDinterq", "DDq95"),
X maxsize=3, maxrow=6, tag=date())
X.RA
X.AG distlist
Xa list of distribution objects to be compared, or else a
Xclass ddoframe object
X.OA
X.AG CDname
Xthe names of the discrepancy functions to be used in the
Xcomparison.
X.AG CL
Xthe confidence level to be plotted in analyses "e" and "E"
X.AG mwspread
Xa specification of how many Mann-Whitney analyses to do in
Xanalysis "m", i.e., should we do them even for ddos of very
Xdifferent ranks (for details, see code).
X.AG analysis
XA single character string which should contain characters of
Xeach desired analysis, from among the six characters D, d, e, E, f, m.
X.AG coscale
Xscaling groups used by analysis f(see gpar.plotdensity()).
X(default: .DDO$coscale)
X.AG DDElist
Xthe types of DD measures to be plotted in analysis E.
X.AG maxsize
Xlargest # of frame rows and columns per page for D, d analyses
X.AG maxrow
Xlargest number of rows per page for f, F analyses
X.AG tag
Xwhat should be stamped in the lower right corner of
Xeach plot
X.RT
XNone.
X.SE
XEntire output is through plots.
X.DT
XThe initial ordering for analyses d, f, and m is obtained from
Xfunction rankddo(). The default method is to pool all
Xof each kind of discrepancy (e.g., CD1) across all class ddo
Xobjects and calculate average ranks. These average ranks are
Xused to give the class ddo objects a ranking from 1 through nobj,
Xwhere nobj is the number of class ddo objects. A different ranking
Xis done for each kind of discrepancy, and then these 1-to-nobj
Xranks are averaged across the discrepancy types, to give each class
Xddo object an average rank somewhere in the range from 1 to nobj.
XThese average ranks are used for the orderings produced; This
Xmethod has many theoretical flaws but is relatively simple
Xcomputationally and helps the user's eye sort out which designs
Xare of interest.
X.SH REFERENCES
XSee reference in print(ddo) for the theoretical construction
Xof the CA, CD, DA, and DD measures.
X.SA
X plot.ddo()
X ddoframe
X.EX
X compareddo(list(ddo1,ddo2))
X #compare class ddo objects ddo1, ddo2
X compareddo(fref)
X #compare objects in class ddoframe object fref
X compareddo(fref[1:3],anal="e")
X #compare efficiencies of first three
X #class ddo objects in ddoframe fref
X.KW ~keyword
X.WR
END_OF_FILE
if test 3078 -ne `wc -c <'.Data/.Help/compareddo'`; then
echo shar: \"'.Data/.Help/compareddo'\" unpacked with wrong size!
fi
# end of '.Data/.Help/compareddo'
fi
if test -f '.Data/.Help/ddo' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'.Data/.Help/ddo'\"
else
echo shar: Extracting \"'.Data/.Help/ddo'\" \(2687 characters\)
sed "s/^X//" >'.Data/.Help/ddo' <<'END_OF_FILE'
X.BG
X.FN ddo
X.TL
XDiscrepancy Distribution Objects
X.DN
Xddo has three meanings (the letters stand for "discrepancy distribution object"),
Xonly the first of which is fully described in the current help file:
X
X (1) ddo() is a function which views or resets default parameter values
X(described in the current help file, help(ddo))
X
X (2) ddo is the name of the shareware system for nonlinear calibration design;
X(basic information may be obtained by typing "ddo "; an
Xindex of functions is contained in help(index.ddo))
X
X (3) ddo is a class of objects which has several associated methods
X(described in help(object.ddo))
X
XThis help() file describes the ddo() function, which operates very much like
Xthe par() function in its ability to view or reset parameters.
X.CS
Xddo(..., initialize=NULL, menu=F)
X.RA
X.AG ...
XIf unnamed, these arguments must be character strings which are names of
Xdefault parameters.
XIf named, these arguments must be the names of (default) parameters
Xand their values are
Xthe new values of these (default) parameters.
X.OA
X.AG initialize
XUsed by the ddobegin() function when it initializes the entire setup.
XUnneeded by most users.
X.AG menu
XIf TRUE, requests that menu parameters be printed. These are lists
Xwhich allow selection of ranges of interest and coscaling for density plots.
X.RT
XIf no arguments are present, ddo() prints out all default parameter values (except
Xfor menu parameters).
XIf the arguments are unnamed and are character strings, the function returns
Xa list of the named parameter values.
XIf the arguments are named, they are requests to reset the corresponding
Xdefault parameters.
X.SE
XMay reset default parameter values.
X.DT
XAll default parameters are stored in object .DDO.
XIf called with an "initialize=" argument, this will
Xinitialize values for the application with the name given
Xin the argument (e.g., initialize="phcal").
XThere must be a ddo.[application name] for each application
Xwhich outputs a list of the parameters needed for that
Xapplication. Thus when ddobegin() calls ddo(initialize="phcal"),
Xddo() in turn calls ddo.phcal which returns the original default
Xparameter values, which are then stored in .DDO.
X.SH NOTE
XThis function is saved as a class "ddoorientation" object so that when
Xits name is typed, a print method for "ddoorientation" prints out the
Xbasic orientation to the shareware system.
X.SA
X ddobegin()
X ddomake()
X make()
X.EX
X ddo(menu=T)
X #prints default parameter values, including menus.
X ddo(n=100)
X #sets the default number of simulations to 100.
X ddo(roi="fall94")
X #sets range of interest to the "fall94" choice from
X #roimenu list (menus only available for certain
X #parameters)
X.KW ~keyword
X.WR
END_OF_FILE
if test 2687 -ne `wc -c <'.Data/.Help/ddo'`; then
echo shar: \"'.Data/.Help/ddo'\" unpacked with wrong size!
fi
# end of '.Data/.Help/ddo'
fi
if test -f '.Data/.Help/ddobegin' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'.Data/.Help/ddobegin'\"
else
echo shar: Extracting \"'.Data/.Help/ddobegin'\" \(909 characters\)
sed "s/^X//" >'.Data/.Help/ddobegin' <<'END_OF_FILE'
X.BG
X.FN ddobegin
X.TL
XFunction to initialize a particular application (e.g., phcal)
X.DN
XInitializes default parameters for a particular application.
XShould be called whenever installing ddo system in a particular
Xdirectory, or when shifting from one application (e.g., the phcal
Xexample that comes with the shareware)
Xto another (e.g., your own).
X.CS
Xddobegin(application="phcal")
X.RA
XNone
X.OA
X.AG application
XThe name of the application as a character string
X.RT
XNone
X.SE
XInitializes default parameter values in object .DDO
X.DT
XObtains default parameter values from function ddo.[application name]
Xand then stores in object .DDO. The function ddo.[application name]
Xmust be supplied by user (other than for phcal example).
X.SA
X application
X ddo
X.EX
X ddobegin() #initializes pH calibration example
X ddobegin("userap") #initializes "userap" application.
X #User must supply routine ddo.userap
X.KW ~keyword
X.WR
END_OF_FILE
if test 909 -ne `wc -c <'.Data/.Help/ddobegin'`; then
echo shar: \"'.Data/.Help/ddobegin'\" unpacked with wrong size!
fi
# end of '.Data/.Help/ddobegin'
fi
if test -f '.Data/.Help/ddoframe' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'.Data/.Help/ddoframe'\"
else
echo shar: Extracting \"'.Data/.Help/ddoframe'\" \(5010 characters\)
sed "s/^X//" >'.Data/.Help/ddoframe' <<'END_OF_FILE'
X.BG D
X.FN ddoframe
X.TL
XObject class ddoframe
X.PP
XThe ddoframe class of objects is a way of organizing class ddo objects so that
Xthey are easier to create, manipulate, and compare. There is an intentional
Xparallel to the data frame class in Splus. Separate class ddo objects
X(created from separate designs, in general) are like rows, and the various
Xelements of the class ddo objects which can be referenced using []
X(such as n, etc.) are like the columns.
XBut since class ddo objects are
Xlarge, class ddoframe objects are created as "virtual" objects which REFER
Xto class ddo objects which are stored in memory. Class ddoframe objects
Xare designed with the idea that usually they will be first generated by
Xthe make() function (see help(make)).
XWhen the name of a class ddoframe object is typed, the following four types of
Xinformation are printed:
X.IP (1)
XThe application (such as phcal) to which all the objects
Xbelong (since make() only operates within one application at a time). The number
Xof class ddo objects is also printed, and the number of elements which
Xcan be referenced by bracket subscripting [] in
Xthe ddo objects.
X.IP (2)
XThe names under which the class ddo objects are stored in memory. When a ddoframe
Xobject has been created by make(), these names will be of the form prefix+(sequential numeric
Xsuffix), as described in help(make).
X.IP (3)
XThe elements which can be referenced in the class ddo objects by []. There will usually be
Xabout 30 elements and they include all the evaluated arguments to the ddomake() command
Xwhich produced the class ddo objects, as well as the outputs from ddomake() such as CDall,
XDDall, beta, and rseedend. The output nerr is also referencable (it is zero when no errors
Xoccurred) but the other error information is not referencable by [].
X.IP (4)
XThe time at which the ddoframe object was created.
X.PP
XClass ddoframe objects (and their constituent class ddo objects)
Xcan then be manipulated by using several functions and methods. The most powerful
Xof these functions is the subscripting [] function, which is listed first. Examples
Xare here given of their use upon class ddoframe objects named f and a shorter (i.e.,
Xfewer class ddo member objects) class ddoframe object g, and class ddo objects ddo1, ddo2, ddo3:
X
X.IP [.ddoframe
X.PP
XThis function takes either 1 or 2 main arguments, along with the possibility
Xof secondary arguments such as virtual, drop, and simplify. The meaning
Xof the syntax varies somewhat depending upon which arguments are included,
Xas seen below:
X f[1:3] #length three ddoframe object, the first 3 class ddo objects in f.
X f[1:3,virtual=F] #a LIST of three class ddo objects.
X f[1] #a class ddo object, the first ddo object in f
X f[1,drop=F] #a class ddoframe object (see "[" for parallel use of drop)
X f[1:3,"n"] #a length three vector of the number of simulations
X f[1:3,"n",simplify=F] #a length 3 list of n from the first three
Xclass ddo objects in f (see sapply for parallel use of simplify).
X f[1:3,"xdesign"] #if all designs of same length, a length*3 dimension
Xmatrix; otherwise a length 3 list of vectors.
X f[,"main"] #a vector of length(f) of the titles of each
Xclass ddo object in f.
X f[,elements(f)] #returns a list of length equal to length(elements(f)).
XEach element of the list is named and corresponds to f[,elements(f)[i]] as
Xi ranges in 1:length(elements(f)).
X
X.IP names.ddoframe
X.PP
Xnames(f) #tells the names of the class ddo objects constituting f.
X.IP elements
X.PP
Xelements(f) #tells the elements of the class ddo objects in f (similar
Xto names(object.ddo), described below).
X.IP dimnames.ddoframe
X.PP
Xdimnames(f) #equivalent to list(names(f),elements(f))
X.IP length.ddoframe
X.PP
Xlength(f) #equivalent to length(names(f))
X.IP dim.ddoframe
X.PP
Xdim(f) #equivalent to c(length(f),length(elements(f)))
X.IP names<-.ddoframe
X.PP
Xnames(f)[1:length(g)] <- names(g) #replaces the first length(fref) virtual
Xelements in f by the elements of g
X.IP c.ddo
X.PP
Xc(ddo1, ddo2, ddo3) #combines ddo objects into a ddoframe
X.IP c.ddoframe
X.PP
Xc(f,g) #creates a new ddoframe of length equal to length(f)+length(g)
X.IP [<-.ddoframe
X.PP
Xf[2] <- ddo1 #makes second element of ddoframe f equal to class ddo object ddo1.
X.IP print.ddoframe
X.PP
X f #prints a summary of ddoframe object f as described above
X for(i in length(f)) print(f[i]) #prints out each class ddo member of f
X.IP names.ddo
X.PP
Xnames(ddo1) #the elements accessible by [] in class ddo object ddo1 (see help(object.ddo))
X.IP [.ddo
X.PP
X ddo1["CDall"] #the matrix of calibration-level discrepancies values from class ddo object ddo1
X.PP
X ddo1[c("n","xdesign")] #a length two named list of the n and xdesign elements of ddo1
X
X.LP
XClass ddoframe objects are accepted as input to compareddo(), (and plot() which
Xcalls compareddo()).
X.SA
X object.ddo()
X make()
X compareddo()
X.KW sysdata
X.WR
X
X
END_OF_FILE
if test 5010 -ne `wc -c <'.Data/.Help/ddoframe'`; then
echo shar: \"'.Data/.Help/ddoframe'\" unpacked with wrong size!
fi
# end of '.Data/.Help/ddoframe'
fi
if test -f '.Data/.Help/ddomake' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'.Data/.Help/ddomake'\"
else
echo shar: Extracting \"'.Data/.Help/ddomake'\" \(5004 characters\)
sed "s/^X//" >'.Data/.Help/ddomake' <<'END_OF_FILE'
X.BG
X.FN ddomake
X.TL
XFunction to create repeated simulations of a single design.
X.DN
X~brief description
X.CS
Xddomake(xdesign, mu=.DDO$mu,
X roixy=roimake(name = .DDO$roi, mu = mu, fam0 = fam0),
X measurement.system=.DDO$measurement.system,
X Xconceivable=.DDO$Xconceivable,
X CDname.store=.DDO$CDname.store,
X droperr=.DDO$droperr, n=.DDO$n, nforce=.DDO$nforce,
X fam0=.DDO$fam0, family=.DDO$family,
X extend.info=list(DDget = T),
X randomseed=.Random.seed,
X main=paste("UNTITLED (", date(), ")", sep = ""),
X tag=T, warn=T, ...)
X.RA
X.AG xdesign
X~Describe xdesign here
X.OA
X.AG mu
XThe "imagined true" parameter values used as the basis of the simulation
X.AG roixy
XA list of range of interest specifications corresponding to an
Xoutput of roixy(). See help(roixy).
X.AG measurement.system
XThe type of complete measurement system used to backcalculate values for Xmeasured
Xwhen the simulated Ymeasured is aberrant. Currently supported options are
X"truncation" and "repeat.measurements".
X.AG Xconceivable
XThe conceivable range of X values. Necessary for dealing with aberrant measurements.
X.AG CDname.store
XThe names of the calibration level discrepancies to be calculated and stored in the returned
Xclass ddo object.
X.AG droperr
XIf FALSE, design-level quantile-based discrepancies (such as
Xmedian CD) are calculated assuming that there is
Xinfinite discrepancy for fittings that failed to converge.
X.AG n
XThe number of simulations to perform.
X.AG nforce
XIf TRUE, routine keeps simulating until there are n successful simulations;
Xif FALSE, only n total attempts are made to simulate.
X.AG fam0
XThe family of functions to be used as the basis for the simulation, e.g.,
Xthe "imagined true" family of calibration curves.
X.AG family
XThe family of functions to be fitted (by least squares) to obtain the calibration curve.
XFor example, the phcal example allows family="quadratic" to permit investigation
Xof how well such a fitting performs.
X.AG extend.info
XReceives information from the extend() function when it calls ddomake. Not of
Xinterest to most users.
X.AG randomseed
XThe random seed at which to begin the simulation
X(default: current .Random.seed).
X.AG main
XThe title of the design (used later by printing and plotting routines).
X.AG tag
XIf TRUE, type date and "FINISHED" when finished executing (default=T).
X.AG warn
Xif TRUE, will emit run-time warnings as they happen (default=T).
X.AG ...
XAdditional arguments which can be passed on to the (application-specific)
Xsimulation routines. They may include
X starting values (for fittings)
X imagined "true" residual SD
X imagined "true" other variance components, if they exist
X number of batches
X.RT
XA class ddo object (see help(object.ddo)).
X.SE
XNone.
X.DT
XInput parameters can be conceptually grouped as follows:
X DESIGN: xdesign
X EVALUATION: roixy, measurement.system, Xconceivable, CDname.store, droperr.
X SIMULATION CONTROL (GENERIC): mu, n, nforce, fam0, family, extend.info, randomseed.
X SIMULATION CONTROL (APPLICATION-SPECIFIC): ...
X OUTPUT CONTROL: main, tag, warn.
X.PP
XThe output class ddo object is structured as a list, but the list is usually
Xinvisible to the user because of the print.ddo method. See help(object.ddo)
Xfor more information. The elements of the list are:
X beta - the results of the simulation:
X CDall - a matrix of Calibration level discrepancy values
X DDall - a matrix of a few Design level discrepancies
X evalcall - an evaluated list of the elements (almost
Xall of them) used in the call or generated by default.
X call - the call to the function
X rseedend - the ending random seed
X error - exists only if failures to converge occurred.
X time - the amount of time taken for (i) simulation and (ii) evaluation.
X.PP
XIf failures to converge when fits are attempted, the error message is stored
Xand a table of the frequency of each error messages is assembled and stored
Xin the class ddo object which is output.
XIf nforce=T and more than n/2 failures to converge occur,
Xddomake() will give a progress report
Xas (number of failures to converge)/n passes each half-integer.
X.PP
XThe routine calls on several other routines for much of its work.
XIt calls on simulation.[application] to perform the simulations, then on
XdiscrepancyC.[application] (or discrepancyC.default) to obtain the
Xcalibration-level discrepancies and on discrepancyD.[application]
X(or discrepancyD.[default]) to obtain the design-level discrepancies.
X.SA
X make()
X extend()
X.EX
Xddomake(seq(6,8,length=14),n=10)
X #10 simulations of 14 points spaced
X #evenly over range from 6 to 8.
Xbegin("phcal")
Xddomake(seq(6,8,length=7),batch=2,main="two batches")
X #In "phcal" application, which includes
X #a "batch" input variable
X #(see simulate.phcal), do default number of simulations
X #(to query, type ddo("n")) of two batches of
X #7 points spread evenly over the range from 6 to 8.
X.KW ~keyword
X.WR
X
X
END_OF_FILE
if test 5004 -ne `wc -c <'.Data/.Help/ddomake'`; then
echo shar: \"'.Data/.Help/ddomake'\" unpacked with wrong size!
fi
# end of '.Data/.Help/ddomake'
fi
if test -f '.Data/.Help/extend' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'.Data/.Help/extend'\"
else
echo shar: Extracting \"'.Data/.Help/extend'\" \(2262 characters\)
sed "s/^X//" >'.Data/.Help/extend' <<'END_OF_FILE'
X.BG
X.FN extend
X.TL
XExtend or modify an existing class ddo object
X.DN
XAllows extension or modification of existing class ddo object by (1) changing n, the
Xnumber of simulations; (2) Changing the range of interest (in roixy - see below);
X(3) changing CDname.store, the set
Xof calibration discrepancies which are computed and stored;
X(4) changing nforce, whether all simulations are forced to
Xobtain convergence (see help(ddomake)); (5) changing the
Xmeasurement system (as expressed in measurement.system and Xconceivable);
X(6) changing droperr, whether or not errors are dropped from Design-level
Xdiscrepancy calculations. When changes are not specified in any of these arguments,
Xthe default is that their values are drawn from the existing object.
X.CS
Xextend(ddo, n=ddo$evalcall$n, nforce=ddo$evalcall$nforce,
X roixy=ddo$evalcall$roixy,
X CDname=ddo$evalcall$CDname.store,
X measurement.system=ddo$evalcall$measurement.system,
X Xconceivable=ddo$evalcall$Xconceivable,
X droperr=ddo$evalcall$droperr, newseed=F)
X.RA
X.AG ddo
XThe class ddo object to be extended or modified.
X.OA
X.AG n
XDesired number of simulations for new object. Must be equal or greater than existing number.
X.AG nforce
XWhether n successful (convergent) simulations should be forced, or only n trials made.
X.AG roixy
XThe range of interest, an output of roimake() (see help(roimake)).
X.AG CDname.store
XThe calibration-level discrepancies to be calculated and stored.
X.AG measurement.system
XThe measurement system (usually \"truncation\" or \"repeat.measures\").
X.AG Xconceivable
XThe conceivable range of X values used in measurement system.
X.AG droperr
XIf FALSE, errors are dropped from calculation of design-level discrepancies.
X.AG newseed
XIf FALSE, the random seed will be set to resume where the previous object left off.
XIf TRUE, whatever random seed is in memory currently will be used.
X.RT
XA class ddo object with the specified values.
X.SE
XNone.
X.SA
X ddomake()
X make()
X roimake()
X.EX
X candidate1aug <- extend(candidate1,n=candidate1["n"]+10)
X #adds 10 more simulations to class ddo obj. candidate1
X candidate1["CDall"] #examine the old
X candidate1aug["CDall"] #examine the new, see the overlap
X.KW ~keyword
X.WR
X
X
END_OF_FILE
if test 2262 -ne `wc -c <'.Data/.Help/extend'`; then
echo shar: \"'.Data/.Help/extend'\" unpacked with wrong size!
fi
# end of '.Data/.Help/extend'
fi
if test -f '.Data/.Help/focus' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'.Data/.Help/focus'\"
else
echo shar: Extracting \"'.Data/.Help/focus'\" \(1115 characters\)
sed "s/^X//" >'.Data/.Help/focus' <<'END_OF_FILE'
X.BG
X.FN focus
X.TL
XFocus on one discrepancy distribution object (ddo).
X.DN
XFunction to "focus in" on one particular class ddo object (the result of repeated
Xsimulations of one design) and plot a large batch of information about
Xthis object.
X.CS
Xfocus(ddo, ...)
X.RA
X.AG ddo
XThe discrepancy distribution object upon which we focus.
X.OA
X.AG ...
XAny arguments to be passed on to subordinate routines (showddo, plotdensity,
Xpairddo, or plotstrata).
X.RT
XNone.
X.SE
XProduces plots.
X.DT
XSequentially calls four routines:
X showddo() -plot the design
X plotdensity() -plot the densities of the
X calibration-level discrepancies
X pairddo() -scatterplots of pairs of
X calibration-level discrepancies
X plotstrata() -plots calibration curves, grouped by
X quantiles of calibration-level
X discrepancies
X.SA
X showddo()
X plotdensity()
X pairddo()
X plotstrata()
X.EX
Xfocus(candidates1,main="First Candidate")
X #get much information for class ddo
X #object candidates1;
X #Send main="First Candidate" argument to showddo()
X.KW ~keyword
X.WR
END_OF_FILE
if test 1115 -ne `wc -c <'.Data/.Help/focus'`; then
echo shar: \"'.Data/.Help/focus'\" unpacked with wrong size!
fi
# end of '.Data/.Help/focus'
fi
if test -f '.Data/.Help/index.ddo' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'.Data/.Help/index.ddo'\"
else
echo shar: Extracting \"'.Data/.Help/index.ddo'\" \(2889 characters\)
sed "s/^X//" >'.Data/.Help/index.ddo' <<'END_OF_FILE'
X.BG D
X.FN index.ddo
X.ce
XINDEX OF NONLINEAR CALIBRATION DESIGN FUNCTIONS
X
X.ce
XI. FUNCTIONS OF GENERAL INTEREST
X.PP
XThe following functions are all of interest to general users.
XEach has a help file.
X.IP ddobegin
X.PP
XInitialize a particular application
X.IP ddo
X.PP
XGeneric function to reset parameter values
X.IP intrinsic
X.PP
XCalculates and plots intrinsic aberrance, a linear approximation
Xfor pointwise discrepancy,
Xand other quantities of interest which are intrinsic to
Xthe entire measurement system.
X.IP ddomake
X.PP
XCreate an estimated distribution, a "discrepancy distribution object
X.IP roimake
X.PP
XCreate a range of interest as expressed using either menu name,
Xx vector, or y vector; weights optional
X.IP compareddo
X.PP
XCompares several distribution objects [and called by plot(many args)].
X.IP plotddo
X.PP
XCombines showddo() and plotdensity() [and called by plot(one arg)].
X.IP showddo
X.PP
XPlots a design and some summary statistics from a ddo.
X.IP focus
X.PP
XRuns several analyses on one design: plotdensity, plotstrata, pairddo
X.IP plotdensity
X.PP
XPlots discrepancy distribution from discrepancy distribution object
X.IP plotstrata
X.PP
XPlot discrepancy curves, stratified by quantiles
X.IP pairddo
X.PP
XPlots correlations of pairs of discrepancy distributions
X.IP mwtest
X.PP
XPerforms a Mann-Whitney rank sum test comparison between
Xtwo discrepancy distribution objects
X.IP extend
X.PP
XExtends an existing ddo to more trials, more discrepancy
Xfunctions, or a new range of interest
X.IP print.ddo
X.PP
X[can be accessed by print() or by typing name of object]
XFunction to print a discrepancy distribution object in
Xa readable output
X.IP plot.ddo [can be accessed by plot()]
X.PP
XCalls either plotddo (if one ddo argument) or compareddo(two
Xor more ddo arguments)
X.IP make
X.PP
XTakes a list of (named) designs and creates a group of ddo's.
X
X
X.ce
X OTHER ONLINE HELP
X
X.IP print(ddo)
X.PP
XTyping "print(ddo)" or simply "ddo" prints out an
Xoverall introduction to the Nonlinear Calibration Design
Xsoftware and includes a reference to literature
Xdescribing the theoretical basis for the measures
Xused by the software.
X
X.ce
XOther help files include the following:
X.IP application
X.PP
XDescribes how to construct a new application such as the
XpH calibration application described in help(phcal) and
Xin the original paper (see print(ddo)).
X.IP ddoframe
X.PP
XDescribes the class ddoframe objects and the
Xassociated data management functions.
X.IP object.ddo
X.PP
XDescribes the contents of a class ddo object
X.IP phcal
XDescribes the phcal application.
X
X
X.ce
XOTHER DOCUMENTATION
X
X.PP
XThere is substantial documentation in the comments in the
Xshareware. Much of this commentary is embedded in the
Xroutines but some is contained in the header to
Xthe entire file in ascii form.
X.KW sysdata
X.WR
X
END_OF_FILE
if test 2889 -ne `wc -c <'.Data/.Help/index.ddo'`; then
echo shar: \"'.Data/.Help/index.ddo'\" unpacked with wrong size!
fi
# end of '.Data/.Help/index.ddo'
fi
if test -f '.Data/.Help/intrinsic' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'.Data/.Help/intrinsic'\"
else
echo shar: Extracting \"'.Data/.Help/intrinsic'\" \(1765 characters\)
sed "s/^X//" >'.Data/.Help/intrinsic' <<'END_OF_FILE'
X.BG
X.FN intrinsic
X.TL
XFunction to display Intrinsic Aberrance and Linearly Approximated Pointwise Discrepancy
X.DN
XPlots graphs of Point-level Intrinsic Aberrance versus X for both absolute
Xand truncation-aberrance, also displaying the maxima of these curves
Xwithin range of interest, which is the Calibration-level Intrinsic Aberrance.
XFor Linearly Approximated Pointwise Discrepancy, plots PD versus X, and
Xdisplays ratio of max to min over the range of interest.
X.CS
Xintrinsic(
X xcon = list(absolute = c( - Inf, Inf),
X truncation = .DDO$Xconceivable),
X mu = .DDO$mu, roi = .DDO$roi,
X roixy = roimake(name = roi, mu = mu, fam0 = fam0),
X xplotlim = <>, xplot = <>,
X epsilonpart = 0.001, sigma = .DDO$sigma0,
X fam0 = .DDO$fam0, doplot = T)
X.RA
X.OA
X.AG xcon
XThe conceivable range (used for calculation of truncation aberrance)
X.AG mu
XThe imagined true curve
X.AG roi
XA name for the range of interest
X.AG roixy
XA range of interest list, an output of roimake().
X.AG xplotlim
XThe X plotting limits desired
X.AG xplot
XThe vector of X values to be plotted
X.AG epsilonpart
XThe size of the increments used for calculating the linear approximation
X.AG sigma
XThe imagined true Y measurement error
X.AG fam0
XThe imagined true family of curves
X.AG doplot
XIf FALSE, a plot will not be done (used by pairddo())
X.RT
XNone unless doplot=F, in which case a list with two elements is returned:
XPA = intrinsic aberrance, a vector same length as xplot;
XXse = linearly approximated pointwise discrepancy of Xmeasured, also a
Xvector same length as xplot.
X.SE
XMakes plots as described.
X.SA
X pairddo()
X.EX
X intrinsic() #makes the two graphs described
X intrinsic(xplotlim=c(4,15)) # same graphs but with
X #specified plotting limits
X.KW ~keyword
X.WR
END_OF_FILE
if test 1765 -ne `wc -c <'.Data/.Help/intrinsic'`; then
echo shar: \"'.Data/.Help/intrinsic'\" unpacked with wrong size!
fi
# end of '.Data/.Help/intrinsic'
fi
if test -f '.Data/.Help/make' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'.Data/.Help/make'\"
else
echo shar: Extracting \"'.Data/.Help/make'\" \(3164 characters\)
sed "s/^X//" >'.Data/.Help/make' <<'END_OF_FILE'
X.BG
X.FN make
X.TL
XFunction to create simulations for a list of designs. Creates a ddoframe object.
X.DN
XTakes designs given in second (makelist) argument, performs simulations on them
Xas specified in third (n) argument (and in default parameters specified in ddo()),
Xand stores the resulting class ddo objects in a ddoframe object with name
Xspecified in the first argument (prefix). The individual class ddo objects are stored
Xin objects with names which are the first argument (prefix)
Xpasted to the counting numbers (i.e., if prefix="pr",
Xthen objects are pr1, pr2, pr3, ... and the ddoframe class object which connects them is pr).
X.CS
Xmake(prefix=., makelist=., n=.DDO$n,
X override=F, where=1, ...)
X.RA
X.AG prefix
XContains the name of the resulting ddoframe object. It is also the
Xprefix to each of the resulting class ddo objects which contain the simulations.
X.AG makelist
XContains a list of the designs to be simulated. The simulations are performed by
Xcalling ddomake() for each design. The elements of names(makelist)
Xare used as the title (main=) argument to ddomake. Each element of makelist can be either
Xa vector or a list in which all elements except the first must be named.
XIf it is a vector, the vector is sent as the xdesign argument
Xto ddomake(). If it is a named list, the first argument should be a vector to be used
Xas the xdesign argument to ddomake() and the subsequent named elements should
Xbe other arguments to be sent to ddomake().
X.OA
X.AG n
XThe number of simulations to be used for each design (unless overridden by a
Xnamed list as an element in makelist). Passed as n argument to ddomake().
X.AG override
XIf FALSE, make() will refuse to overwrite preexisting objects.
X.AG where
XThe data frame in which to store the resulting class ddo and class ddoframe objects.
XPassed as an argument to assign().
X.AG ...
XAny other arguments to be passed on to ddomake().
X.RT
XNone.
X.SE
XPrints name, date and time of creation of each created object, and stores
Xcreated objects in data base specified by where.
X.DT
XIf the function receives inadequate arguments, it will store an NA in
Xthe corresponding class ddo object.
X.SH WARNING
XCurrently does not check to see if there is enough space in memory.
X.SA
X ddomake() #creates a single class ddo object
X extend() #updates a class ddo object
X ddoframe #class ddoframe objects consist of many
X #ddo objects assembled together
X ddo() #sets default parameter values used by make()
X.EX
Xmake("candidates",list(twobatch=seq(from=7,to=8,length=14),
X threebatch=seq(from=7,to=8,length=21)),n=10)
X #creates class ddoframe object candidates and
X #class ddo objects candidates1 and candidates2,
X #each the result of 10 simulations
Xmake("canmore",list("5 batches"=seq(from=7,to=8,length=35),
X newcon=list(candidates1["xdesign"],
X Xconceivable=c(5,9))))
X #creates class ddoframe object canmore and
X #class ddo objects canmore1 and canmore2.
X #The design for canmore2 is the same as
X #the design for candidates1;
X #the default number of repetitions is
X #used (see help(ddo)).
X.KW ~keyword
X.WR
X
END_OF_FILE
if test 3164 -ne `wc -c <'.Data/.Help/make'`; then
echo shar: \"'.Data/.Help/make'\" unpacked with wrong size!
fi
# end of '.Data/.Help/make'
fi
if test -f '.Data/.Help/mwtest' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'.Data/.Help/mwtest'\"
else
echo shar: Extracting \"'.Data/.Help/mwtest'\" \(1145 characters\)
sed "s/^X//" >'.Data/.Help/mwtest' <<'END_OF_FILE'
X.BG
X.FN mwtest
X.TL
XFunction to make Mann-Whitney comparison between two class ddo objects.
X.DN
XPerforms Mann-Whitney comparison between class ddo objects dist1 and
Xdist2 for Calibration-level discrepancies named in CDname.
X.CS
Xmwtest(dist1, dist2, CDname=.DDO$CDname.store, droperr=F)
X.RA
X.AG dist1
XThe first class ddo object
X.AG dist2
XThe second class ddo object
X.OA
X.AG CDname
XThe names of the Calibration-level discrepancies to be compared
X.AG droperr
XIf TRUE, errors (nonconvergences) are ignored; if FALSE, they are
Xtreated as having infinite discrepancy.
X.RT
XA list with three entries:
X.IP mw
Xa vector of length=length(CDname) with each entry being
XP[X'.Data/.Help/object.ddo' <<'END_OF_FILE'
X.BG D
X.FN object.ddo
X.TL
XClass ddo objects.
X.PP
XClass ddo objects are the core of the ddo (Discrepancy
XDistribution Object) software system for calibration
Xdesign. There are routines for (I) creating them singly
X(ddomake) or in groups (make), for (II) viewing them singly
X(plot, focus, showddo, plotdensity, plotstrata, pairddo)
Xor in groups (plot, compareddo), or (III) printing out information
Xsingly (print, names, "[") or for pairs (mwtest). These objects may
Xbe combined together into ("virtual") class ddoframe objects
Xfor convenient manipulation and viewing (c).
X.PP
XTheir elements
Xmay also be accessed individually by methods for "[" and "[<-".
XThe elements of a class ddo object fall into the following major
Xgroups (some of which only contain one element):
X.IP beta
X.PP
XA matrix containing the simulated fitted parameter values.
X.IP CDall
X.PP
XA matrix containing the Calibration-level discrepancies for each simulation.
X.IP DDall
X.PP
XA matrix of Design-level summary statistics.
X.IP rseedend
X.PP
XThe value of the random seed when the simulations were finished.
X.IP call
X.PP
XThe call which produced the object (if it was updated, a list of calls).
X.IP evalcall
X.PP
XAll the arguments which produced the simulation, as they
Xwere evaluated when received in the call. The evaluated argument,
Xof course, may simply be a
Xcopy of the relevant default parameter (from .DDO)
Xat the time of the call to ddomake().
XThe arguments which are stored in the class ddo object
Xinclude the following parameters common to all applications:
Xxdesign, mu, roixy, measurement.system, Xconceivable, CDname.store,
Xdroperr, n, nforce, fam0, family, randomseed, and main.
XAlso stored are any application-specific parameters.
XFor the phcal application, these application-specific
Xparameters are batch, N, sigma, sigmabatch, and start
X[see help(application) and help(phcal)].
X.IP error
X.PP
XExists only if failures to converge occurred. Contains
Xinformation about the failures to converge - generally
Xa table giving counts of
Xwhat messages were sent regarding the type of failure to
Xconverge
X.IP time
X.PP
XThe amount of time taken for (i) simulation and (ii)
Xevaluation. Somewhat system-dependent.
X.KW sysdata
X.WR
END_OF_FILE
if test 2199 -ne `wc -c <'.Data/.Help/object.ddo'`; then
echo shar: \"'.Data/.Help/object.ddo'\" unpacked with wrong size!
fi
# end of '.Data/.Help/object.ddo'
fi
if test -f '.Data/.Help/pairddo' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'.Data/.Help/pairddo'\"
else
echo shar: Extracting \"'.Data/.Help/pairddo'\" \(1007 characters\)
sed "s/^X//" >'.Data/.Help/pairddo' <<'END_OF_FILE'
X.BG
X.FN pairddo
X.TL
XFunction to do pairwise scatterplots of simulated calibration discrepancies.
X.DN
XAll scatterplots are on one page, as with the output of the pairs() function in Splus.
X.CS
Xpairddo(ddo,...)
X.RA
X.AG ddo
XThe class ddo object.
X.OA
X.AG ...
XSeveral other named arguments exist, but are of little interest. CAintrinsic is defined
Xas an expression which computes intrinsic aberrance. Other arguments are inherited from pairs().
X.RT
XNone.
X.SE
XCreates pairwise scatterplots.
X.DT
XThis function is an adaptation of the pairs() function in Splus. It does pairwise
Xscatterplots of the calibration discrepancies (including calibration-level
Xaberrance) in a class ddo object (the discrepancy values are contained
Xin object.ddo["CDall"]. It also plots the intrinsic aberrance so that the
Xcentering of the calibration-level aberrance around intrinsic aberrance
Xcan be monitored.
X.SA
X focus()
X plotstrata()
X.EX
Xpairddo(candidate1)
X #plots scatterplots for class ddo object candidate1
X.KW ~keyword
X.WR
X
END_OF_FILE
if test 1007 -ne `wc -c <'.Data/.Help/pairddo'`; then
echo shar: \"'.Data/.Help/pairddo'\" unpacked with wrong size!
fi
# end of '.Data/.Help/pairddo'
fi
if test -f '.Data/.Help/phcal' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'.Data/.Help/phcal'\"
else
echo shar: Extracting \"'.Data/.Help/phcal'\" \(2135 characters\)
sed "s/^X//" >'.Data/.Help/phcal' <<'END_OF_FILE'
X.BG D
X.FN phcal
X.TL
XThe "phcal" application for the Nonlinear Calibration Design software
X.PP
XThe Nonlinear Calibration Design software was developed
Xat U.C. Berkeley in
Xconjunction with an optometric application. This calibration design
Xproblem is described more extensively in the reference given in
Xprint(ddo) and involved calibration using a three-parameter
Xlogistic curve to relate pH level (the X variable)
Xwith fluorescein intensity ratio, or fir (the Y variable).
X.PP
XThe phcal application is the only example of an application that is
Xbundled with the software, and it is the default application when
Xddobegin() is invoked. Users wishing to implement their own
Xapplications should read help(application) and may want to examine
Xthe phcal functions and objects. Some of the functions are the
Xphcal versions of functions required for all applications, some
Xare optional (but interfaced with the shareware), and some were
Xsimply autonomous and idiosyncratic to the phcal application. For example,
Xdesignphcal() was used for generating design points for exploration
Xof heuristic notions for experimental design (See Appendix I
Xof the Master's thesis by Doug Oman referenced in print(ddo)).
X
X.ce
XPH-CALIBRATION SPECIFIC FUNCTIONS AND OBJECTS:
X
X.ce
X(REQUIRED FUNCTIONS)
X.IP ddo.phcal
X.PP
Xfunction to initialize (or reset) phcal values
X.IP yofx.phcal
X.PP
Xfind fir as a function of ph, beta
X.IP xofy.phcal
X.PP
Xfind ph as a function of fir, beta
X.IP simulate.phcal
X.PP
Xfunction to simulate a calibration experiment
X
X.ce
X(OPTIONAL FUNCTIONS)
X.IP print.DDO.phcal
X.PP
Xaugments printDDO by printing out application-specific
Xdefault parameters.
X.IP printdetail.phcal
X.PP
Xaugments generic print.ddo function with phcal-specific
Xmethods.
X
X.ce
X(USED BY OTHER PHCAL-SPECIFIC; or AUTONOMOUS)
X.IP sigmaiidphcal
X.PP
Xcalculate the overall residual sigma assuming iid.
X(used by discrepancy.phcal)
X.IP designphcal
X.PP
Xinitializes ph-calibration design tools
X
X.ce
X(OBJECT)
X.IP phrab.bestfit1
X.PP
Xhelped supply default parameters for mu, sigma0, roi
X.KW sysdata
X.WR
X
END_OF_FILE
if test 2135 -ne `wc -c <'.Data/.Help/phcal'`; then
echo shar: \"'.Data/.Help/phcal'\" unpacked with wrong size!
fi
# end of '.Data/.Help/phcal'
fi
if test -f '.Data/.Help/plot.ddo' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'.Data/.Help/plot.ddo'\"
else
echo shar: Extracting \"'.Data/.Help/plot.ddo'\" \(1144 characters\)
sed "s/^X//" >'.Data/.Help/plot.ddo' <<'END_OF_FILE'
X.BG
X.FN plot.ddo
X.FN plot.ddoframe
X.TL
XPlotting methods for class ddo and class ddoframe objects
X.DN
XFor one class ddo argument, plot basic information with plotddo(); for
Xmultiple class ddo objects, do comparison analyses through compareddo().
X.CS
Xplot(...), plot.ddo(...) or plot.ddoframe(...)
X.RA
XAt least one class ddo object,
Xperhaps contained in a class ddoframe object.
X.OA
X.AG ...
XAuxiliary arguments to be passed along to plotddo() or compareddo().
X.RT
XNone.
X.SE
XPlots as described below.
X.DT
XThe plot method for class ddo or ddoframe is the same.
XAll arguments are examined and separated into those comprising
Xclass ddo or ddoframe objects, and all others. Then either
Xplotddo() or compareddo() is called, depending upon whether there
Xis only one (total) ddo object, or many (including all those contained
Xin ddoframe objects). The other
Xarguments are then passed along as auxiliary arguments to
Xplotddo() or compareddo().
X.SA
X plotddo()
X compareddo()
X.EX
X plot(ddo1) #same as plotddo(ddo1)
X plot(CDname="CD1",c(ddo1,ddo2),anal="f")
X #plots CD1 densities for ddo1 and ddo2
X.KW ~keyword
X.WR
X
END_OF_FILE
if test 1144 -ne `wc -c <'.Data/.Help/plot.ddo'`; then
echo shar: \"'.Data/.Help/plot.ddo'\" unpacked with wrong size!
fi
# end of '.Data/.Help/plot.ddo'
fi
if test -f '.Data/.Help/plotddo' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'.Data/.Help/plotddo'\"
else
echo shar: Extracting \"'.Data/.Help/plotddo'\" \(1192 characters\)
sed "s/^X//" >'.Data/.Help/plotddo' <<'END_OF_FILE'
X.BG
X.FN plotddo
X.TL
XFunction to plot basic summary information for a class ddo object
X.DN
XPlots the design and discrepancy densities of a class
Xddo object, and also some summary numerical information.
X.CS
Xplotddo(...)
X.RA
XThe first argument should contain the class ddo object.
X.OA
X.AG ...
XIt can be supplied with auxiliary arguments used by plotdensity() or showxddo().
X.RT
XNone.
X.SE
XA single page plot with the design and numerical information in a large
Xplot in the upper left corner, and the Calibration-level discrepancy
Xplots arranged below and to the right.
X.DT
XThe plotddo() function is essentially a combination of showddo() and
Xplotdensity(), which it calls in sequence. The plot of the design and
Xthe associated numerical summary information (from ddo["DDall"],
Xwhere ddo is the class ddo object) is generated by showddo() and is
Xmade larger than the other plots so that the numerical information will
Xbe readable.
X.SA
X plot.ddo()
X plotdensity()
X.EX
X plotddo(ddo1) #displays information for
X #class ddo object ddo1
X plotddo(ddo1,CDname=c("CD1","CA"))
X #displays design, summary info,
X #but densities only for CD1 and CA.
X.KW ~keyword
X.WR
X
END_OF_FILE
if test 1192 -ne `wc -c <'.Data/.Help/plotddo'`; then
echo shar: \"'.Data/.Help/plotddo'\" unpacked with wrong size!
fi
# end of '.Data/.Help/plotddo'
fi
if test -f '.Data/.Help/plotstrata' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'.Data/.Help/plotstrata'\"
else
echo shar: Extracting \"'.Data/.Help/plotstrata'\" \(1900 characters\)
sed "s/^X//" >'.Data/.Help/plotstrata' <<'END_OF_FILE'
X.BG
X.FN plotstrata
X.TL
XPlot simulated calibration curves from a class ddo object, stratified by discrepancy level.
X.DN
XFor each type of calibration discrepancy (CD) stored, stratifies calibration curves stored
Xin the class ddo object by quantile of CD. Then outputs a plot of POINTWISE DISCREPANCY
Xat X (that is, median(Xmeasured) minus Xtrue) or POINTWISE ABERRANCE
Xversus X for each curve in the stratum.
X.CS
Xplotstrata(ddo,
X CDname=ddo$evalcall$CDname.store, PDname=<>,
X qua=c(0.05, 0.1, 0.25, 0.5, 0.75, 0.9, 0.95),
X xplotlim=<>, xplot=<>,
X xcon=<>,...)
X.RA
X.AG ddo
XA class ddo object
X.OA
X.AG CDname
XNames of the Calibration Discrepancy strata to be plotted, if possible.
X.AG PDname
XA character vector the same length as the CDname vector which
Xspecifies whether pointwise discrepancy (PD) or pointwise aberrance (PA)
Xshould be plotted.
X.AG qua
XThe quantiles which form the cutpoints for constructing the CD strata.
X.AG xplotlim
XThe min and max of the X points to be plotted. A default exists.
X.AG xplot
XThe vector of X points to be plotted. A default exists.
X.AG xcon
XThe conceivable ranges for X. Should be a list of form
Xlist(absolute=c(xmin1,xmax1),truncation=c(xmin2,xmax2)) with
Xa named element for each measurement system of interest.
XA default exists.
X.AG ...
XExtra arguments are ignored.
X.RT
XNone.
X.SE
XPlots a separate page for each CD, with a separate frame for each stratum.
X.DT
XOften there seems to be strong similarities between the stratifications from
Xdifferent CDs (e.g., CD1, CD2, CDinf), as would be expected.
X.SA
X focus()
X pairddo()
X.EX
X plotstrata(candidate1)
X #Plots stratifications by all stored
X #Calibration-level discrepancies.
X plotstrata(candidate1,CDname="CD1")
X #Plots stratifications only by CD1 (mean
X #absolute deviation).
X.KW ~keyword
X.WR
X
X
END_OF_FILE
if test 1900 -ne `wc -c <'.Data/.Help/plotstrata'`; then
echo shar: \"'.Data/.Help/plotstrata'\" unpacked with wrong size!
fi
# end of '.Data/.Help/plotstrata'
fi
if test -f '.Data/.Help/print.ddo' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'.Data/.Help/print.ddo'\"
else
echo shar: Extracting \"'.Data/.Help/print.ddo'\" \(1273 characters\)
sed "s/^X//" >'.Data/.Help/print.ddo' <<'END_OF_FILE'
X.BG
X.FN print.ddo
X.TL
XPrint method for a class ddo (Discrepancy Distribution) object
X.DN
XPrints most elements in a class ddo object which have names given
Xby names(object.ddo), excepting the large ones such as
Xbeta and CDall. Gives a short description of each element
Xand categorizes according to conceptual type (simulation control,
Xevaluation control, etc). Also prints some summary statistics,
Xincluding a confidence level for
Xthe frequency of errors (i.e., nonconvergences)
X.CS
Xprint.ddo(ddo, CLerr=0.95, quote, prefix)
X.RA
X.AG ddo
XThe class ddo object to be printed
X.OA
X.AG CLerr
XConfidence level for the confidence interval to be printed for
Xthe chance of nonconvergence
X.AG quote
X.AG prefix
XThese two arguments are needed
Xso that print.list won't freak out when a ddo is an
Xelement of a list.
X.RT
XThe class ddo object itself.
X.SE
XPrints the ddo.
X.DT
XCalls upon print.detail.[application name] if it exists,
Xin order to print out application-specific parameters and
Xarguments.
X.SH NOTE
XCLerr is not strictly a confidence level because
Xpointwise convergence does not occur near 0 and 1.
X.SA
X application
X.EX
X ddo1 #prints class ddo object ddo1
X print(ddo1,CLerr=.99) #Prints object with
X #.99 "Confidence Levels"
X #for chance of nonconvergence.
X.KW ~keyword
X.WR
END_OF_FILE
if test 1273 -ne `wc -c <'.Data/.Help/print.ddo'`; then
echo shar: \"'.Data/.Help/print.ddo'\" unpacked with wrong size!
fi
# end of '.Data/.Help/print.ddo'
fi
if test -f '.Data/.Help/roimake' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'.Data/.Help/roimake'\"
else
echo shar: Extracting \"'.Data/.Help/roimake'\" \(2506 characters\)
sed "s/^X//" >'.Data/.Help/roimake' <<'END_OF_FILE'
X.BG
X.FN roimake
X.TL
XFunction to construct a list containing complete information
Xabout the range of interest (ROI).
X.DN
XFunction roimake() constructs a list containing the following four elements:
X(1) name of the range of interest (2) ROI expressed on X scale (3) ROI
Xexpressed on Y scale (via Y=F_mu(X)) (4) weights of each of the points ROI.
XThe latter three elements are numeric vectors which are used for numerical
Xintegrations over the range of interest. The output of roimake() is used as
Xan input to ddomake().
X.CS
Xroimake(mu=.DDO$mu, fam0=.DDO$fam0, ...)
X.RA
XRequired input may be given in three possible forms. The
Xthree forms, and the arguments to which they are supplied,
Xare as follows:
X.AG name
Xlooks up arguments in .DDO$roimenu corresponding
Xto the name. These arguments are appended to
Xthe argument list received by this function.
X(And, since an x vector is sufficient to generate
Xthe entire object, no other arguments need be given
Xif, as is usually the case, the menu item is of the
Xform name=list(x=c(...)) ).
X.AG "x or y"
XIf either x or y vector missing, it will be calculated
Xfrom the other using mu and fam0. At least one of
Xx or y must be supplied (either directly or though using .DDO$roimenu).
X.AG nwxy
XA complete roimake() output object which will simply be re-output.
X.OA
X.AG w
XRelative weights. Repeated to make vector the
Xsame length as x or y. Default is w=1. Is then normalized
Xto sum to 1 before being put in output object.
X.AG mu
XThe parameters of the "imagined true" calibration curve.
X.AG fam0
XThe "imagined true" family of calibration curves.
X.RT
XA list with the following names and meanings:
X.IP name
Xname argument, or else (the default value) a constructed name of the
Xform "n points X in [xmin,xmax]".
X.IP x
Xroi expressed on X scale
X.IP y
Xroi expressed on Y scale
X.IP w
Xrelative weights of different objects (normalized to sum to 1).
X.DT
XOther than mu and fam0, all arguments are input through a "..." argument. Note also
Xthat the function uses default values for mu and fam0 from default parameters;
Xthus, evaluating the same roimake() expression when these values have changed
Xresults in a different output.
X.SA
X ddomake()
X ddo()
X.EX
X begin("phcal")
X roimake(name="fall94") #gets ROI specified in the
X #menu specified in ddo.phcal
X roimake(x=seq(6.8,8.0,length=30),w=1+((1:30)>15))
X #uses 30 points in [6.8,8.0] with the
X #larger half of range weighted double
X.KW ~keyword
X.WR
X
X
END_OF_FILE
if test 2506 -ne `wc -c <'.Data/.Help/roimake'`; then
echo shar: \"'.Data/.Help/roimake'\" unpacked with wrong size!
fi
# end of '.Data/.Help/roimake'
fi
if test -f '.Data/.Help/showddo' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'.Data/.Help/showddo'\"
else
echo shar: Extracting \"'.Data/.Help/showddo'\" \(535 characters\)
sed "s/^X//" >'.Data/.Help/showddo' <<'END_OF_FILE'
X.BG
X.FN showddo
X.TL
XFunction to plot the design of a class ddo object, and basic information
X.DN
XSame as plotddo(), but no Calibration-level discrepancy densities are plotted.
X.CS
Xshowddo(ddo, ...)
X.RA
X.AG ddo
XThe class ddo object to be examined.
X.OA
X.AG ...
XUninteresting to most users; passed on to showxddo()
X.RT
XNone.
X.SE
XGenerates plot described
X.DT
XPrints the numerical information in ddo["DDall"] on the face of a plot
Xof the design.
X.SA
X plotddo()
X compareddo()
X.EX
X showddo(ddo1) #shows class ddo object ddo1
X.KW ~keyword
X.WR
END_OF_FILE
if test 535 -ne `wc -c <'.Data/.Help/showddo'`; then
echo shar: \"'.Data/.Help/showddo'\" unpacked with wrong size!
fi
# end of '.Data/.Help/showddo'
fi
echo shar: End of shell archive.
echo "Now get into S and type source(''DDOCODE'')"
exit 0