Nonparametric Kernel Smoothing Methods. The sm library in Xlisp-Stat .

In this paper we describe the Xlisp-Stat version of the sm library, a software for applying nonparametric kernel smoothing methods. The original version of the sm library was written by Bowman and Azzalini in S-Plus , and it is documented in their book Applied Smoothing Techniques for Data Analysis (1997). This is also the main reference for a complete description of the statistical methods implemented. The sm library provides kernel smoothing methods for obtaining nonparametric estimates of density functions and regression curves for different data structures. Smoothing techniques may be employed as a descriptive graphical tool for exploratory data analysis. Furthermore, they can also serve for inferential purposes as, for instance, when a nonparametric estimate is used for checking a proposed parametric model. The Xlisp-Stat version includes some extensions to the original sm library, mainly in the area of local likelihood estimation for generalized linear models. The Xlisp-Stat version of the sm library has been written following an object-oriented approach. This should allow experienced Xlisp-Stat users to implement easily their own methods and new research ideas into the built-in prototypes. plug-in method of selecting in nonparametric


Installing and loading the sm library into Xlisp-Stat
Assuming that a version of Xlisp-Stat or Arc is already installed on your computer, the following steps are needed for installing the sm library: 1. get the sm.zip file; 2. for inflating the archive in Windows OS you may use Winzip, while in Unix/Linux you need to type > unzip -a sm.zip Stuffit Expander may be used in Macintosh OS.
3. copy the files extracted from the zipped archive in a directory you have created, for example one called sm; 4. copy the file sm-library.lsp in the directory where you run Xlisp-Stat or Arc. In the latter case, you may want to copy the file in the Extras sub-directory so it will be automatically loaded at the start-up; 5. if needed modify the *sm-home-directory* at the beginning of the file sm-library.lsp. This variable should contain the directory to be added to your search paths. By default, *sm-homedirectory* is set to the directory sm relative to your working directory.
For example, suppose you are working on Unix/Linux system, and you run Xlisp-Stat or Arc from the directory xls. Then, by default the variable *sm-home-directory* is set to "sm/". However, if the source sm lisp files are located in a different directory, say " user/xls/library/sm", then you must edit the file sm-library.lsp as follows: (defparameter *sm-home-directory* #+macintosh "Extras:sm:" #+unix "library/sm/" #+msdos "Extras\\sm\\") So, we simply changed the second row from "sm/" to "library/sm/". For other systems the procedure is analogous, provided that the right row is modified depending on your system.
Once you have installed the sm library, you may simply run Xlisp-Stat or Arc and load the library as follows: > (load "sm-library") If you are using Arc the sm library is automatically loaded at start-up if the file sm-library.lsp is located on the: Windows: Extras sub-directory of your Arc directory; Macintosh: Extras sub-folder of your Arc folder; Unix/Linux: Extras sub-directory of the directory from which you start Arc.
Computational efficiency may be speed-up compiling the library. The function > (sm-library-compile) will do the job for you, and the resulting .fsl files will be saved on the directory specified by *sm-homedirectory*.

Getting on-line help
A very basic on-line help system can be invoked through the command > (sm-help) This function opens a window like that on Figure 1. From this window you may select an item and click on the Print button (or double-clicking on the item itself) to see the corresponding documentation on the listener. The standard lisp procedure can also be used to obtain the documentation associated with a function. For example, suppose you want to display the help for the hnorm function. You need to type the following in the text windows command line following the prompt:  are the joint smoothing parameters, and the kernel function is obtained as the product of two univariate kernels.

The sm-density function
The function sm-density allows nonparametric density estimation in one or two dimensions: sm-density Args: (x &key h (hmult 1) (h-weights nil) (weights nil) (binning nil) (nbins nil) eval-points positive delta (model "none") (display "estimate") if 1D or (display "persp") if 2D (add nil) (props (list 75 50 25)) xlab ylab xlim ylim yht (ngrid *sm-ngrid*) (points-on-plot T) (color 'black) (symbol 'cross) if 1D or (symbol 'disk) if 2D (type 'solid) (width 1)) The only required argument is the data x, which can be a list in the one-dimensional case, or a list of lists or a matrix with two columns in the two-dimensional case. The keyword argument :h provides the smoothing parameter (or two smoothing parameters in the two-dimensional case, one for each dimension), and if this parameter is omitted a normal optimal smoothing parameter is used (in the two-dimensional case the default smoothing parameter is computed separately for each dimension).
In large datasets a procedure called binning can be used for increasing the computational speed. Setting the keyword argument :binning to T, the kernel density estimate is computed on the binned data with number of bins provided by :nbins, or if not provided it is computed by the nbins function. For further details on binning see Section 4.4.
Variable bandwidths may be used throught the keyword argument :h-weights, a list of weights which multiply the smoothing parameter used in the kernel function at each observation (for details on the methodology see Bowman and Azzalini, 1997, p. 17, and the example on script 1.11).
A list of weights may be assigned to each observation using the keyword argument :weights. Weights are used, for example, in the computation of density estimate from binned data.
The points at which the density estimator should be evaluated can be set using the keyword argument :eval-points. If not provided, a regular grid is used for computing the estimates with number of points given by :ngrid. By default :ngrid provided by the global variable *sm-ngrid*: for onedimensional density estimate this gives the number of points in the regular grid used, whereas for twodimensional data, approximately *sm-ngrid*/2 number of points for each axis are used.
Data with positive bounded support can be handled setting :positive to T. In this case, a log transformation is applied to the data before construction of a density estimate. The result is then transformed back to the original scale (the method is described in Bowman and Azzalini, 1997, pp. 14-16).
The graphical output is controlled by the keyword argument :display. The setting "none" or NIL will prevent any graphical output from being produced. In one dimension, the default setting "estimate" will produce the density estimate, while the setting "se" will in addition produce a variability band, showing the variability, but not the bias, of the estimate. In two dimensions, the default setting "persp" will produce a perspective plot of the density estimate, while the setting "slice" will produce a slice or contour plot. Observed points are added to the plot if :points-on-plot is set to T, unless a perspective plot is drawn.
In one dimension a reference band, indicating where a density estimate is likely to lie when the data are normally distributed, may be superimposed on the plot setting the keyword argument :model to "Normal". The default "none" assumes no reference model.
By default any plot is drawn on a new window, but a plot may be added to an existing one if you provide a graph object to the keyword argument :add.
The function sm-density returns a sm-density-proto object as its result. Hence, the user may assign the returned object to a variable and then send messages to it.
Example 2.1 (aircraft span data) Consider the data on aircraft technology described in Section 1.2 of Bowman and Azzalini (1997). This dataset can be analyzed in Arc loading the file sm-aircraft.lsp. Suppose we want to obtain a density estimate for the aircraft span on a log scale for the period 1956-1984. The following expressions must be typed in the listener window: > (load "sm-aircraft") > (def y3 (log (select span (which (= period 3))))) > (def sm (sm-density y3 :xlab "Log Span")) After loading the data, a variable containing the log-span for the third period is defined. Then, the last line produces a graph of the density estimate and returns a sm-density-proto object which is assigned to the variable sm (but you may want to use a different name as well). Slot values can be retrieved from this object by sending the messages summarized in Table 1. The results from the density estimation are stored in the slots shown in Table 2.
Information from the resultant object can be obtained by sending to it the :help message:  ADD ADD-METHOD ADD-SLOT BAND COLOR COMPUTE DELETE-DOCUMENTATION DELETE-METHOD  DELETE-SLOT DELTA DENSITY-EVAL-1D DENSITY-EVAL-2D DENSITY-POSITIVE-1D  DENSITY-POSITIVE-2D DIAGNOSTIC-PLOTS DISPLAY DOC-TOPICS DOCUMENTATION EST-CDF  ESTIMATE EVAL-POINTS GET-METHOD GRAPH H H-WEIGHTS HAS-METHOD HAS-SLOT HELP  HMULT INTERNAL-DOC ISNEW LOWER  Help on a specific topic can be obtained as follows: > (send sm :help :plot) PLOT Args: () Plots the density estimate unless the display slot is set to "none" or NIL. if a band for normality is computed (T) or not (NIL) :display controls the graphical output :points-on-plot logical flag which controls whether or not observed points are drawn :props list defining the proportions of the data to be included within each contour in a slice plot for two-dimensional data :add object address which contains the graphical output if required, NIL otherwise As an example of their use, suppose we want to add to the previous plot a density estimate for the other two periods. This is very easy, since we need only to get the graph address from the sm-density object created, and then call the function sm-density specifying that we want to add the density estimate to an existing graph.
Example 2.2 (tephra data) As a further example, script 2.6 contained in the file sm-script.lsp shows how to obtain reference bands for normality using the tephra data. These data record the percentages of aluminium oxide found in samples from a tephra layer resulting from a volcanic eruption in Iceland around 3500 years ago. Since the variable Al203 is a percentage, before computing the density estimate a logit transformation is applied.
In the previous estimate we did not provide a bandwidth, so by default a normal optimal bandwidth has been used (see 4.2, and Bowman and Azzalini, 1997, pp. 31-32). It is known that bandwidths chosen using this criterion tend to oversmooth if the underlying density is not normal. Perhaps a better approach would be to select the bandwidth based on some other more general criterion, as for example by cross-validation (see Section 4.2, and Bowman and Azzalini, 1997, pp. 32-34). The following instrunctions allow to use the cross-validation criterion to select the bandwidth, and then compute the density estimate.
> (hnorm logit) 0.0263566 > (setf hcv (hcv logit)) 0.017678 > (sm-density logit :h hcv :model "Normal" :xlab "logit(Al203)") The hnorm function returns the normal optimal smoothing parameter, which is larger than the bandwidth choosen by cross-validation using the hcv function. The density estimate based on the latter bandwidth is shown on the right panel of Figure 3. Bowman and Azzalini (1997, p. 40, script 2.6) uses the Sheather-Jones plug-in criterion to select the bandwidth, and this provides a value of 6 very close to that selected by cross-validation.
Figure 3: Density estimates of the tephra.lsp data with reference bands for normality. On the left panel a normal optimal smoothing parameter has been used, whereas a smoothing parameter selected by cross-validation on the right panel. Table 3: Methods for computing and graphing density estimates Method Description :compute starts computing density estimates based on the information stored in the slots of the current object :plot draws a plot of the current density estimates Table 3 may be of interest to experienced Xlisp-Stat users. The two methods shown are responsible for computing and plotting the density estimates. This may be helpful if we want to re-compute or re-graph the density estimate based on the same data but changing one or more parameters. For instance, in some circumstances we would like to investigate the sensitivity of our nonparametric estimate to bandwidth size changing. Thus, we may define different bandwidths and plotting the results on the same graph to allow comparisons. Two approaches may be adopted. The first approach defines the list of values for the bandwidth, computes the first density estimate, and then adds the successive estimates to the existing plot. Taking a different approach, we may define only one sm-density-proto object, and then repeatedly modify the parameter hmult to control the amount of smoothing applied. > (def sm (sm-density logit :xlab "logit(Al203)")) > (send sm :h) get the default bandwidth used 0.0263566 > (send sm :hmult) get the default multiplier of h 1 > (def hmult (list 0.5 0.75 1.25 1.5)) define a list of values for hmult > (send sm :add (send sm :graph)) set the plot for adding the density estimates > (dotimes (i (length hmult)) starts the loop (send sm :hmult (select hmult i)) ... re-defines the value of the slot hmult (send sm :compute) ... re-computes the estimates (send sm :plot)) ... graphs the results

Example 2.3 (aircraft span data)
Data on the aircraft span based on the first two principal components can be used as an example of density estimation in two dimensions (see also script 1.4). After loading the dataset, a perspective plot and a contour plot using the default normal bandwidth can be obtained as follows: > (load "sm-aircraft-pc") > (def pc3 (list (select Comp1 (which (= period 3))) (select Comp2 (which (= period 3))))) > (sm-density pc3 :xlab "Comp1" :ylab "Comp2") > (sm-density pc3 :display "slice" :xlab "Comp1" :ylab "Comp2")  This function allows a set of univariate density estimates to be compared, both graphically and formally in a bootstrap hypothesis test of equality (see Bowman and Azzalini, 1997, Section 6.2). The data x are grouped based on the argument group, and the same smoothing parameter h is used for each group. If :model is "none" comparison is restricted to plotting only. The alternative value "equal" produces a permutation hypothesis test of equality. For two groups, an appropriate reference band is displayed, unless :band is set to NIL. The number of permutation simulations is controlled by the argument :nsim.

Nonparametric density estimation of stationary time series data
sm-ts-pdf Args: (x &key h lags (display "persp") (ngrid *sm-ngrid* set-ngrid) (points-on-plot T) (varname "x")) This function estimates the density function of a time series x, assumed to be stationary (see Bowman and Azzalini, 1997, Section 7.2). The univariate marginal density is estimated in all cases. Bivariate densities of pairs of lagged values are estimated depending on the keyword argument :lags, a list of values for the lags to consider in the joint distribution estimation. The bandwidth can be set by the argument :h, and if missing a normal optimal bandwidth is used. The :display argument controls the bivariate density plot, with possible values "persp" (default) and "slice".

Nonparametric regression estimation
A general formulation for a nonparametric regression model takes the form denotes the response variable, ¢ the explanatory variable and ¦ an independent error term with zero mean and variance equal to¨H . The aim of a nonparametric regression is to provide an estimate of the smooth function § 8 . A simple kernel estimator of § 8 is the local mean estimator (see Bowman and Azzalini, 1997, p. 49) where the kernel function is a smooth positive function, wich gives weights that decrease monotonically as D i ncreases in size. As for density estimation, a normal kernel function centered at the evaluation point with mean 0 and standard deviation 6 is used. Therefore, 6 acts as smoothing parameter.  Bowman and Azzalini, 1997, p. 50). Extension to two dimensions is straightforward (see Bowman and Azzalini, 1997, p. 52-53).

The sm-regression function
The function sm-regression can be used for nonparametric regression estimation in one or two dimensions: The required arguments are the covariate(s) x (list of values in the one covariate case, or a list of lists or a matrix with two columns in the two covariates case), the response variable y (a list of values), and the smoothing parameter h (a list of one or two values, for the one-and two-dimensional case respectively).
The keyword argument :poly-index controls whether local mean estimator (0) or a local linear estimator (1) is used.
The graphical output is controlled by the keyword argument :display. The setting "none" or NIL will prevent any graphical output from being produced. With one covariate, the default setting "lines" will produce the regression estimate, while the setting "se" will in addition produce a variability band, showing the variability, but not the bias, of the estimate (see Bowman and Azzalini, 1997, pp. 75-77). With two covariates, the default setting "persp" will produce a perspective plot of the regression estimate, while the setting "contour" will produce a contour plot with levels specified by the argument :levels. This must be list of heights of contour lines, and by default 5 levels are used to drawn contours which cover the range of the smooth estimate of y.
The keyword argument :model defines the reference model. The values "none", "no effect" and "linear" are possible. When a model is provided, a reference band for the reference model is plotted on the graph. Moreover, if :test is T a formal test is produced using the reference model as the null hypothesis (see Bowman and Azzalini, 1997, Sections 5.2, 5.3).
Most of the remaining arguments have the same meaning we have already discussed for the smdensity function. Details can be obtained throught the on-line help.
The function sm-regression returns a sm-regression-proto object as its result. Then, the user may assign the returned object to a variable and send messages to it.
The slot values can be retrieved sending the messages summarized on Table 4, and the results from the nonparametric regression are stored in the slots reported on Table 5. As we saw for the density estimation case, slot values can be modified and the needed computations invoked sending the message :compute, while a new plot may be obtained through the :plot method (see Table 3.1). weights associated to each observation :poly-index the estimator used, 0 for local mean estimator and 1 for local linear estimator :model reference model :band when a reference model is provided, this slot is set to T and a band will be plotted on the graph (only for one covariate); setting this slot to NIL will prevent any band to be plotted :display controls the graphical output :points-on-plot logical flag which controls whether or not observed points are drawn :add object address which contains the graphical output if required, NIL otherwise Example 3.1 (brain weight data) The following example considers the data on brain weight (BrainWt) in grams and body weight BodyWt in kilograms for sixty-two species of mammals. This dataset is described by Cook and Weisberg (1999, Section 5.1) and can be analyzed in Arc by loading the file brains.lsp. A linear regression model for BrainWt on BodyWt, both expressed on a log scale, with constant variance function seems reasonable. We may further investigate this linear relationship comparing the estimated parametric model with a nonparametric regression model (see Bowman and Azzalini, 1997, Section 5.3). The cross-validation criterion may be used to select the smoothing parameter, then a plot of the estimated smooth function together with reference bands for linearity allows to perform a graphical check. This can be obtained typing the following commands:  address of the last graph drawn from the object :regression-test computes p-value for a hypothesis test that checks a reference parametric regression model with the estimated nonparametric model. The function hcv returns a list giving the value of 6 which minimizes the cross-validatory criterion and a plot showing the criterion function plotted over the search grid of smoothing parameters (see top-left panel of Figure 5).
The nonparametric regression analysis is then performed through the function sm-regression, which returns an object assigned to the variable sm (again, you may want to use a different name as well). The graphical output shown in the bottom panel of Figure 5 is also produced. From this plot we can see that the estimated smooth function lies inside the reference band over almost all the range of the predictor. A formal hypothesis test suggest that there is no strong evidence against the reference linear model. Since this test depends on the smoothing parameter selected, it is possible to obtain the observed significance for a range of smoothing parameters as follows: The significance trace plot, shown on the top-right panel of Figure 5, suggests no strong evidence against the null hypothesis (for further details on the methodology see Bowman and Azzalini, 1997, p. 93). equivalence ratio (E), a measure of the richness of the air and fuel mix in the engine, and of compression ratio of the engine (C). This dataset has been studied extensively, among others, by Cleveland (1993). We start by loading the data and using the function hcp for selecting the bandwidths (see Section 5.4), both separately and jointly for the two predictors. We also provide the corresponding approximate span computed as defined in Section 4.1. The selected bandwidth for E is equal to 0.0205786, which corresponds to a span of approximately 0.22. The same function is then applied to C, but the criterion does not converge even after re-adjusting the seach grid. The problem is due to the fact that NOx is almost constant as C varies, so the criterion tends to select the largest bandwidth value available in the search grid. This difficulty also emerges if we use both predictors: the bandwidths selected are extremely large. Adopting an heuristic approach for bandwidth selection, we computed the values corresponding approximately to a span of 0.5, and then we used them for smoothing data in two dimensions.
:xlab "E" :zlab "C" :ylab "NOx") > (def sm (sm-regression (list E C) NOx :h (list 0.058 0.875) :display "contour" :levels '(3 2 1 0) :xlab "E" :zlab "C" :ylab "NOx")) > ( The first fitting produces the perspective plot shown in the left panel of Figure 6, while the latter gives the contour plot in the right panel of Figure 6, where levels are computed at selected value of NOx estimated. From such graphs the underlying relationship between the response and the two predictors can be easily assessed.

Nonparametric analysis of covariance
A nonparametric analysis of covariance model (see Bowman and Azzalini, 1997, Sections 6.4, 6.5) may be written as observations within each group. Furthermore, each group has its own smoothing function . This general nonparametric analysis of covariance model may be adapted to reflect the assumption of parallel nonparametric regression functions, that is we may write where there is a single smooth function § , but different intercepts © ' for each group. The sm library can be used to perform a nonparametric analysis of covariance using the function sm-ancova Args: (x y &key group h (display "lines") (model "none") (band T) (test T) (h.alpha nil) weights (binning nil) (nbins nil) (add nil) (ngrid *sm-ngrid*) eval-points xlab ylab (points-on-plot T) color symbol type width) The required arguments are the covariate x, the response variable y, the grouping variable group, and the smoothing parameter h.
The reference model can be set using the keyword argument :model, which can take the values "none", "equal" and "parallel". Setting :band to T a reference band for the reference model will be produced (only in the case of two groups). A formal test can also be obtained setting :test to T.

Nonparametric estimation of the autoregression function
The dependence structure of a process ¡ ¡ G can be studied nonparametrically analyzing the conditional mean function or autoregression function. This may be written as (see Bowman and Azzalini, 1997, Section 7 As usual, the argument x provide the data, in this case a list of time series values, while h the bandwidth used for kernel smoothing. The keyword argument :d provides the number of past observations used for conditioning (1, the default, or 2), whereas :lags needs to be a list containing the lags considered for conditioning. If :se is set to T pointwise confidence bands are computed of approximate ¡ level.

Nonparametric analysis of repeated measurements data
A repeated measurements model can be formulated in a nonparametric context as follows (see Bowman and Azzalini, 1997, Section 7 is a function of time, and the errors ¦ ' are correlated within each individual and independent between individuals (a further assumption often made is that of stationarity).
Such a model can be fitted using the sm library through the function sm-rm Args: (y &key Time This function estimates nonparametrically the mean profile for repeated measurements (i.e. longitudinal data) from a matrix y, with each rows associated to individuals and columns associated to observation times. The smoothing parameter is computed using the modified Rice criterion, which selects the bandwidth taking into account the estimated autocorrelation function (see Bowman and Azzalini, 1997, p. 139).
Further details are provided in the on-line help. Scripts 7.3 and 7.4 show some statistical analyses on repeated measurements data.

Nonparametric regression with autocorrelated errors
A nonparametric regression model with autocorrelated errors can be stated as follows is a smooth function which depends on single covariate observed over time ( ! , and the errors ¦ are generated by a stationary stochastic process with mean 0 (see Bowman and Azzalini, 1997, Section 7.5).
The function in the sm library which fits a nonparametric regression model with autocorrelated errors is the following sm-regression-autocor Args: (y hfirst &key (x nil) (minh 0.5) (maxh 10) (hngrid 50) This function estimates nonparametrically the regression function of y on x when the error terms are serially correlated. The keyword :method specifies the optimality criterium adopted: possible values are "no.cor" "direct" (default), and "indirect".

Other functions in the sm library
The file sm-functions.lsp contains the definition of global variables used by the sm library (see Table 7), and several functions used in conjunction with methods for nonparametric density and regression estimation.

Controlling the span and the bandwidth of the smoothing parameter
The amount of smoothing applied to the data is controlled by the smoothing parameter 6 , which in the sm library corresponds to the standard deviation of a gaussian kernel function with zero mean. Often, since the bandwidth is expressed in absolute value, a more easily interpretable measure is computed, namely the span, which is the proportion of sample data providing positive weights. Since we are using a normal kernel with standard deviation 6 , we may reasonably assume that data points outside a range of ¡ £ ¢ times the standard deviation receive virtually zero weight. Hence, there exists the following approximate relations Default value Description *sm-band* T If T reference bands produced by sm-density and sm-regression are plotted as a shaded area, otherwise as lines. *sm-ngrid* 100 The number of points in a regular grid used to plot the kernel estimate. For two-dimensional data approximately, *sm-ngrid*/2 number of points are used along each dimension. *sm-glm-ngrid* 25 The number of points in a regular grid used for fitting local likelihood models. *sm-glm-aic* NIL If T the corrected version of AIC is used in local likelihood model (see Section 6.8) *sm-plot-size* '(600 450) Sets the size of the graphs produced through the sm library.
Note that the above approximations turns out to be exact, except near the boundaries, when data points along the x-axis are evenly spaced. Hence, the approximate value of the span is really the proportion of the observed range of ¢ which receives positive weights for a given evaluation point. In this respect, it provides a simple descriptive statistic of the amount of smoothing applied to the data.
Two functions are available for computing the span for a given bandwidth and viceversa.

sm-bandwidth
Args: (x span &key (nnbr nil)) This function calculates the bandwidth for the data x (a list, a list of lists or two-columns matrix) corresponding to the span provided by the argument span. If nnbr is T a list of weights based on nearest neighbour distances is calculated; this allows to use variable bandwidths. In two dimensions, the same calculations are applied to each dimension. If nnbr is NIL only the bandwidth is returned, otherwise a list with the bandwidth and a list of weights. These weights are computed as is the th nearest neighbour distance of the covariate value ¢ ' (computed with the function (nnbr x k) for k=round(n*span)) and for an overall bandwidth 6 (see Bowman and Azzalini, 1997, pp. 17-19). Weights based on nearest neighbour distances can be provided as argument to :h-weigths in both the sm-density and sm-regression functions.

Args: (x h)
This function returns the span for the data x (a list, a list of lists or two-column matrix) corresponding to the bandwidth h.
Example 4.1 (follicle data) The follicle.lsp dataset can be used to show how variable bandwidths may be fitted with the sm library. The data concern the number of ovarian follicles counted from sectioned ovaries of women of various ages. The scatterplot of log [Count] vs Age shows that points are sparse along the x-axis, with data denser at older ages. Using a span of 0.5 we compute a nonparametric regression curve with fixed bandwidth and one with variable bandwidth. For comparison purposes also a lowess estimate is computed and then added to the plot. The last command provides a legend. > (load "follicle") > (def h (sm-bandwidth Age 0.5 :nnbr t)) > h

Functions related to density estimation
Nearest neighbour distances from data in one or two dimensions nnbr Args: (x k) This function calculates the nearest neighbour distance from each value in ¢ to the remainder of the data. In two dimensions, Euclidean distance is used after standardizing the data to have unit variance in each component.

Mean integrated squared error for density estimation with normal data
nmise Args: (sd n h) This function evaluates the mean integrated squared error of a density estimate which is constructed from data which follow a normal distribution (see Bowman and Azzalini, 1997, pp. 26-17). This function evaluates the smoothing parameter which is asymptotically optimal for estimating a density function when the underlying distribution is normal (see Bowman and Azzalini, 1997, pp. 31-32).

Cross-validatory choice of smoothing parameter
hcv Args: (x &key (y nil) (h-weights nil) (ngrid 8) (hstart nil) (hend nil) (display "none") (add nil)) This function uses the technique of cross-validation to select a smoothing parameter suitable for constructing a density estimate or nonparametric regression curve in one or two dimensions (see Bowman and Azzalini, 1997, pp. 32-34, 77-79). This function computes a cross-validatory criterion, based on integrated squared error, for use in selecting a smoothing parameter in nonparametric density estimation. This function uses the Sheather-Jones plug-in method for selecting a bandwidth which is suitable for constructing a density estimate in the one-dimensional case (see Bowman and Azzalini, 1997, pp. 34-36). This function computes a criterion associated with the Sheather-Jones plug-in method of selecting a bandwidth in nonparametric density estimation.

Sheather-Jones criterion for nonparametric density estimation
Integrated squared error between a density estimate and a normal density nise Args: (y &optional (hmult 1) (nbins (nbins y))) This function evaluates the integrated squared error between a density estimate constructed from a standardized version of the univariate data y and a standard normal density function (see Bowman and Azzalini, 1997, pp. 38-40). The two functions return, respectively, the weighted mean and the weighted variance of x with weights given by w, where x and w can be a list or a vector.

A significance trace for a hypothesis test comparing a reference parametric model with a nonparametric one
sig-trace Args: (smreg hvec &key (plot T)) This functions creates a significance trace for a hypothesis test based on a nonparametric smoothing procedure. The argument smreg must be a sm-regression object and hvec a list of smoothing parameters for which the test will be performed. It :plot is T the p-value of the test is plotted against a range of smoothing parameters (see Bowman and Azzalini, 1997, pp. 89, 93).

sm-weight
Args: (x h &key (eval-points x) (hmult 1) (h-weights (repeat 1 (length x))) (poly-index 1) (cross nil) (weights (repeat 1 (length x)))) For the one-covariate case, this function computes a length(eval-points) by length(x) smoother matrix such that the local polynomial smooth of y on x with bandwidth h and weights hweights is given by y.
An analogue function called sm-weight2 allows to obtain the smoothing matrix in nonparametric regressions with two covariates.

sm-sigma
Args: (x y &key (diff-ord 2) (devs nil) (weights (repeat 1 (length y)))) This function uses ideas of local differencing to estimate the standard deviation of the errors in a nonparametric regression model with one covariate. Simple first-order differencing of pairs of neighbouring observations, or a method based on pseudo-residuals constructed from three neighbouring observations, may be used (see Bowman and Azzalini, 1997, pp. 73-75).

Binning
When the number of data points is high, the method of binning (see Wand and Jones, 1995, Appendix D) may be used to avoid lengthy computations and use of very large matrices. If binning is used in density estimation, the data are tabulated with the use of the function binning, and the outcome is passed with the frequencies computed as weights.
binning Args: (x &key (y nil) (breaks nil) (nbins nil)) Given the data in x, this function constructs a frequency table associated to appropriate intervals covering the range of x. This can be a list, a vector, or a matrix with one column in the one dimensional case, whereas a list of lists, or a matrix with two columns in the two dimensional case. If you want to use binning in nonparametric regression you must provide a list or a one-column matrix of values for the response variable with :y. The argument :breaks, either a list or a matrix with two columns, assigns the division points of the axis, or the axes in the matrix case. If :breaks is not given, it is computed by dividing the range of x into :nbins intervals for each of the axes. The argument :nbins gives the number of intervals on the x axis (in the univariate case), or a list of two elements with the number of intervals on each axes of x (in the bivariate case). If :nbins is not given, a value is computed as (round (+ (log (length x) 2) 1)) or using a similar expression in the bivariate case. The function returns a list of several elements, and for details see the on-line help.
A sort of "optimal" number of bins is returned by the following function: For the list or list of lists x, given the cutoff value k, returns the number of bins to use according to the following rule: where n is the number of observations and ndim is the dimension of x.

Extra functions for nonparametric density and regression
The file sm-addson contains code not provided in the original sm library and implements methods discussed mainly by Loader (1999).

Diagnostics for density estimation
Once we have obtained a density estimate for the data at hand, we may ask whether such estimate fit the data or not. This question is not easily answered since there is no natural definition of residuals as for regression problems, and no reference model for comparison, such as a saturated or null model. Loader (1999, pp. 87-90) suggests two graphical checks: 1. compare the integral of the density estimate and use the latter (which is the MLE for the true cumulative distribution function) as a reference for comparing the distribution function based on kernel estimation. . Also in this case, any departure from a straight line indicates lack of fit.

use PP-plot and QQ-plot for
The library allows, only for the one-dimensional case, to obtain a plot of the estimated cdf vs the emprical distribution function, and the PP-plot for the estimated cdf.
Example 5.1 (aircraft span data) Recalling the aircraft.lsp dataset, we may get a density estimate for the log span on the third period and these dignostic plots as follows: > (load "sm-aircraft") > (def y3 (log (select span (which (= period 3))))) > (def sm (sm-density y3 :xlab "Log Span")) > (send sm :diagnostic-plots) The method :diagnostic-plots produces the plots shown in Figure 8. The estimated CDF is computed by the :est-cdf method.  , can be retrieved sending the message :residuals. A plot of the residuals vs the explanatory variable (one-predictor case) or the fitted values (two-predictors case) is produced by the method :plot-residuals.
The approximate "number of parameters" used by the smoothing procedure is given by the degrees of freedom of a local fit. In analogy with the definitions used in parametric linear models (see Hastie and Tibshirani, 1990, Sec. 3 These results form the basis for some diagnostic procedures (see Loader, 1999, pp. 27-29). Assume the evaluation points are the actual observed values of the explanatory variables, so is a % © % smoothing matrix. The diagonal elements of , denoted as ¤ ' ¢ ' , may be used as a measure of the sensitivity of the estimated curve to the individual data points. Another useful quantity is the variance reducing factor £ ¢ H which measures the reduction in variance due to the local regression. It can be shown that The two extreme cases correspond, respectively, to a smoothed curve which is equal to the sample average of the response variable, and a curve which interpolates exactly the observed points.
Diagnostic plots can be obtained through the method :diagnostic-plots. As a result a graph with the following functions is drawn: Influence function: a plot of diagonal elements of vs the explanatory variable (one-predictor case) or the fitted values (two-predictors case); Variance function: a plot of variance reducing factors vs the explanatory variable (one-predictor case) or the fitted values (two-predictors case).
Finally, a printed summary for the fitted nonparametric model is provided when the evaluation points are set equal to the observed points. In any case, this summary may be obtained sending the message :summary to a sm-regression object.

Generalized cross-validatory choice of smoothing parameter for nonparametric regression
The generalized cross-validation criterion can be used to select a smoothing parameter for a nonparametric regression curve in one or two dimensions. The GCV statistic is an estimate of the predictive mean squared error (PMSE), and it is given by is equal to the trace of the smoothing matrix, so % 1 ¤ £ is an estimate of the degrees of freedom already discussed.
In the sm library the GCV criterion can be invoked using the following function: hgcv Args: (x y &key (h-weights nil) (ngrid 10) (hstart nil) (hend nil) (nnbr nil) (display "none") (add nil) (color 'black) (type 'solid) ) The x argument provides the explanatory variable (a list, a vector, or a single column matrix) in the onedimensional case, or the explanatory variables (a list of lists or a two-column matrix) in the two-dimensional case. y is a list of response values. If the argument :nnbr is set to T a nearest neighbour bandwidth is used. The argument :display controls the graphical output: any character setting other than "none" will cause the criterion function to be plotted over the search grid of smoothing parameters. The particular value "log" will use a log scale for the grid values, whereas "df" will use the fitted degrees of freedom on the x-axis. By default fitted df are computed using ¤ # in Section 5.2, whereas ¤ H is used if :display is set to "df2". The remaining arguments were already discussed in the hcv function (see Section 4.2).

Mallows' CP criterion for selecting the smoothing parameter for nonparametric regression
Mallows' CP criterion can be used to select a smoothing parameter for nonparametric regression in one or two dimensions. The CP criterion is an estimate of the squared estimation error, and is given by An estimate of¨H is required. Cleveland and Devlin (1988) suggested to use the normalized residual sum of squares evaluated at the smallest bandwidth considered (i.e. where the bias is the smallest): The CP criterion can be applied using the following function: hcp Args: (x y &key (h-weights nil) (ngrid 10) (hstart nil) (hend nil) (nnbr nil) (display "none") (add nil) (color 'black) (type 'solid) ) All the arguments have the same meaning already discussed in the previous section for the hgcv function.
Example 5.2 (motorcycle data) Considers the mcycle.lsp dataset (Silverman, 1985), which records data from an experiment to test crash helmets. A series of measurements of head acceleration after a simulated motorcycle crash were recorded (accel), together with the time in milliseconds since impact (times). Since observations are sparse for large values of the covariate, we might think to use a variable bandwidths based on nearest neighbour distances.
> (load "mcycle.lsp") > (def hgcv (hgcv times accel :nnbr T :ngrid 20 :display T)) (1.53183 #<Object: ca4160, prototype = SCATTERPLOT-PROTO>) > (def hcp (hcp times accel :nnbr T :ngrid 20 :display T)) (1.57266 #<Object: c9dfa0, prototype = SCATTERPLOT-PROTO>) > (sm-span times 1.53) 0.166304 The smoothing parameter selected by the GCV criterion and the CP criterion are almost the same (see also Figure 9), and corresponds to a span of about 0.166. Since we are using variable bandwidths, we first need to get the weights for any observation at the selected span. Then, we fit a local linear regression model and we produce some diagnostic plots.  The top panel of Figure 10 shows both the smooth regression function obtained using variable bandwidths (solid line) and the smooth function obtained with fixed bandwidth as follows:

Local likelihood estimation for generalized linear models
Local likelihood fitting was proposed by Tibshirani (1984) and Hastie and Tibshirani (1987). The basic idea is that whenever a likelihood function is available we can maximize it locally in order to obtain an estimate of a response variable at some point, say ¢ ¡ . Data points in a neighbourhood of ¢ ¡ are weighted according to a smooth function which gives weights that decrease as ¢ ' increases in distance from ¢ ¢ . A technical discussion of the local likelihood approach is given by Fan et al. (1995), while a full account of the topic is provided by Loader (1999).
This approach allows to fit a smooth curve on a scatterplot of data, taking directly into account the nature of the response variable. In the case of a binomial response, and even more for a bernoulli variable, smoothers developed for a continuous variable are not satisfactory because they do not use the information available. For instance, if we smooth a scatterplot of a binary vs a continuous predictor with a smoother for continuous data, a clear drawback is that we are not taking into account the binary nature of the response, so the fitted smooth curve could lie outside the range ¥ ¥ # . On the contrary, using a local likelihood approach we can model directly this information and obtain a smooth curve that fit the data in a more sensible way. The same arguments apply for bounded variables, such as nonnegative variables, or discrete variables, such as those expressing counts.
The sm library provided by Bowman and Azzalini (1997) allows to smooth binomial and Poisson data. The corresponding Xlisp-Stat functions are described in Section 6.10. A more general approach suitable for any GLM family should be included in the further release of the sm library for S-Plus. In the current release of the sm library for Xlisp-Stat, we developped several functions for applying the local likelihood approach to any GLM, and these are described in the following sections.

A review of the theory of local likelihood for GLMs
Assume that a response variable ¦ has a distribution in the exponential family, taking the form where ¥ , the canonical parameter, and @ , the dispersion parameter, are scalar values. § , , are known functions, which determine the specific parametric family of distributions (see McCullagh and Nelder, 1989,  . Maximum likelihood estimates for are obtained by setting the first derivative of the log-likelihood with respect to each parameter equal to zero and then solving the system of equations. Since this is usually a set of nonlinear equations, estimates are obtained numerically through the Fisher scoring algorithm or the Newton-Raphson method (Nelder and Wedderburn, 1972).

Local likelihood function
Suppose we have a single predictor, and we denote the log-likelihood as 4 ¡ . The local likelihood approach essentially applies a weighting mechanism to the log-likelihood, so we can write is a weighting function that determines the relative importance of each contribution Only the observations that lie in this neighbourhood will contribute to the log-likelihood. The proportion of points that lie in it is called the span. Thus, the smoothing parameter controls the degree of smoothing applied to the data.
The local likelihood in equation (1)  , we need to define the weights to use on each where the generic ¤ ¥ ¡ -th element is given by . Therefore, the ¦ matrix has elements equal to the difference between the point of estimation and the observed value of the predictor. The smoothing weights are then obtained by applying a smoothing function, which depends on have approximately zero weight and their contribution to the local likelihood will be nearly zero.
Once we have obtained the weights we can use the IRWLS procedure for fitting a generalized linear model using as prior weights , and repeating the fitting procedure times, one for each evaluation point

Variability bands
As in all statistical models, it is useful to have a measure of the uncertainty associated with the estimated values. Usually this is done providing standard errors for the estimates and/or confidence intervals for a given confidence level. An informal approach can be pursued providing the so-called variability bands (see Bowman and Azzalini, 1997, pp. 75-76)

Approximate degrees of freedom
As we saw in Section 5.2 the approximate fitted degrees of freedom expresses the "number of parameters" used by the smoothing procedure. For generalized linear models, the number of coefficients estimates can be computed through the hat matrix defined as is the deign matrix and ¢ is the matrix of known weights obtained at the end of the iterative fitting algorithm (see McCullagh and Nelder, 1989, p. 405). In analogy with the linear regression models, the trace of the hat matrix equals the number of parameters.
In a local likelihood GLM with evaluation points equal to the observed points, for each the corresponding leverage is given by is the Residual degrees of freedom are easily obtained as ¡ ¢ £ ¢ ¥ ¤ ¦ ¢ ! % 1 ¥ ¥

Goodness of fit measures
A general measure of goodness of fit used in GLM is the deviance. This is defined as twice the difference between the log-likelihood of a saturated model (which has as many parameters as the number of observations) and the log-likelihood for the current model. The exact form depends on the exponential family distribution under study, and it can be expressed as the sum of

Residuals
Residuals for local likelihood are defined similarly to residuals for GLM based on the full likelihood. Some possible choices are: is the variance function evaluated at the ¡ th observation.

Choosing the smoothing parameter
As for linear nonparametric regression models, the choice of the smoothing parameter is a crucial step in local likelihood models. An extension of the Akaike Information Criterion (AIC) can be used to guide this selection. The AIC, which is an estimate of the expected Kullback-Leibler information (Burnham and Anderson, 1998), is defined as The slot values can be retrieved sending the messages summarized on Table 8, and the results from the local likelihood fitting are stored in the slots reported on Table 9. As we saw for density and regression smoothing, slot values can be modified and the needed computations invoked sending the message :compute, while a new plot may be obtained through the :plot method (see Table 6.9). where the arguments have the same meaning as in the hcv, hgcv, and hcp functions for bandwidth selection already discussed. Furthermore, arguments required by the sm-glm function must be provided as well.

Local likelihood for logistic and Poisson regression and PRLT test
In the sm library for S-Plus two functions are available for fitting binomial and Poisson local likelihood models with canonical link function. These functions are also available in the Xlisp-Stat version of the library. Actually, they merely provide a quick way of calling the sm-glm function with the required arguments but without specifying the family and link function to use.

sm-logit
Args: (x y &key (trials (repeat 1 (length y))) (h nil) (ngrid *sm-glm-ngrid*) (eval-points nil) (display "estimate") (verbose t) (add nil) (xlab "x") (ylab "y") (points-on-plot T) (color 'black) (symbol 'disk) (type 'solid) (width 1)) This function estimates the logistic regression curve using the local likelihood approach for the observations on a binomial response and an associated set of covariate values. The arguments have the same meaning discussed for the sm-glm function (see Section 6.9). This function estimates the regression curve using the local likelihood approach for a vector of Poisson observations and an associated vector of covariate values. Also in this case, the arguments have the same meaning discussed for the sm-glm function (see Section 6.9). A comparison between a parametric GLM and a nonparametric model may be performed through the so-called Pseudo-Likelihood Ratio Test (see Azzalini et al. 1989, Azzalini andBowman, 1997, pp. 98-104). This bootstrap goodness-of-fit test allows to test a parametric model vs a nonparametric model. The functions for applying such tests are: sm-logit-bootstrap Args: (x y &key (trials (repeat 1 (length y))) (h nil) (nboot 100) (degree 1) (phi 1) (plot nil) (xlab nil) (ylab nil)) This function performs a Pseudo-Likelihood Ratio Test for the goodness-of-fit of a standard parametric logistic regression. The argument nboot gives the number of bootstrap samples (by default is set equal to 100). The parametric model is fitted using a polynomial in x with degree given by :degree (default 1). A dispersion parameter can be provided through the argument :phi. The function returns a list of lists containing: the plot object representing the bootstrap samples; the observed value of the Pseudo-Likelihood Ratio Test statistic; the observed p-value estimated via the bootstrap method.
A similar function is also provided for a Poisson regression model:

sm-poisson-bootstrap
Args: (x y &key (offset nil) (h nil) (nboot 100) (degree 1) (scaled-dev nil) (plot nil) (xlab nil) (ylab nil)) The arguments take the same meaning as in the previous function, except :scaled-dev which if set to T allows the differences in deviance to be scaled by an estimate of the dispersion parameter. The function returns a list as in the previous case.
Example 6.1 (mortality rate) The dataset from Henderson and Sheppard (1919) in the file mortable.lsp reports the observed number of deaths (D) for ages from 55 to 99 (Age), and the corresponding total number of patients in each class (N) . We could think to model this data as they come from a binomial distribution where the response variable is the number of death people and the trials are the number of subjects in the class. Suppose we want to estimate a function that smooths the data and gives a summary of the mortality rates. To apply a local likelihood fitting, we must choose the smoothing parameter. The AIC criterion applied to this dataset gives the following results: > (load "mortable") > (def h (haic-loclik Age D :trials N :family "binomial" :display T :ngrid 20)) AIC criterion for local likelihood models ---------------------------------------   Figure 11 shows the smooth curve obtained applying the local likelihood estimation with smoothing parameter 6 ! © § . We also plot the variability bands for the estimated smooth curve. The death rate increases as age increases up to age 88 approximately, then becomes constant and rapidly increases after 95 years. However, we must be careful in the interpretation of this latter shape, since there are few observations after age 95, and for age 99 just one case, which is a death. This uncertainty is revealed also by the width of the variability bands, which is quite large for older and younger ages.
Depending on the context we could stop here the analysis. Anyway, we proceed a little further, assuming that, based on the nonparametric local fitting, we want to seek for a suitable parametric model that captures the main features on the data. We fitted a polynomial logistic model and Table 11 shows that a linear and a quadratic term in Age seem to be required. A way of testing the adequacy of a parametric model compared to a nonparametric one obtained through local likelihood may be obtained via the pseudo-LRT:  The observed significance seems to suggest that a second degree polynomial on Age should suffice, so providing a further confirmation of what discussed previously.
Example 6.2 (birds in paramo vegetation) A study was conducted to investigate the number of bird species in isolated islands of paramo vegetation in the northern Andes. The aim was to investigate how the number of species is related to some covariates. The dataset is contained in the file paramo.lsp and we focus on the relationship between the number of species and the area of the island, so a Poisson model seems to be reasonable.
Again, the AIC criterion can be used for bandwidth selection. Since the number of observations is small, we may adopt the corrected version of AIC (see Section 6.8) by setting the global variable *SM-GLM-AIC* to T. Then, we fit the Poisson local likelihood model.  The fitted curve is shown on the right panel of Figure 12.

Nonparametric regression with survival data
The function sm-survival creates a smooth, nonparametric estimate of the quantile of the distribution of survival data as a function of a single covariate. A weighted Kaplan-Meier survivor function is obtained by smoothing across the covariate scale. A small amount of smoothing is then also applied across the survival time scale in order to achieve a smooth estimate of the quantile (for details on the methodology see Bowman and Azzalini, 1997, pp. 59-62).

sm-survival
Args: (x y status h &key (hv 0.05) (p 0.5) (status.code 1) (eval-points nil) (ngrid 50) (display "lines") (add nil) (xlab "x") (ylab "y") (color 'black) (type 'solid) (widht 1)) The required arguments are x, a list of covariate values, y, a list of survival times, and status, a list of indicator values defining the complete or censored survival times (by default an uncensored observation is identified by the status given through the :status.code argument). Also required is h, the smoothing parameter applied to the covariate scale. A normal kernel function is used and h is its standard deviation.
The smoothing parameter applied to the weighted Kaplan-Meier functions derived from the smoothing procedure in the covariate scale is provided by :hv. This ensures that a smooth estimate is obtained. The quantile to be estimated at each covariate value is given by :p, by default the median is estimated. The remaining arguments have the same meaning already discussed for the functions sm-density, smregression, etc.
This function returns a list whose elements are: the estimated quantiles at the covariate evaluation values; the covariate evaluation values; the smoothing parameter used for quantile estimation; the smoothing parameter used for the weighted Kaplan-Meier; the graph address of the plot.
Example 7.1 (Stanford data) As an example consider the data on the survival times of patients of different ages from the Stanford Heart Transplant Study. Thus, the response variable is taken to be the survival time on a log scale (Log.time) and the covariate is given the patient age in years (Age). The first call to the function sm-survival estimates the median of the survival times, while the following the 25th and 10th percentile respectively which are added to the first plot. A legend is also added to the plot, which is shown in Figure 13.

Extra functions accompanying the sm library
The file sm-miscfun.lsp provides several functions which can be of interest to Xlisp-Stat programmers and users. Most of them originated from existing S-Plus functions. They have been written during the development of the Xlisp-Stat version of the sm library, but hopefully they can result useful per se. In the following we briefly describe some of the most user-oriented, while the remaining functions can be read off directly from the source code.

Perspective plot
persp.s Args: (x1 x2 y &key (add nil) (spline 3) (axes nil) (frame T) (type 'solid) (color 'black) (variable-labels nil)) Creates a perspective plot given a matrix that represents heights on an evenly spaced grid (similar to a S-Plus function). x1 = vector containing x1 coordinates of the grid over which y is evaluated.
The values should be in sorted order. x2 = vector containing x2 coordinates of the grid over which y is evaluated.
The values should be in sorted order. y = matrix of size length(x1) by length(x2) giving the surface height at grid points, i.e. y [i,j] is evaluated at (x1[i], x2[j]). The rows of y are indexed by x1, and the columns by x2. add = if nil a new plot is returned, otherwise the perspective plot is added to the graph object given as argument.
axes = if axis is T then axes are plotted on the graph. frame = if frame is T then a frame box is drawn on the graph. variable-labels = if not NIL they provide the labels for X, Y, and Z axes. It must be a list of three strings.

Contour plot
contour.s Args: (x1 x2 y &key (add nil) (variable-labels nil) (nlevels 5) levels plevels (color 'black) (type 'solid) (width 1)) Creates a contour plot given a matrix that represents heights on an evenly spaced grid (similar to a S-Plus function). x1 = vector containing x1 coordinates of the grid over which y is evaluated. The values should be in sorted order. x2 = vector containing x2 coordinates of the grid over which y is evaluated.
The values should be in sorted order. y = matrix of size length(x1) by length(x2) giving the surface height at grid points, i.e. x1 [i,j] is evaluated at (x1[i], x2[j]). The rows of y are indexed by x1, and the columns by x2. add = if nil a new plot is returned, otherwise the contour is added to the graph object given as argument. variable-labels = if not NIL they provide the labels for X1, Y, and x2 axes.
It must be a list of three strings. nlevels = the approximate number of contour intervals desired. This is not needed if either levels or plevels are specified. levels = list of heights of contour lines. By default, nlevels lines are drawn which cover the range of y. plevels = a list of heights of contour lines expressed as the proportion of the maximum height of the surface. color = color for contour lines (see *colors*) type = line type, either 'solid or 'dashed width = line width, an integer value by default equal to 1 (normal line)

Adds or draws polygons to a plot
polygon Args: (x y &key (add nil) (frame nil) (color 'red) (title "Polygon Plot") variable-labels) Adds or draws a polygon with the specified vertices to a graph object.
X, Y = lists of values providing the coordinates of the vertices of a polygon. Each list must be ordered, i.e. the ith point is connected to the i+1st. ADD = if nil a new plot is returned, otherwise the polygon is added to the graph object given as argument. FRAME = if T a border is drawn around the polygonal region. COLOR = specifies the color of the area within the polygon. TITLE = a string for the plot title. VARIABLE-LABELS = a list of two strings providing labels for the axes.
The function returns a plot object.
:set-axis Meesage args: (axis &optional (min nil) (max nil) (n.ticks nil) (n.minors nil) (format nil)) This method allows formatting axes labels and tick marks. axis = 0 for x axis and 1 for y axis. min = the minimum of the axis max = the maximum of the axis n.ticks = the number of major tick marks n.minors = the number of minor tick marks format = a string specifying the format for labeling the major tick marks. Standard lisp format parameters can be used. Any values which are not to be changed should be specified as NIL.
A graph-proto method for adding a legend to a graph :legend Message args: (x y legend &key color type width symbol (frame T) (outer nil)) Adds a legend to the plot at the location (x,y), the top left corner of the rectangle. The content of the legend is specified by the argument legend, a list of string to be associated with color, line types, plotting symbols and width. If frame is T then the a box surrounding the legend is drawn.
To remove a legend send the :remove-legend message to a graph.

Functions related to matrices
var Args: (x &key (print nil) (varname nil)) If x is a matrix the function returns the sample variance-covariance matrix, if x is list or a vector it returns the sampling variance. If print is t then a nice printing of the covariance matrix is returned, else the matrix in lisp format. varname may be a list of strings with the name of the variables.
corr Args: (x &key (print nil) (varname nil) (r-squared nil)) Computes the correlation matrix of x, which can be a matrix or a list of lists. If r-square is T the squared correlation matrix is returned. If print is t then a nice printing of the correlation matrix is returned, else the matrix in lisp format. varname may be a list of strings with the name of the variables.
diag Args: (x) If x is a matrix, returns the diagonal of x. If x is a sequence, returns a diagonal matrix of rank equal to (length x) with diagonal elements equal to