; XLISP-STAT code for the *QV Calculator*
;
; (including the online QV calculator at
; http://www.stats.ox.ac.uk/~firth/qvcalc/)
;
; Copyright David Firth
; Version 1.1, May 2000
(provide "qvcalc")
;
; Set up a GLIM model prototype for the "lognormal" criterion function used in
; #'MAKE-QV-OBJECT to compute quasi-variances
;
(require "glim")
;;;; Normal-error model with exp link prototype
(defproto normal-exp-proto () () glim-proto)
;;;; LINK FUNCTION (exponential)
(defmeth normal-exp-proto :initialize-search (&optional eta)
(call-next-method (if eta eta (exp (send self :y)))))
(defmeth normal-exp-proto :fit-means
(&optional (eta (send self :eta))) (log eta))
(defmeth normal-exp-proto :fit-link-derivs (mu)
(exp mu))
;;;; VARIANCE FUNCTION and DEVIANCE (constant, and squared error)
(defmeth normal-exp-proto :fit-variances (mu)
(repeat 1 (length mu)))
(defmeth normal-exp-proto :fit-deviances (mu)
(flet ((log+ (x) (log (if-else (< 0 x) x 1)))) ; to prevent log of zero
(let* ((y (send self :yvar))
(pw (send self :pweights))
(raw-dev (^ (- y mu) 2)))
(if pw (* raw-dev pw) raw-dev))))
;;;; Model constructor function
(defun normal-exp-model (x y &rest args)
(apply #'send normal-exp-proto :new :x x :y y args))
(send normal-exp-proto :epsilon 0.00000001) ; use strict
(send normal-exp-proto :epsilon-dev 0.00000001) ; convergence criteria
;
; Define an object prototype that will hold quasi-variances and
; associated methods
;
(defproto qv-proto ; the quasi-variance object prototype
'(model-name factor-name levels level-names response-name
covariance-matrix quasi-variances simple-contrast-labels
simple-contrast-variances approximated-simple-contrast-variances)
()
*object*
"Quasi-variances for the coefficients associated with a factor
in a regression model")
;
; Accessor methods for all slots in QV-PROTO
;
(defmeth qv-proto :model-name (&optional (model-name nil set))
"Message args: (&optional (model-name nil set))
Sets or returns a description of the model to which the
quasi-variances relate."
(if set (setf (slot-value 'model-name) model-name))
(slot-value 'model-name))
(defmeth qv-proto :factor-name (&optional (factor-name nil set))
"Message args: (&optional (factor-name nil set))
Sets or returns a description of the factor to which the
quasi-variances relate."
(if set (setf (slot-value 'factor-name) factor-name))
(slot-value 'factor-name))
(defmeth qv-proto :levels (&optional (levels nil set))
"Message args: (&optional (levels nil set))
Sets or returns the number of levels of the factor to which the
quasi-variances relate."
(if set (setf (slot-value 'levels) levels))
(slot-value 'levels))
(defmeth qv-proto :level-names (&optional (level-names nil set))
"Message args: (&optional (level-names nil set))
Sets or returns the level labels (a list of strings) for the factor
to which the quasi-variances relate."
(if set (setf (slot-value 'level-names) level-names))
(slot-value 'level-names))
(defmeth qv-proto :response-name (&optional (response-name nil set))
"Message args: (&optional (response-name nil set))
Sets or returns the name of the response variable to which the
quasi-variances relate."
(if set (setf (slot-value 'response-name) response-name))
(slot-value 'response-name))
(defmeth qv-proto :covariance-matrix (&optional (covariance-matrix nil set))
"Message args: (&optional (covariance-matrix nil set))
Sets or returns the covariance matrix from which the quasi-variances
are derived."
(if set (setf (slot-value 'covariance-matrix) covariance-matrix))
(slot-value 'covariance-matrix))
(defmeth qv-proto :quasi-variances (&optional (quasi-variances nil set))
"Message args: (&optional (quasi-variances nil set))
Sets or returns the quasi-variances as a numeric list."
(if set (setf (slot-value 'quasi-variances) quasi-variances))
(slot-value 'quasi-variances))
(defmeth qv-proto :simple-contrast-variances
(&optional (simple-contrast-variances nil set))
"Message args: (&optional (simple-contrast-variances nil set))
Sets or returns the variances of all simple contrasts, as a numeric list."
(if set (setf (slot-value 'simple-contrast-variances)
simple-contrast-variances))
(slot-value 'simple-contrast-variances))
(defmeth qv-proto :approximated-simple-contrast-variances
(&optional (approximated-simple-contrast-variances nil set))
"Message args: (&optional (approximated-simple-contrast-variances nil set))
Sets or returns the QV-approximated variances of
all simple contrasts, as a numeric list."
(if set (setf (slot-value 'approximated-simple-contrast-variances)
approximated-simple-contrast-variances))
(slot-value 'approximated-simple-contrast-variances))
(defmeth qv-proto :simple-contrast-labels
(&optional (simple-contrast-labels nil set))
"Message args: (&optional (simple-contrast-labels nil set))
Sets or returns a list of labels for all simple contrasts,
in the order used for slots SIMPLE-CONTRAST-VARIANCES and
APPROXIMATED-SIMPLE-CONTRAST-VARIANCES."
(if set (setf (slot-value 'simple-contrast-labels)
simple-contrast-labels))
(slot-value 'simple-contrast-labels))
(defmeth qv-proto :summary (&optional (as-string nil))
"Message args: (&optional (as-string nil))
If AS-STRING is non-nil, returns a string which summarizes the
quasi-variance approximation to the covariance matrix of coefficients
associated with a factor. Otherwise prints that summary to standard output
and returns NIL."
(let* ((relerrs ; relative errors for simple contrasts
(- (sqrt
(/ (send self :approximated-simple-contrast-variances)
(send self :simple-contrast-variances)))
1))
)
(format (not as-string)
(format-output
(send self :model-name)
(send self :response-name)
(send self :factor-name)
(send self :covariance-matrix)
(send self :levels)
(send self :level-names)
(send self :quasi-variances)
(* 100 (min relerrs))
(* 100 (max relerrs))
(* 100 (- (first (send self :worst-errors)) 1))
(* 100 (- (second (send self :worst-errors)) 1))))))
(defmeth qv-proto :worst-errors ()
"Message args ()
Computes the worst error incurred by the quasi-variance approximation,
in the space of all possible contrasts. The result is a 2-element list
containing the smallest and largest ratios of approximate to exact
standard error."
(let* ((cvmat (send self :covariance-matrix))
(qvmat (diagonal (send self :quasi-variances)))
(r-form (reduced-form cvmat qvmat))
(r-cvmat (first r-form))
(r-qvmat (second r-form))
(inverse-sqrt (inverse (first (chol-decomp r-cvmat))))
(evalues (eigenvalues
(matmult inverse-sqrt r-qvmat (transpose inverse-sqrt))))
)
(sqrt (list (min evalues) (max evalues)))))
(defun reduced-form (cvmat qvmat)
; called by method :WORST-ERRORS
(let* ((levels (first (array-dimensions cvmat)))
(first-row (coerce (first (row-list cvmat)) 'list))
(ones (repeat 1 levels))
(j (outer-product ones ones))
(not-zero (iseq 1 (- levels 1)))
(r-cvmat (select
(+ cvmat (* (first first-row) j)
(- (+ (outer-product first-row ones)
(outer-product ones first-row))))
not-zero not-zero))
(qv1 (select qvmat 0 0))
(r-qvmat (select (+ qvmat (* qv1 j)) not-zero not-zero))
)
(list r-cvmat r-qvmat)))
(defmeth qv-proto :contrast-variance (c &key (print t))
"Message args: (c &key (print t))
Computes the actual and approximated contrast variance for any
specified contrast, C.
Prints a summary if :print is non-nil, and returns a 2-element list."
(cond ((not (eq (length c) (send self :levels)))
(error (format nil
"Argument to method :CONTRAST-VARIANCE should be a list of ~
length ~,2d" (send self :levels))))
((not (< (abs (sum (/ c (sum (abs c))))) 0.01))
(error (format nil
"Argument to method :CONTRAST-VARIANCE does not appear ~
to be a contrast (its elements do not sum to 0).")))
(t (let* ((cvmat (send self :covariance-matrix))
(qvmat (diagonal (send self :quasi-variances)))
(result (list (matmult c cvmat c) (matmult c qvmat c)))
)
(cond
(print
(princ "Contrast coefficients: ")
(princ c) (terpri)
(princ "Contrast variance: ")
(format t "~5,3g" (first result)) (terpri)
(princ "QV approximation: ")
(format t "~5,3g" (second result)) (terpri)
result)
(t result))))))
(defun qv-objectp (symbol)
"Args: (symbol)
Determines whether SYMBOL is bound to a QV-PROTO type object."
(cond ((boundp symbol)
(cond ((objectp (eval symbol))
(cond ((eq (send (eval symbol)
:slot-value 'proto-name) 'qv-proto) t)
(t nil)))
(t nil)))
(t nil)))
(defun qv-objects (&optional (prefix-string nil))
"Args: (&optional (prefix-string nil))
Returns a list of the names of all def'd objects that are instances
of QV-PROTO. If PREFIX-STRING is given, only those names begining
with that string are listed."
(cond ((or (stringp prefix-string) (not prefix-string))
(let* ((all-qv-objects (select (variables)
(which (mapcar #'qv-objectp (variables)))))
)
(cond (prefix-string (select all-qv-objects (which
(mapcar #'(lambda (symbol) (begins-with-p prefix-string
(symbol-name symbol)))
all-qv-objects))))
(t all-qv-objects))))
(t (error "Argument to (QV-OBJECTS) must be a string or NIL."))))
;
; Next define functions which compute quasi-variances from Lisp-Stat
; model objects. Such model objects need to have available a method named
; :COEF-COVARIANCE-MATRIX, which delivers the variance-covariance matrix
; of estimated parameters.
;
; To make these functions work on all linear, nonlinear and generalized
; linear regression models which inherit from REGRESSION-MODEL-PROTO:
(defmeth regression-model-proto :coef-covariance-matrix ()
"Method args: ()
Returns the estimated variance-covariance matrix of the
regression parameter estimates."
(* (send self :xtxinv) (^ (send self :sigma-hat) 2)))
; For GEE or Cox regression model objects, for example, use one
; of the following as appropriate:
;
; (defmeth GEE-proto :coef-covariance-matrix ()
; "Method args: ()
; Returns the estimated variance-covariance matrix of the
; regression parameter estimates."
; (send self :covariance-matrix))
;
; (defmeth cox-regression-proto :coef-covariance-matrix ()
; "Method args: ()
; Returns the estimated variance-covariance matrix of the
; regression parameter estimates."
; (send self :covariance-matrix))
;
; [For more information on GEE and Cox regression model objects, see
; Thomas Lumley's web site at
; http://www.biostat.washington.edu/~thomas/]
;
(defun get-qvs
(model-object subset &key (reference-level 1) (factor-name nil)
(model-name nil))
"Args: (model-object subset &key (reference-level 1) (factor-name nil)
(model-name nil))
MODEL-OBJECT - a model object (which inherits from
REGRESSION-MODEL-PROTO, or is an instance of
GEE-PROTO or COX-REGRESSION-PROTO, for example)
SUBSET - either a predictor-name prefix string or a list of
indices corresponding to the coefficients of the
factor of interest
REFERENCE-LEVEL - the level (if any) whose coefficient has been set to
zero to identify the parameters. This is level 1 by
default. Set this to NIL for a parametrization where
no coefficient has been constrained to zero and the
variance-covariance matrix of all p coefficients is
available, or to another value if a different
coefficient has been set to zero.
FACTOR-NAME - an optional string as a reminder of which factor the
quasi-variance structure refers to.
MODEL-NAME - an optional string as a reminder of which model the
quasi-variance structure refers to.
Computes quasi-variances from the covariance matrix of a subset of
estimated coefficients in a model. The coefficients concerned should
correspond to the levels of a categorical variable (factor).
Result is a quasi-variance object, from prototype QV-PROTO."
(cond
((and reference-level (not (numberp reference-level)))
(error "REFERENCE-LEVEL not a number"))
((< reference-level 1) (error "REFERENCE-LEVEL cannot be 0"))
((and (not (stringp subset)) (not (listp subset)))
(error "SUBSET must be a list or a string"))
((not (and (send model-object :has-method :coef-covariance-matrix)
(send model-object :has-method :predictor-names)
(send model-object :has-method :response-name)
(send model-object :has-method :intercept)))
(error (format nil
"Model object lacks one or more required methods, most likely ~
:coef-covariance-matrix. See the QV Calculator manual.")))
(t (let*
((subset (if (stringp subset)
(let* ((indices
(which (mapcar #'(lambda (name)
(begins-with-p subset name))
(send model-object :predictor-names))))
)
(if (send model-object :intercept)
(+ indices 1) indices))
subset))
(levels (if reference-level (+ 1 (length subset)) (length subset)))
(covmat (effect-covariance-submatrix model-object subset
:reference-level reference-level))
(level-names (effect-level-names model-object subset
:reference-level reference-level))
(qv-holder (make-qv-object covmat))
)
(send qv-holder :factor-name (cond (factor-name factor-name)
((stringp subset) subset)
(t nil)))
(send qv-holder :model-name (cond (model-name model-name)
(t nil)))
(send qv-holder :response-name (send model-object :response-name))
(send qv-holder :level-names level-names)
qv-holder))))
(defun effect-level-names (model-object subset &key (reference-level 1))
"Args: (model-object subset &key (reference-level 1))
Extracts the predictor names corresponding to a specified subset of
the coefficients in a model.
MODEL-OBJECT -- a Lisp-Stat model object
SUBSET -- a list of indices, or a prefix string
REFERENCE-LEVEL -- has value 1 by default. Otherwise set this to a different
value if a different level has its coefficient
constrained at zero, or to NIL if all the levels are
represented in the list of model coefficients."
(let* ((intercept (send model-object :intercept))
(predictor-names (send model-object :predictor-names))
(level-names
(if intercept (select predictor-names (- subset 1))
(select predictor-names subset)))
)
(cond ((eq reference-level 1)
(combine "reference category" level-names))
((eq reference-level levels)
(combine level-names "reference category"))
(reference-level (combine
(select level-names (iseq 0 (- reference-level 2)))
"reference category"
(select level-names
(iseq (- reference-level 1)
(- (length subset) 1)))))
(t level-names))))
(defun effect-covariance-submatrix
(model-object subset &key (reference-level 1))
"Args: (model-object subset &key (reference-level 1))
Extracts the covariance matrix of a subset of the estimated coefficients
in a model. 'Subset' is either a list of coefficient indices
(with origin zero), or a predictor-name prefix string. 'Reference level'
is level 1 by default; this can be set to a different value, or to NIL
for a parametrization where no coefficient has been constrained to zero and
the variance-covariance matrix of all p coefficients is available."
(let* ((proto-name (symbol-name (send model-object
:slot-value 'proto-name)))
(the-whole-covariance-matrix
(send model-object :coef-covariance-matrix))
(submat (select the-whole-covariance-matrix subset subset))
(size (length subset))
)
(cond ((not reference-level) submat)
(t (flet ((up-iseq (i j)
(if (> i j) nil (iseq i j)))
)
(let* ((levels (remove-if-not #'(lambda (x) x)
(combine
(up-iseq 1
(- reference-level 1))
0
(up-iseq reference-level
size))))
)
(select
(bind-rows (repeat 0 (+ 1 size))
(bind-columns (repeat 0 size) submat))
levels levels)))))))
(defun begins-with-p (prefix-string target-string)
"Args: (prefix-string target-string)
Tests whether TARGET-STRING begins with PREFIX-STRING"
; called by #'GET-QVS
(let* ((l (length prefix-string))
(reduced-target-string
(cond ((eq l 0) "")
((> l (length target-string))
nil)
(t (select target-string (iseq 0 (- l 1))))))
)
(equalp prefix-string reduced-target-string)))
;
; Now define the basic driver function for the QV calculator for use
; with a user-supplied covariance matrix, either via the web page or
; directly in Xlisp-Stat
;
(defun qvcalc (levels style data &optional (web nil))
"Args: (levels style data &optional (web nil))
LEVELS - the number of levels of the factor of interest (including
any reference level)
STYLE - the way in which the variance-covariance information will
be supplied. This must be one of
'VC -- the whole variance-covariance matrix, either as a
list or as a Lisp-Stat matrix representation
'VCL -- the lower triangle of the variance-covariance matrix,
as a list
'CORR -- coefficient standard errors (including a zero for any
reference level) and lower triangle of the
correlation matrix (excluding the diagonal), as a
single list.
(Only styles 'VCL and 'CORR are available from the web input
form.)
DATA - a list of the variance-covariance information in the style
chosen
WEB - indicates whether this function is being called from the web
page. If so, the function returns a string summarizing the
quasi-variance approximation, rather than a quasi-variance
object.
Example of use:
(qvcalc 5 'VCL (list
0
0 0.0533
0 0.0428 0.1830
0 0.0390 0.0384 0.1427
0 0.0404 0.0411 0.0387 0.0940))
for the ship-type example given on the QV Calculator web page
(http://www.stats.ox.ac.uk/~firth/qvcalc/); see also file 'example.lsp'.
The result is a quasi-variance object, unless WEB is non-nil in which
case it is a nicely-formatted summary string. The summary string can be
obtained by sending the :SUMMARY message to the quasi-variance object."
(cond (web
(let* ((error (web-error-check levels data))
)
(if (not error)
(let* ((cvmatrix (make-cvmatrix levels style data))
)
(send (make-qv-object cvmatrix) :summary t))
error)))
(t (make-qv-object (make-cvmatrix levels style data)))))
;
; The next function is used only when #'QVCALC is called with WEB non-nil
;
(defun web-error-check (levels data)
(cond
((not levels)
(format nil "Input error: Number of levels not specified"))
((not (integerp levels))
(format nil "Input error: Number of levels should be a ~
single integer"))
((eq levels 1)
(format nil "Input error: A one-level variable does not ~
make sense"))
((eq levels 2)
(format nil "Input error: Quasi-variances are unnecessary ~
for this case. (Think about it: with just 2 levels ~
there's only one contrast and you already have its ~
standard error.)"))
((> levels 20)
(format nil "Input error: That many levels? You must ~
be joking! We're offering free computation here, ~
but resources are not unlimited."))
((not (every #'numberp data))
(format nil "Input error: The data are not all numbers"))
((not (eq (* 2 (length data)) (* levels (+ 1 levels))))
(format nil "Input error: Data length does not match ~
number of levels: if the number of levels is p, ~
you need to supply p*(p+1)/2 numbers for the ~
variance-covariance matrix."))
(t nil)))
;
; The next few functions are used by #'QVCALC in computing quasi-variances
;
(defun make-cvmatrix (levels style data)
; converts the data as input into a the full covariance matrix for the effect
; of interest
(cond ((equalp style 'VC)
(cond ((matrixp data) data)
((listp data) (make-array (list levels levels)
:initial-contents data))))
((equalp style 'VCL) (make-cvmatrix-from-triangle levels data))
((equalp style 'CORR) (make-cvmatrix-from-corrs levels data))))
(defun make-cvmatrix-from-triangle (levels data)
; makes a levels by levels symmetric matrix from data which
; is a list of the contents of the lower triangle
(do* ((i 1 (+ 1 i))
(contents nil (append contents (select data (iseq (- i 1)))
(repeat 0 (- levels (- i 1)))))
(data data (cond ((eq (- i 1) levels) nil) (t
(select data (iseq (- i 1) (- (length data) 1))))))
)
((eq (- i 1) levels) (triangle-to-symmat levels contents))))
(defun make-cvmatrix-from-corrs (levels data)
; makes a symmetric var-cov matrix from a list of standard errors
; followed by correlations in lower triangle form
(let* ((sterrs (select data (iseq levels)))
(corrs (select data (iseq levels (- (length data) 1))))
)
(do* ((i 1 (+ 1 i))
(contents (append (list 1) (repeat 0 (- levels i)))
(append contents
(select corrs (iseq (- i 1)))
(list 1)
(repeat 0 (- levels i))))
(corrs corrs
(cond ((eq i levels) nil)
(t (select corrs
(iseq (- i 1) (- (length corrs) 1))))))
)
((eq i levels)
(* (triangle-to-symmat levels contents)
(outer-product sterrs sterrs))))))
(defun triangle-to-symmat (levels contents)
; called by #'MAKE-CVMATRIX-FROM-TRIANGLE and #'MAKE-CVMATRIX-FROM-CORRS
; converts (list of) matrix contents with zeros above diagonal into
; the symmetric matrix whose contents are the lower triangle
(let* ((matrix (make-array (list levels levels) :initial-contents contents))
(transpose (transpose matrix))
(matrix (+ matrix transpose))
(j (make-array (list levels levels)
:initial-contents (repeat 1 (length contents))))
(i (diagonal (repeat 1 levels)))
)
(* matrix (- j (/ i 2)))))
(defun make-qv-object (cvmatrix)
; computes the quasi-variances and returns a QV-PROTO type object
(let* ((levels (first (array-dimensions cvmatrix)))
(ij-pairs (simple-contrasts levels))
(contrast-variances
(combine (mapcar #'(lambda (ij-pair)
(simple-contrast-variance ij-pair cvmatrix))
ij-pairs)))
(qv-model (normal-exp-model (make-qv-design levels)
(log contrast-variances)
:intercept nil :print nil :verbose nil))
(qv (send qv-model :coef-estimates))
(fitvars (send qv-model :fit-values))
)
(send qv-proto :new
:levels levels
:quasi-variances qv
:covariance-matrix cvmatrix
:simple-contrast-variances contrast-variances
:approximated-simple-contrast-variances fitvars
:simple-contrast-labels ij-pairs)))
(defun simple-contrast-variance (ij-pair cvmatrix)
; computes the variance of a simple contrast
(let* ((i (first ij-pair))
(j (second ij-pair))
(vi (select cvmatrix (- i 1) (- i 1)))
(vj (select cvmatrix (- j 1) (- j 1)))
(cov (select cvmatrix (- i 1) (- j 1)))
)
(list (- (+ vi vj) (* 2 cov)))))
(defun ravel (array)
; unravels the elements of an array into a list, in row-major order
; (as "," in APL)
; called by #'SIMPLE-CONTRASTS
(let* ((length (array-total-size array))
(vector (make-array length :displaced-to array))
)
(coerce vector 'list)))
(defun simple-contrasts (levels)
; returns a list of simple-contrast pairs
(let* ((i (iseq 1 levels))
(ij-pairs (ravel (outer-product i i #'list)))
)
(select ij-pairs (which
(mapcar #'(lambda (pair) (> (second pair) (first pair)))
ij-pairs)))))
(defun qv-design-row (ij-pair k)
; makes a row of the qvcalc design matrix
; IJ-PAIR is a list of 2 integers between 1 and k
(if-else
(mapcar
#'(lambda (x) (or (= x (first ij-pair)) (= x (second ij-pair))))
(iseq 1 k))
1 0))
(defun make-qv-design (k)
; puts together rows to make the QVcalc design matrix
(let* ((i (iseq 1 k))
(ij-pairs (simple-contrasts k))
(rows (mapcar #'(lambda (ij-pair) (qv-design-row ij-pair k))
ij-pairs))
)
(reduce #'bind-rows rows)))
;
; Functions used to format the printed summary of a QV-PROTO object
;
(defun format-cvmatrix (matrix)
; pretty-print the covariance matrix
(let* ((string (with-output-to-string (s)
(print-matrix matrix s :float-digits 2)))
(string (remove #\) string))
(string (remove #\( string))
(string (select string (iseq 4 (- (length string) 6))))
)
string))
(defun format-output (model-name response-name factor-name
matrix levels level-names qv
max-negerr-simple max-poserr-simple
max-negerr-overall max-poserr-overall)
; pretty-print the quasi-variances and related quantities
(setf level-names (if level-names level-names (repeat "" levels)))
(concatenate 'string
(format nil "~%QUASI-VARIANCE SUMMARY ~%")
(if model-name (format nil "~%Model name: ~a" model-name))
(if response-name (format nil "~%Response variable: ~a" response-name))
(if factor-name (format nil "~%Factor name: ~a" factor-name))
(format nil "~%Covariance matrix: ~%~%")
(format-cvmatrix matrix)
(format nil "~%~%")
(format nil
"Quasi-variances (and corresponding quasi-se's) computed from~%~
the above matrix by the QV Calculator: ~%~%")
(cond ((every #'(lambda (x) (not (< x 0))) qv)
(let* ((qse (sqrt qv))
(full-header-string (format nil "~a~%~a~%"
" Level Quasi-variance Quasi-s.e."
" ----- -------------- ----------"))
)
(do* ((i 0 (+ i 1))
(string full-header-string
(concatenate 'string string
(format nil "~10,d ~18,3g ~16,3g ~a~%"
i (select qv (- i 1)) (select qse (- i 1))
(concatenate 'string " "
(select level-names (- i 1))))))
)
((= i levels) string))))
(t (let* ((part-header-string (format nil "~a~%~a~%"
" Level Quasi-variance"
" ----- --------------"))
)
(do* ((i 0 (+ i 1))
(string part-header-string
(concatenate 'string string
(format nil "~10,d ~18,3g ~a~%"
i (select qv (- i 1))
(concatenate 'string " "
(select level-names (- i 1))))))
)
((= i levels) string)))))
(format nil "~%")
(format nil "~a~%~a~%~a~3,1f% and ~3,1f%."
"For the standard error of a simple contrast, i.e., of a difference"
"between two levels, the error incurred by the quasi-variance"
"approximation is between " max-negerr-simple max-poserr-simple)
(format nil "~%~%")
(format nil "~a~%~a~3,1f% and ~3,1f%."
"In the set of *all* possible contrasts the approximation error"
"is between "
max-negerr-overall max-poserr-overall)))