1 PAST PAleontological STatistics Version 3.02 Reference manual Øyvind Hammer Natural History Museum University of Oslo ohammer@nhm.uio.no 1999-2014 2 Contents Welcome to the PAST! .......................................................................................................................... 10 Installation............................................................................................................................................. 11 The spreadsheet and the Edit menu ..................................................................................................... 12 Entering data..................................................................................................................................... 12 Selecting areas................................................................................................................................... 12 Moving a row or a column ................................................................................................................ 12 Renaming rows and columns ............................................................................................................ 13 Increasing the size of the array ......................................................................................................... 13 Cut, copy, paste................................................................................................................................. 13 Remove.............................................................................................................................................. 13 Row colors and symbols.................................................................................................................... 13 Selecting datatypes for columns, and specifying groups .................................................................. 13 Remove uninformative rows/columns.............................................................................................. 14 Transpose .......................................................................................................................................... 14 Grouped columns to multivar ........................................................................................................... 14 Grouped rows to multivar................................................................................................................. 14 Observations to contingency table ................................................................................................... 15 Stack grouped rows into columns..................................................................................................... 15 Value pairs to matrix ......................................................................................................................... 15 Samples to events (UA to RASC)........................................................................................................ 15 Events to samples (RASC to UA)........................................................................................................ 15 Loading and saving data.................................................................................................................... 15 Importing data from Excel................................................................................................................. 16 3 Reading and writing Nexus files ........................................................................................................ 16 Importing text files............................................................................................................................ 17 Counter.............................................................................................................................................. 17 Transform menu.................................................................................................................................... 18 Logarithm .......................................................................................................................................... 18 Subtract mean................................................................................................................................... 18 Remove trend.................................................................................................................................... 18 Row percentage................................................................................................................................. 18 Box-Cox.............................................................................................................................................. 18 Remove size from distances.............................................................................................................. 19 Abundance to presence/absence NOT YET IN PAST 3 ...................................................................... 20 Landmarks, Procrustes fitting............................................................................................................ 20 Landmarks, Bookstein fitting............................................................................................................. 21 Project to tangent space NOT YET IN PAST 3 .................................................................................... 21 Remove size from landmarks NOT YET IN PAST 3............................................................................. 21 Transform landmarks ........................................................................................................................ 21 Regular interpolation ........................................................................................................................ 22 Evaluate expression........................................................................................................................... 22 Plot menu .............................................................................................................................................. 24 Graph................................................................................................................................................. 24 XY graph............................................................................................................................................. 25 XY graph with error bars ................................................................................................................... 26 Histogram.......................................................................................................................................... 27 Bar chart/box plot ............................................................................................................................. 28 Percentiles......................................................................................................................................... 29 Missing values are deleted................................................................................................................ 29 Normal probability plot..................................................................................................................... 30 4 Ternary .............................................................................................................................................. 31 Bubble plot ........................................................................................................................................ 32 Matrix ................................................................................................................................................ 33 Surface............................................................................................................................................... 34 Statistics menu ...................................................................................................................................... 35 Univariate.......................................................................................................................................... 35 Two-sample tests .............................................................................................................................. 39 t test and related tests for equal means ....................................................................................... 39 F test for equal variances .............................................................................................................. 41 Mann-Whitney test for equal medians......................................................................................... 42 Kolmogorov-Smirnov test for equal distributions......................................................................... 43 Coefficient of variation (Fligner-Kileen test) ................................................................................. 45 Two-sample paired tests (t, sign, Wilcoxon) ..................................................................................... 47 Several-sample tests.......................................................................................................................... 49 One-way ANOVA ........................................................................................................................... 49 Kruskal-Wallis ................................................................................................................................ 52 Several-samples repeated measures tests........................................................................................ 54 Tukey’s pairwise post-hoc tests .................................................................................................... 55 Friedman test ................................................................................................................................ 55 Two-way ANOVA ............................................................................................................................... 58 Two-way ANOVA without replication ............................................................................................... 61 One-way ANCOVA ............................................................................................................................. 62 Correlation table ............................................................................................................................... 63 Intraclass correlation......................................................................................................................... 66 Normality tests.................................................................................................................................. 68 Contingency table (chi2 etc.) ............................................................................................................. 71 Risk/odds........................................................................................................................................... 73 5 Survival analysis (Kaplan-Meier curves, log-rank test etc.)............................................................... 75 Combine errors.................................................................................................................................. 77 Multivar menu....................................................................................................................................... 78 Principal components........................................................................................................................ 78 Principal coordinates......................................................................................................................... 82 Non-metric MDS................................................................................................................................ 83 Detrended correspondence analysis................................................................................................. 87 Canonical correspondence................................................................................................................ 88 Seriation ............................................................................................................................................ 89 CABFAC factor analysis, NOT YET IN PAST 3...................................................................................... 90 Discriminant analysis......................................................................................................................... 91 Two-block PLS.................................................................................................................................... 93 Cluster analysis.................................................................................................................................. 94 Neighbour joining.............................................................................................................................. 96 K-means clustering............................................................................................................................ 97 Multivariate normality ...................................................................................................................... 98 Box’s M.............................................................................................................................................. 99 MANOVA ......................................................................................................................................... 101 One-way ANOSIM............................................................................................................................ 102 One-way PERMANOVA.................................................................................................................... 103 Two-way ANOSIM............................................................................................................................ 104 Two-way PERMANOVA.................................................................................................................... 105 Mantel test and partial Mantel test................................................................................................ 106 SIMPER ............................................................................................................................................ 108 Paired Hotelling............................................................................................................................... 109 Modern Analog Technique.............................................................................................................. 110 Similarity and distance indices........................................................................................................ 112 6 Genetic sequence stats ................................................................................................................... 118 Model menu........................................................................................................................................ 119 Linear, bivariate............................................................................................................................... 119 Linear, multivariate (one independent, n dependent).................................................................... 123 Linear, multiple (one dependent, n independent).......................................................................... 124 Linear, multivariate multiple (m independent, n dependent) ........................................................ 125 Generalized Linear Model ............................................................................................................... 126 Polynomial regression..................................................................................................................... 128 Nonlinear......................................................................................................................................... 129 Linear........................................................................................................................................... 129 Quadratic..................................................................................................................................... 129 Power........................................................................................................................................... 130 Exponential.................................................................................................................................. 130 Von Bertalanffy............................................................................................................................ 130 Michaelis-Menten ....................................................................................................................... 130 Logistic......................................................................................................................................... 130 Gompertz..................................................................................................................................... 131 Gaussian ...................................................................................................................................... 131 Sinusoidal regression....................................................................................................................... 132 Smoothing spline............................................................................................................................. 134 LOESS smoothing............................................................................................................................. 135 Mixture analysis .............................................................................................................................. 136 Abundance models.......................................................................................................................... 138 Species packing (Gaussian) NOT YET IN PAST 3 .............................................................................. 140 Logarithmic spiral............................................................................................................................ 141 Diversity menu .................................................................................................................................... 142 Diversity indices............................................................................................................................... 142 7 Quadrat richness ............................................................................................................................. 144 Beta diversity................................................................................................................................... 147 Taxonomic distinctness ................................................................................................................... 149 Individual rarefaction ...................................................................................................................... 151 Sample rarefaction (Mao’s tau)....................................................................................................... 153 SHE analysis..................................................................................................................................... 154 Diversity permutation test .............................................................................................................. 155 Diversity t test ................................................................................................................................. 156 Diversity profiles.............................................................................................................................. 157 Time series menu ................................................................................................................................ 158 Simple periodogram........................................................................................................................ 158 REDFIT spectral analysis .................................................................................................................. 159 Multitaper spectral analysis ............................................................................................................ 161 Walsh transform.............................................................................................................................. 162 Short-time Fourier transform.......................................................................................................... 163 Wavelet transform .......................................................................................................................... 164 Autocorrelation ............................................................................................................................... 166 Autoassociation............................................................................................................................... 167 Cross-correlation............................................................................................................................. 169 Mantel correlogram (and periodogram)......................................................................................... 170 Runs test.......................................................................................................................................... 172 ARMA (and intervention analysis) NOT YET IN PAST 3.................................................................... 173 Insolation (solar forcing) model NOT YET IN PAST 3....................................................................... 175 Point events..................................................................................................................................... 176 Markov chain................................................................................................................................... 178 Simple smoothers............................................................................................................................ 179 FIR filter ........................................................................................................................................... 180 8 Date/time conversion...................................................................................................................... 182 Geometrical menu............................................................................................................................... 183 Directions (one sample) .................................................................................................................. 183 Directions (two samples)................................................................................................................. 186 Circular correlation.......................................................................................................................... 188 Spherical (one sample).................................................................................................................... 189 Point pattern analysis - nearest neighbours ................................................................................... 190 Ripley’s K point pattern analysis ..................................................................................................... 192 Kernel density.................................................................................................................................. 194 Point alignments.............................................................................................................................. 196 Spatial autocorrelation (Moran’s I)................................................................................................. 197 Gridding (spatial interpolation)....................................................................................................... 199 Multivariate allometry NOT YET IN PAST 3 ..................................................................................... 202 PCA of 2D landmarks (relative warps)............................................................................................. 203 Thin-plate splines for 2D landmarks................................................................................................ 204 Size from landmarks (2D or 3D) NOT YET IN PAST 3 ....................................................................... 205 Distance from landmarks (2D or 3D) NOT YET IN PAST 3................................................................ 206 All distances from landmarks (EDMA) NOT YET IN PAST 3 ............................................................. 207 Landmark linking NOT YET IN PAST 3.............................................................................................. 208 Fourier shape (2D) NOT YET IN PAST 3............................................................................................ 209 Elliptic Fourier shape analysis ......................................................................................................... 210 Hangle Fourier shape analysis......................................................................................................... 211 Eigenshape analysis NOT YET IN PAST 3.......................................................................................... 212 Coordinate transformation ............................................................................................................. 213 Stratigraphy menu............................................................................................................................... 214 Unitary Associations........................................................................................................................ 214 Ranking-Scaling................................................................................................................................ 219 9 Constrained optimization (CONOP)................................................................................................. 221 Appearance Event Ordination NOT YET IN PAST 3.......................................................................... 222 Diversity curve NOT YET IN PAST 3.................................................................................................. 222 Range confidence intervals NOT YET IN PAST 3.............................................................................. 222 Distribution-free range confidence intervals NOT YET IN PAST 3................................................... 222 Spindle diagram NOT YET IN PAST 3................................................................................................ 222 Cladistics.............................................................................................................................................. 222 Scripting............................................................................................................................................... 222 10 Welcome to the PAST! This program was originally designed as a follow-up to PALSTAT, a software package for paleontological data analysis written by P.D. Ryan, D.A.T. Harper and J.S. Whalley (Ryan et al. 1995). Through continuous development for more than ten years, PAST has grown into a comprehensive statistics package used not only by paleontologists, but in many fields of life science, earth science, engineering and economics. Further explanations of many of the techniques implemented together with case histories are found in the book “Paleontological data analysis” (Hammer & Harper 2005). If you have questions, bug reports, suggestions for improvements or other comments, we would be happy to hear from you. Contact us at ohammer@nhm.uio.no. For bug reports, remember to send us the data used, as saved from PAST, together with a complete description of the actions that lead to the problem. The latest version of Past, together with documentation and a link to the Past mailing list, are found at http://folk.uio.no/ohammer/past We are grateful if you cite PAST in scientific publications. The official reference is Hammer et al. (2001). References Hammer, Ø. & Harper, D.A.T. 2006. Paleontological Data Analysis. Blackwell. Hammer, Ø., Harper, D.A.T., and P. D. Ryan, 2001. PAST: Paleontological Statistics Software Package for Education and Data Analysis. Palaeontologia Electronica 4(1): 9pp. Harper, D.A.T. (ed.). 1999. Numerical Palaeobiology. John Wiley & Sons. 11 Installation The installation of PAST is easy: Just download the file 'Past3.exe' (unzipped) or ‘Past3.zip’ (zipped) and put it anywhere on your harddisk. Double-clicking the file will start the program. Windows will consider this a breach of security, and will ask if you trust the software provider. If you want to use the program, you will have to answer yes. We suggest you make a folder called 'past' anywhere on your harddisk, and put all the files in this folder. The lack of “formal” Windows installation is intentional, and allows installation without administrator privileges. On Mac, you should be able to just download the ‘dmg’ package and click on it to install. 12 The spreadsheet and the Edit menu PAST has a spreadsheet-like user interface. Data are entered as an array of cells, organized in rows (horizontally) and columns (vertically). Entering data To input data in a cell, click on the cell with the mouse and type in the data. The cells can also be navigated using the arrow keys. Any text can be entered in the cells, but most functions will expect numbers. Both comma (,) and decimal point (.) are accepted as decimal separators. Absence/presence data are coded as 0 or 1, respectively. Any other positive number will be interpreted as presence. Absence/presence-matrices can be shown with black squares for presences by ticking the 'Square mode' box above the array. Genetic sequence data are coded using C, A, G, T and U (lowercase also accepted). Missing data are coded with question marks (‘?'). Unless support for missing data is specifically stated in the documentation for a function, the function will not handle missing data correctly, so be careful. The convention in PAST is that items occupy rows, and variables columns. Three brachiopod individuals might therefore occupy rows 1, 2 and 3, with their lengths and widths in columns A and B. Cluster analysis will always cluster items, that is rows. For Q-mode analysis of associations, samples (sites) should therefore be entered in rows, while taxa (species) are in columns. For switching between Q-mode and R-mode, rows and columns can easily be interchanged using the Transpose operation. Selecting areas Most operations in PAST are only carried out on the area of the array which you have selected (marked). If you try to run a function which expects data, and no area has been selected, you will get an error message.  A row is selected by clicking on the row label (leftmost column).  A column is selected by clicking on the column label (top row).  Multiple rows are selected by selecting the first row label, then shift-clicking (clicking with the Shift key down) on the additional row labels.  Multiple columns are similarly marked by shift-clicking the additional column labels.  You can also select disjunct rows or columns by ctrl-clicking.  The whole array can be selected by clicking the upper left corner of the array (the empty grey cell) or by choosing 'Select all' in the Edit menu.  Smaller areas within the array can be selected by clicking and shift-clicking. Moving a row or a column 13 Select the ‘Drag rows/columns’ button in the ‘Click mode’ box. A row or a column can now be moved simply by clicking on the label and dragging to the new position. Renaming rows and columns When PAST starts, rows are numbered from 1 to 99 and columns are labelled A to Z. For your own reference, and for proper labelling of graphs, you should give the rows and columns more descriptive but short names. Select the 'Row attributes' option above the spreadsheet to see an editable column of the row names. Select the ‘Column attributes’ option to see an editable row of the column names. Increasing the size of the array By default, PAST has 99 rows and 26 columns. If you should need more, you can add rows or columns by choosing 'Insert more rows' or 'Insert more columns' in the Edit menu. Rows/columns will be inserted after the marked area, or at the bottom/right if no area is selected. When loading large data files, rows and/or columns are added automatically as needed. Cut, copy, paste The cut, copy and paste functions are found in the Edit menu. You can cut/copy data from the PAST spreadsheet and paste into other programs, for example Word and Excel. Likewise, data from other programs can be pasted into PAST – these need to be in a tab-separated text format. Before pasting, select the top left cell of the spreadsheet area in Past you want to paste into. Take care not to paste into the possibly hidden column and row attribute fields, unless you mean to. Remove The remove function (Edit menu) allows you to remove selected row(s) or column(s) from the spreadsheet. The removed area is not copied to the paste buffer. Row colors and symbols Each row can be given a color and a symbol (dot, cross, square etc.). These will be used in scatter plots and other plots. Select the ‘Row attributes’ option to edit the rows and colors individually, or use the ‘Row colors/symbols’ function to set all selected rows simultaneously (optionally based on the group, see below). Selecting datatypes for columns, and specifying groups Each column can be given a datatype using the ‘Column attributes’ mode: Unspecified (-) This is the default datatype. Ordinal, nominal or binary 14 Specifying one of these types is only required if you wish to use mixed similarity/distance measures. Group In a group column, you can enter identifiers for groups of data. You can use integers, or strings such as ‘males’ and ‘females’ (without the apostrophes). This will allow group-based polygons or ellipses in scatter plots. A group column is also required for many analyses, such as MANOVA. It is recommended to have rows in the same group as consecutive. Some analyses (e.g. two-way ANOVA) require two or even more group columns. Note that unlike previous versions of Past, there are no automatic links between colors, symbols and groups. If you wish to use different colors and/or symbols for different groups, you can set up the group column first and then use the ‘Row colors/symbols’ function in the Edit menu to assign colors/symbols accordingly. Remove uninformative rows/columns Rows or columns can be uninformative especially with respect to multivariate analyses. Such rows and columns should be considered for removal. Several types can be searched for and removed: Rows or columns with only zeroes, rows or columns with only missing data (‘?’), and rows or columns with only one non-zero cell (singletons). Transpose The Transpose function, in the Edit menu, will interchange rows and columns. This is used e.g. for switching between R mode and Q mode in cluster analysis, principal components analysis and seriation. Grouped columns to multivar Converts from a format with multivariate items presented in consecutive groups of N columns to the Past format with one item per row and all variates along the columns. For N=2, two specimens and four variables a-d, the conversion is from a1 b1 a2 b2 c1 d1 c2 d2 to a1 b1 c1 d1 a2 b2 c2 d2 Grouped rows to multivar Converts from a format with multivariate items presented in consecutive groups of N rows to the Past format with one item per row and all variates along the columns. For N=2, two specimens and four variables a-d, the conversion is from 15 a1 b1 c1 d1 a2 b2 c2 d2 to a1 b1 c1 d1 a2 b2 c2 d2 Observations to contingency table Expects two columns of data. Each row is one observation. Each column contains categories coded as numbers, e.g. males=0, females=1 in the first column; European=1, African=3, Asian=5 in the second column. The occurrences of different combinations are counted, giving a contingency table which can be analysed with the ‘Contingency table’ module (Univariate menu). Stack grouped rows into columns Stacks groups horizontally along columns. This can be useful e.g. for performing univariate statistics on pairs of columns across groups. Value pairs to matrix Very similar to “Observations to contingency table”. Expects two columns of data, numbers or strings. Each row is one observation. Each column contains categories, e.g. Europe, Africa, Asia in the first column; Dogs, Cats, Foxes in the second column. The occurrences of different combinations are counted, giving a full data matrix, in this case with localities in columns and taxa in rows. Samples to events (UA to RASC) Given a data matrix of occurrences of taxa in a number of samples in a number of sections, as used by the Unitary Associations module, this function will convert each section to a single row with orders of events (FADs, LADs or both) as expected by the Ranking-Scaling module. Tied events (in the same sample) will be given equal ranking. Events to samples (RASC to UA) Expects a data matrix with sections/wells in rows, and taxa in columns, with FAD and LAD values in alternating columns (i.e. two columns per taxon). Converts to the Unitary Associations presence/absence format with sections in groups of rows, samples in rows and taxa in columns. Loading and saving data 16 The 'Open' function is in the File menu. You can also drag a file from the desktop onto the PAST window. PAST uses a text file format for easy importing from other programs (e.g. Word), as follows: The top left cell must contain a colon (:). Cells are tab-separated. There are two top rows with data types and column names, and three left columns with colors, symbols and row names. Here is an example: : - - - Group Slow Med Fast Species Black Dot North 4 2 3 0 Black Dot South 4 3 7 0 Red Dot West 18 24 33 1 Red Dot East 10 6 7 1 In addition to this format, Past can also detect and open files in the following formats:  Excel (only the first worksheet).  Nexus (see below), popular in systematics.  TPS format developed by Rohlf. The landmark, outlines, curves, id, scale and comment fields are supported, other fields are ignored.  NTSYS. Multiple tables and trees are not supported. The file must have the extension ‘.nts’.  FASTA molecular sequence format, simplified specification according to NCBI.  PHYLIP molecular sequence format. The file must have the extension ‘.phy’.  Arlequin molecular sequence format. For genotype data the two haplotypes are concatenated into one row. Not all options are supported.  BioGraph format for biostratigraphy (SAMPLES or DATUM format). If a second file with the same name but extension '.dct' is found, it will be included as a BioGraph dictionary.  RASC format for biostratigraphy. You must open the .DAT file, and the program expects corresponding .DIC and .DEP files in the same directory.  CONOP format for biostratigraphy. You must open the .DAT file (log file), and the program expects corresponding .EVT (event) and .SCT (section) files in the same directory. Importing data from Excel There are several ways to get data from Excel to Past.  Copy from Excel and paste into PAST. Make sure you click (select) the top left cell where the data should be placed in Past before pasting. This will depend on whether row or column attributes are included in the data.  Open the Excel file from PAST.  Save as tab-separated text in Excel. The resulting text file can be opened in PAST. Reading and writing Nexus files 17 The Nexus file format is used by many systematics programs. PAST can read and write the Data (character matrix) block of the Nexus format. Interleaved data are supported. Also, if you have performed a parsimony analysis and the 'Parsimony analysis' window is open, all shortest trees will be written to the Nexus file for further processing in other programs (e.g. MacClade or Paup). Note that not all Nexus options are presently supported. Importing text files Text files with values separated by white space, tabs or commas can be opened using the ‘Import text file’ function in the File menu. Counter A counter function is available in the Edit menu for use e.g. at the microscope when counting microfossils of different taxa. A single row (sample) must be selected. The counter window will open with a number of counters, one for each selected column (taxon). The counters will be initialized with the column labels and any counts already present in the spreadsheet. When closing the counter window, the spreadsheet values will be updated. Count up (+) or down (-) with the mouse, or up with the keys 0-9 and a-z (only the first 36 counters). The bars represent relative abundance. A log of events is given at the far right - scroll up and down with mouse or arrow keys. An optional auditive feedback has a specific pitch for each counter. 18 Transform menu These routines subject your data to mathematical operations. This can be useful for bringing out features in your data, or as a necessary preprocessing step for some types of analysis. Logarithm The Log function in the Transform menu log-transforms your data using the base-10 logarithm. If the data contain zero or negative values, it may be necessary to add a constant (e.g. 1) before logtransforming (use Evaluate Expression x+1). This is useful, for example, to compare your sample to a log-normal distribution or for fitting to an exponential model. Also, abundance data with a few very dominant taxa may be log-transformed in order to downweight those taxa. Missing data supported. Subtract mean This function subtracts the column mean from each of the selected columns. The means cannot be computed row-wise. Missing values supported. Remove trend This function removes any linear trend from a data set (two columns with X-Y pairs, or one column with Y values). This is done by subtraction of a linear regression line from the Y values. Removing the trend can be a useful operation prior to time series analyses such as spectral analysis, auto- and cross-correlation and ARMA. Missing data supported. Row percentage All values converted to the percentage of the row sum. Missing values supported. Box-Cox The Box-Cox transformation is a family of power transformations with the purpose of making data x more normally distributed. The transformation has a parameter : 19         0ln 0 1     x x y If the smallest input value is zero or negative (which would invalidate the transform), a constant is added to all data such that the minimum input value becomes 1. The default value of the parameter is calculated by maximizing the log likelihood function:      n i ix n L 1 2 ln1ˆln 2   , where σ2  is the variance of the transformed data. This optimal value can be changed by the user, limited to the range -4    4. Missing values supported. Remove size from distances Attempts to remove the size component from a multivariate data set of measured distances (specimens in rows, variables in columns). Three methods are available.  Isometric Burnaby’s method projects the set of measured distances onto a space orthogonal to the first principal component. Burnaby's method may (or may not!) remove isometric size from the data, for further "size-free" data analysis. Note that the implementation in PAST does not center the data within groups - it assumes that all specimens (rows) belong to one group.  Allometric Burnaby’s method will log-transform the data prior to projection, thus conceivably removing also allometric size-dependent shape variation from the data.  Allometric vs. standard estimates allometric coefficients with respect to a standard (reference) measurement L such as overall length (Elliott et al. 1995). This standard variable should be placed in the first column. Each additional column is regressed onto the first column after log-transformation, giving a slope (allometric coefficient) b for that variable. An adjusted measurement is then computed from the original value M as b adj L L MM        . Reference Elliott, N.G., K. Haskard & J.A. Koslow 1995. Morphometric analysis of orange roughy (Hoplostethus atlanticus) off the continental slope of southern Australia. Journal of Fish Biology 46:202-220. 20 Abundance to presence/absence NOT YET IN PAST 3 All positive (non-zero) numbers are replaced with ones. Missing values supported. Landmarks, Procrustes fitting Transforms your measured point coordinates to Procrustes coordinates. There is also a menu choice for Bookstein coordinates. Specimens go in different rows and landmarks along each row. If you have three specimens with four landmarks in 2D, your data should look as follows: x1 y1 x2 y2 x3 y3 x4 y4 x1 y1 x2 y2 x3 y3 x4 y4 x1 y1 x2 y2 x3 y3 x4 y4 For 3D the data will be similar, but with additional columns for z. Landmark data in this format could be analyzed directly with the multivariate methods in PAST, but it is recommended to standardize to Procrustes coordinates by removing position, size and rotation. A further transformation to Procrustes residuals (approximate tangent space coordinates) is achieved by selecting 'Subtract mean' in the Edit menu. You must convert to Procrustes coordinates first, then to Procrustes residuals. The “Rotate to major axis” option places the result into a standard orientation for convenience. The “Keep size” option adds a final step where the shapes are scaled back to their original centroid sizes. A thorough description of Procrustes and tangent space coordinates is given by Dryden & Mardia (1998). The algorithms for Procrustes fitting are from Rohlf & Slice (1990) (2D) and Dryden & Mardia (1998) (3D). It should be noted that for 2D, the iterative algorithm of Rohlf & Slice (1990) often gives slightly different results from the direct algorithm of Dryden & Mardia (1998). Past uses the former in order to follow the “industry standard”. Missing data is supported but only by column average substitution, which is perhaps not very meaningful. References Dryden, I.L. & K.V. Mardia 1998. Statistical Shape Analysis. Wiley. Rohlf, F.J. & Slice, D. 1990. Extensions of the Procrustes method for the optimal superimposition of landmarks. Systematic Zoology 39:40-59. 21 Landmarks, Bookstein fitting Bookstein fitting has a similar function as Procrustes fitting, but simply standardizes size, rotation and scale by forcing the two first landmarks onto the coordinates (0,0) and (1,0). It is not in common use today. Bookstein fitting is only implemented for 2D. Project to tangent space NOT YET IN PAST 3 After Procrustes or Bookstein fitting, some statistical procedures are ideally carried out on tangent space projected coordinates (usually it doesn’t make any difference, but don’t quote us on that!). With d the number of dimensions and p the number of landmarks, the projection is  ccdp XXIXX t  . Here, X is the nxdp matrix of n specimens, X’ is the transformed matrix, I the dpxdp identity matrix, and Xc the mean (consensus) configuration as a dp-element row vector. Remove size from landmarks NOT YET IN PAST 3 The 'Remove size from landmarks' option in the Transform menu allows you to remove size by dividing all coordinate values by the centroid size for each specimen (Procrustes coordinates are also normalized with respect to size). See Dryden & Mardia (1998), p. 23-26. Reference Dryden, I.L. & K.V. Mardia 1998. Statistical Shape Analysis. Wiley. Transform landmarks Allows rotation of the point cloud in steps of 90 degrees, and top-bottom or left-right flipping (mirroring), mainly for plotting convenience. The mirror operation may be useful for reducing a bilaterally symmetric landmark data, by Procrustes fitting the left half to a mirrored version of the right half (and optionally averaging the two). Only for 2D coordinates. 22 Regular interpolation Interpolates an irregularly sampled time series or transect (possibly multivariate) into a regular spacing, as required by many methods for time series analysis. The x values should be in the first selected column. These will be replaced by a regularly increasing series. All additional selected columns will be interpolated correspondingly. The perils of interpolation should be kept in mind. You can either specify the total number of interpolated points, or the new point spacing. Three interpolation methods are available. Evaluate expression This powerful feature allows flexible mathematical operations on the selected array of data. Each selected cell is evaluated, and the result replaces the previous contents. A mathematical expression must be entered, which can inlude any of the operators +, -, *, /, ^ (power), and mod (modulo). Also supported are brackets (), and the functions abs, atan, asin, cos, sin, exp, ln, sqrt, sqr, round and trunc. The following values are also defined:  x (the contents of the current cell)  l (the cell to the left if it exists, otherwise 0)  r (the cell to the right)  u (the cell above, or up)  d (the cell below, or down)  mean (the mean value of the current column)  min (the minimum value)  max (the maximum value)  n (the number of cells in the column)  i (the row index)  j (the column index)  random (uniform random number from 0 to 1)  normal (Gaussian random number with mean 0 and variance 1).  integral (running sum of the current column)  stdev (standard deviation of the current column)  sum (total sum of the current column) In addition, other columns can be referred to using the column name preceded by '%', for example %A. Examples: sqrt(x) Replaces all numbers with their square roots (x-mean)/stdev Mean and standard deviation normalization, column-wise x-0.5*(max+min) Centers the values around zero (u+x+d)/3 Three-point moving average smoothing x-u First-order difference 23 I Fills the column with the row numbers (requires non-empty cells, such as all zeros) sin(2*3.14159*i/n) Generates one period of a sine function down a column (requires non-empty cells) 5*normal+10 Random number from a normal distribution, with mean of 10 and standard deviation of 5. Missing values supported. 24 Plot menu Graph Plots one or more columns as separate graphs. The x coordinates are set automatically to 1,2,3,... There are four plot styles available: Graph (line), points, line with points, and bars. The 'Row labels' options sets the x axis labels to the appropriate row names. The “Log Y” option log-transforms the values to base 10. For values <=0, the log value is set to 0. The sequence can be smoothed with a 3-point moving average. Missing values are disregarded. 25 XY graph Plots one or more pairs of columns containing x/y coordinate pairs. The 'log Y' option log-transforms your Y values (zero or negative values are set to 0). The curve can also be smoothed using 3-point moving average. 95% concentration ellipses can be plotted in most scatter plots in PAST, such as scores for PCA, CA, DCA, PCO and NMDS. The calculation of these ellipses assumes bivariate normal distribution. They estimate a region where 95% of population points are expected to fall, i.e. they are not confidence regions for the mean. Convex hulls can also be drawn in the scatter plots, in order to show the areas occupied by points of different 'colors'. The convex hull is the smallest convex polygon containing all points. The minimal spanning tree is the set of lines with minimal total length, connecting all points. In the XY graph module, Euclidean lengths in 2D are used. Points with missing values in X and/or Y are disregarded. 26 XY graph with error bars As XY graph, but expects four columns (or a multiple), with x, y, x error and y error values. Symmetric error bars are drawn around each point, with half-width as specified. If an error value is set to zero or missing, the corresponding error bar is not drawn. Points with missing values in X and/or Y are disregarded. 27 Histogram Plots histograms (frequency distributions) for one or more columns. The number of bins is by default set to an "optimal" number (the zero-stage rule of Wand 1997):   31 349.1s,min49.3   nIQh where s is the sample standard deviation and IQ the interquartile range. The number of bins can be changed by the user. The "Fit normal" option draws a graph with a fitted normal distribution (Parametric estimation, not Least Squares). Kernel Density Estimation is a smooth estimator of the histogram. PAST uses a Gaussian kernel with range according to the rule given by Silverman (1986):   51 34.1,min9.0   nIQsh . Missing values are deleted. References Silverman, B.W. 1986. Density estimation for statistics and data analysis. Chapman & Hall. Wand, M.P. 1997. Data-based choice of histogram bin width. American Statistician 51:59-64. 28 Bar chart/box plot Bar or box plot for one or several columns (samples) of univariate data. Missing values are deleted. Bar chart For each sample, the mean value is shown by a bar. In addition, “whiskers” can optionally be shown. The whisker interval can represent a one-sigma or a 95% confidence interval (1.96 sigma) for the estimate of the mean (based on the standard error), or a one-sigma or 95% concentration interval (based on the standard deviation). Box plot For each sample, the 25-75 percent quartiles are drawn using a box. The median is shown with a horizontal line inside the box. The minimal and maximal values are shown with short horizontal lines ("whiskers"). If the "Outliers" box is ticked, another box plot convention is used. The whiskers are drawn from the top of the box up to the largest data point less than 1.5 times the box height from the box (the "upper inner fence"), and similarly below the box. Values outside the inner fences are shown as circles, values further than 3 times the box height from the box (the "outer fences") are shown as stars. The quartile methods (rounding or interpolation) are described under “Percentiles” below. Jitter plot Each value is plotted as a dot. To show overlapping points more clearly, they can be displaced using a random “jitter” value controlled by a slider. Bar chart Box plot 29 Percentiles For each percentile p, plots the value y such that p percent of the points are smaller than y. Two popular methods are included. For a percentile p, the rank is computed according to k=p(n+1)/100, and the value that corresponds to that rank taken. In the rounding method, k is rounded to the nearest integer, while in the interpolation method, non-integer ranks are handled by linear interpolation between the two nearest ranks. Missing values are deleted. 30 Normal probability plot Plots a normal probability (normal QQ) plot for one column of data. A normal distribution will plot on a straight line. For comparison, an RMA regression line is given, together with the Probability Plot Correlation Coefficient. (three groups were given in this example) Missing values are deleted. The normal order statistic medians are computed as N(i) = G(U(i)), where G is the inverse of the cumulative normal distribution function and U are the uniform order statistic medians:            ni nini inU iU n1 5.0 1,...3,2365.03175.0 1),(1 )( 31 Ternary Ternary plot for three columns of data, normally containing proportions of compositions. If a fourth column is included, it will be shown using either a bubble representation or as a color/grayscale map. Rows with missing value(s) in any column are deleted. When using the color map option, rows with only the fourth variable missing are included in the plot but do not contribute to the map. 32 Bubble plot Plotting 3D data (three columns) by showing the third axis as size of disks. Negative values are not shown. Select "Subtract min" to subtract the smallest third axis value from all values - this will force the data to be positive. The "Size" slider scales the bubbles relative to unit radius on the x axis scale. Rows with missing value(s) in any column are deleted. 33 Matrix Two-dimensional plot of the data matrix, using a grayscale with white for lowest value, black for highest, or a colour scale. Includes contouring. Use to get an overview over a large data matrix. Missing values are plotted as blanks (allowing holes and non-square boundaries). 34 Surface Three-dimensional landscape plot of a data matrix of elevation values. Colors are assigned according to height, and/or the surface can be gray-shaded using a lighting model. Vertical exaggeration is adjustable. Missing values are replaced with the average. The data in the example below are the same as for the matrix plot above. 35 Statistics menu Univariate This function computes a number of basic descriptive statistics for one or more samples (columns) of univariate data. Each sample must have at least 3 values. The columns can have different numbers of values. The following numbers are shown for each sample: N: The number of values n in the sample Min: The minimum value Max: The maximum value Mean: The estimate of the mean, calculated as n x x i Std. error: The standard error of the estimate of the mean, calculated as n s SEx  where s is the estimate of the standard deviation (see below). 36 Variance: The sample variance, calculated as      22 1 1 xx n s i . Stand. dev.: The sample standard deviation, calculated as      2 1 1 xx n s i . Median: The median of the sample. For n odd, the given value such that there are equally many values above and below. For n even, the average of the two central values. 25 prcntil: The 25th percentile, i.e. the given value such that 25% of the sample is below, 75% above. The “interpolation” method is used (see Percentile plot above). 75 prcntil: The 75th percentile, i.e. the given value such that 75% of the sample is below, 25% above. The “interpolation” method is used (see Percentile plot above). Skewness: The sample skewness, zero for a normal distribution, positive for a tail to the right. Calculated as     3 2 3 1 1 1)2)(1(                xx n xx nn n G i i . Note there are several versions of this around – Past uses the same equation as SPSS and Excel. Slightly different results may occur using other programs, especially for small sample sizes. Kurtosis: The sample kurtosis, zero for a normal distribution. Calculated as               32 1 3 1 1321 1 2 4 2 4 2                    nn n xx n xx nnn nn G i i . Again Past uses the same equation as SPSS and Excel. Geom. mean: The geometric mean, calculated as   n nxxx /1 21  . Logarithms are used internally. Coeff.var: Coefficient of variation, or ratio of standard deviation to the mean, in percent:   100 1 1 100 2      x xx n x s CV i 37 Bootstrapping Selecting bootstrapping will compute lower and upper limits for 95% confidence intervals, using the specified number of bootstrap replicates. Confidence intervals for the min and max values are not given, because bootstrapping is known to not work well for these statistics. Three different bootstrap methods are available (cf. Davison and Hinkley 1997): Simple (basic): The statistic estimated from the original sample is t. The simulated estimates from R bootstrap replicates are t1 * , t2 * , …, tR * . For a 95% CI, we set the one-tailed error =0.025. The simple (or basic) bootstrapped CI is then      * 1 * 11 2,2    RR tttt . To ensure integer-valued subscripts, values for R such as 999, 9999 or 99999 are convenient. Percentile: An even simpler estimate:      * 11 * 1 ,   RR tt . BCa (adjusted percentile method): This is a complex method, but somewhat more accurate than the simple and percentile bootstrap. Estimate a bias correction factor (called z in some texts):               1 * 1 R tt w r , where  is the cumulative normal function and || is the number of elements in the set. Note we use strictly less than, unlike some sources. Then calculate a skewness correction factor:     2 3 1 2 1 3 6                n i i n i i tt tt a , where t-i is the statistic computed with value i removed (jackknifed), and t is the mean of the jackknifed values. With these values for w and a, compute adjusted CI endpoint values 38           96.11 96.1 1 wa w wa           96.11 96.1 2 wa w wa , where 1.96 is the approximate quantile for the normal distribution corresponding to a 95% CI (the actual value used is 1.959964). The bootstrapped confidence interval is     * 1 * 1 21 , aRaR tt  . No interpolation is used if the index is not an integer. Missing data: Supported by deletion. 39 Two-sample tests A number of classical statistics and tests for comparing two univariate samples, as given in two columns. It is also possible to specify the two groups using a single column of values and an additional Group column. Missing data are disregarded. t test and related tests for equal means Sample statistics Means and variances are estimated as described above under Univariate statistics. The 95% confidence interval for the mean is based on the standard error for the estimate of the mean, and the t distribution. Normal distribution is assumed. With s the estimate of the standard deviation, the confidence interval is             n s tx n s tx nn 1,21,2 ,  . Here, t has n-1 degrees of freedom, and 1-α = 0.95 for a 95% confidence interval. The 95% confidence interval for the difference between the means accepts unequal sample sizes:     DdfDdf styxstyx ,2,2 ,   , 40 where    22   yyxxSSE ii    11 21  nndf dfSSEMSE / 21 11 2 nn nh   h D n MSE s 2  . The confidence interval is computed for the larger mean minus the smaller, i.e. the center of the CI should always be positive. The confidence interval for the difference in means is also estimated by bootstrapping (simple bootstrap), with the given number of replicates (default 9999). t test The t test has null hypothesis H0: The two samples are taken from populations with equal means. The t test assumes normal distributions and equal variances. From the standard error sD of the difference of the means given above, the test statistic is Ds yx t   . Unequal variance t test The unequal variance t test is also known as the Welch test. It can be used as an alternative to the basic t test when variances are very different, although it can be argued that testing for difference in the means in this case is questionable. The test statistic is 21 )Var()Var( nynx yx t    . The number of degrees of freedom is 41     1 )Var( 1 )Var( )Var()Var( 2 2 2 1 2 1 2 21            n ny n nx n y n x df . Monte Carlo permutation test The permutation test for equality of means uses the absolute difference in means as test statistic. This is equivalent to using the t statistic. The permutation test is non-parametric with few assumptions, but the two samples are assumed to be equal in distribution if the null hypothesis is true. The number of permutations can be set by the user. The power of the test is limited by the sample size – significance at the p<0.05 level can only be achieved for n>3 in each sample. Exact permutation test As the Monte Carlo permutation test, but all possible permutations are computed. Only available if the sum of the two sample sizes is less than 27. F test for equal variances The F test has null hypothesis 42 H0: The two samples are taken from populations with equal variance. Normal distribution is assumed. The F statistic is the ratio of the larger variance to the smaller. The significance is two-tailed, with n1 and n2 degrees of freedom. Monte Carlo and exact permutation tests on the F statistic are computed as for the t test above. Mann-Whitney test for equal medians The two-tailed (Wilcoxon) Mann-Whitney U test can be used to test whether the medians of two independent samples are different. It is a non-parametric test and does not assume normal distribution, but does assume equal-shaped distribution in both groups. The null hypothesis is H0: The two samples are taken from populations with equal medians. For each value in sample 1, count the number of values in sample 2 that are smaller than it (ties count 0.5). The total of these counts is the test statistic U (sometimes called T). If the value of U is smaller when reversing the order of samples, this value is chosen instead (it can be shown that U1+U2=n1n2). The program computes an asymptotic approximation to p based on the normal distribution (twotailed), which is only valid for large n. It includes a continuity correction and a correction for ties: 43  112 5.02 33 21 21              nn ffnnnn nnU z g gg where n=n1+n2 and fg is the number of elements in tie g. A Monte Carlo value based on the given number of random permutations (default 9999) is also given – the purpose of this is mainly as a control on the asymptotic value. For n1+n2<=30 (e.g. 15 values in each group), an exact p value is given, based on all possible group assignments. If available, always use this exact value. For larger samples, the asymptotic approximation is quite accurate. Kolmogorov-Smirnov test for equal distributions The Kolmogorov-Smirnov test is a nonparametric test for overall equal distribution of two univariate samples. In other words, it does not test specifically for equality of mean, variance or any other parameter. The null hypothesis is H0: The two samples are taken from populations with equal distribution. In the version of the test provided by Past, both columns must represent samples. You can not test a sample against a theoretical distribution (one-sample test). 44 The test statistic is the maximum absolute difference between the two empirical cumulative distribution functions:    xSxSD NN x 21 max  The algorithm is based on Press et al. (1992), with significance estimated after Stephens (1970). Define the function         1 21 22 12 j jj KS eQ   . With Ne = N1N2/(N1+N2), the significance is computed as   DNNQp eeKS 11.012.0  . The permutation test uses 10,000 permutations. Use the permutation p value for N<30 (or generally). References Press, W.H., Teukolsky, S.A., Vetterling, W.T. & Flannery, B.P. 1992. Numerical Recipes in C. 2nd edition. Cambridge University Press. Stephens, M.A. 1970. Use of the Kolmogorov-Smirnov, Cramer-von Mises and related statistics without extensive tables. Journal of the Royal Statistical Society, Series B 32:115-122. 45 Coefficient of variation (Fligner-Kileen test) This module tests for equal coefficient of variation in two samples. The coefficient of variation (or relative variation) is defined as the ratio of standard deviation to the mean in percent, and is computed as:   100 1 1 100 2      x xx n x s CV i . The 95% confidence intervals are estimated by bootstrapping (simple bootstrap), with the given number of replicates (default 9999). The null hypothesis if the statistical test is: H0: The samples were taken from populations with the same coefficient of variation. If the given p(normal) is less than 0.05, equal coefficient of variation can be rejected. Donnelly & Kramer (1999) describe the coefficient of variation and review a number of statistical tests for the comparison of two samples. They recommend the Fligner-Killeen test (Fligner & Killeen 1976), as implemented in Past. This test is both powerful and is relatively insensitive to distribution. The following statistics are reported: 46 T: The Fligner-Killeen test statistic, which is a sum of transformed ranked positions of the smaller sample within the pooled sample (see Donnelly & Kramer 1999 for details). E(T): The expected value for T. z: The z statistic, based on T, Var(T) and E(T). Note this is a large-sample approximation. p: The p(H0) value. Both the one-tailed and two-tailed values are given. For the alternative hypothesis of difference in either direction, the two-tailed value should be used. However, the Fligner-Killeen test has been used to compare variation within a sample of fossils with variation within a closely related modern species, to test for multiple fossil species (Donnelly & Kramer 1999). In this case the alternative hypothesis might be that CV is larger in the fossil population, if so then a one-tailed test can be used for increased power. The screenshot above reproduces the example of Donnelly & Kramer (1999), showing that the relative variation within Australopithecus afarensis is significantly larger than in Gorilla gorilla. This could indicate that A. afarensis represents several species. References Donnelly, S.M. & Kramer, A. 1999. Testing for multiple species in fossil samples: An evaluation and comparison of tests for equal relative variation. American Journal of Physical Anthropology 108:507- 529. Fligner, M.A. & Killeen, T.J. 1976. Distribution-free two sample tests for scale. Journal of the American Statistical Association 71:210-213. 47 Two-sample paired tests (t, sign, Wilcoxon) Three statistical tests (one parametric, two non-parametric) for two samples (columns) of univariate data. The data points are paired, meaning that the two values in each row are associated. For example, the test could be the for length of the left vs. the right arm in a number of people, or the diversity in summer vs. winter at a number of sites. Controlling for a “nuisance factor” (person, site) in this way increases the power of the test. The null hypothesis is: H0: The mean (t test) or median (sign test, Wilcoxon test) of the difference is zero. All reported p values are two-tailed. t test Testing for mean difference equal to zero using the standard one-sample t test on the differences. With di=xi-yi , we have      2 1 1 dd n s i , ns d t  . There are n-1 degrees of freedom. This test assumes normal distribution of the differences. 48 The exact version of the test calculates all possible group reassignments within pairs. It is only computed for less than 26 pairs. Sign test The sign (binomial) test simply counts the number of cases n1 where xi>yi and n2 where yi>xi.The number max(n1, n2) is reported. The p value is exact, computed from the binomial distribution. The sign test may have lower power than the other paired tests, but makes few assumptions. Wilcoxon signed rank test A non-parametric rank test that does not assume normal distribution. The null hypothesis is no median shift (no difference). All rows with zero difference are first removed by the program. Then the absolute values of the differences |di| are ranked (Ri), with mean ranks assigned for ties. The sum of ranks for pairs where di is positive is W+ . The sum of ranks for pairs where di is negative is W. The reported test statistic is W = max(W+ , W- ) (note that there are several other, equivalent versions of this test, reporting other statistics). For large n (say n>10), the large-sample approximation to p can be used. This depends on the normal distribution of the test statistic W: 4 )1( )(   nn WE      4824 121 3      g gg ff nnn WVar . The last term is a correction for ties, where fg is the number of elements in tie g. The resulting z is reported, together with the p value. The Monte Carlo significance value is based on 99,999 random reassignments of values to columns, within each pair. This value will be practically identical to the exact p value. For n<26, an exact p value is computed, by complete enumeration of all possible reassignments (there are 2n of them, i.e. more than 33 million for n=25). This is the preferred p value, if available. Missing data: Supported by deletion of the row. 49 Several-sample tests One-way ANOVA and Kruskal-Wallis tests for equality of means or medians between several univariate samples, given in separate columns. It is also possible to specify the groups using a single column of values and an additional Group column. Missing data are supported by deletion. One-way ANOVA One-way ANOVA (analysis of variance) is a statistical procedure for testing the null hypothesis that several univariate samples are taken from populations with the same mean. The samples are assumed to be close to normally distributed and have similar variances. If the sample sizes are equal, these two assumptions are not critical. If the assumptions are strongly violated, the nonparametric Kruskal-Wallis test should be used instead. ANOVA table The between-groups sum of squares is given by:    g Tgg xxn 2 bgSS , where ng is the size of group g, and the means are group and total means. The between-groups sum of squares has an associated dfbg , the number of groups minus one. 50 The within-groups sum of squares is    g i gi xx 2 wgSS where the xi are those in group g. The within-groups sum of square has an associated dfwg, the total number of values minus the number of groups. The mean squares between and within groups are given by bg bg bg df SS MS  wg wg wg df SS MS  Finally, the test statistic F is computed as wg bg MS MS F The p value is based on F with dfbg and dfwg degrees of freedom. Omega squared The omega squared is a measure of effect size, varying from 0 to 1: wgtotal wgbgbg2 MSSS MSdfSS    . Levene's test Levene's test for homogeneity of variance (homoskedasticity), that is, whether variances are equal as assumed by ANOVA, is also given. Two versions of the test are included. The original Levene's test is based on means. This version has more power if the distributions are normal or at least symmetric. The version based on medians has less power, but is more robust to non-normal distributions. Note that this test can be used also for only two samples, giving an alternative to the F test for two samples described above. Unequal-variance (Welch) ANOVA If Levene's test is significant, meaning that you have unequal variances, you can use the unequalvariance (Welch) version of ANOVA, with the F, df and p values given. 51 Analysis of residuals The “Residuals” tab shows properties of the residuals, in order to evaluate some assumptions of ANOVA such as normal and homoskedastic distribution of residuals. The Shapiro-Wilk test for normal distribution is given, together with several common plots of residuals (normal probability, residuals vs. group means, and histogram). Tukey’s pairwise post-hoc tests If the ANOVA shows significant inequality of the means (small p), you can go on to study the given table of "post-hoc" pairwise comparisons, based on the Tukey-Kramer test. The Studentized Range Statistic Q is given in the lower left triangle of the array, and the probabilities p(equal) in the upper right. n XX Q SL wgMS   , where XL is the larger and XS the smaller mean of the two samples being compared. If the sample sizes are not equal, their harmonic mean is used for n. Its significance is estimated according to Lund & Lund (1983), with dfwg degrees of freedom. 52 Kruskal-Wallis The Kruskal-Wallis test is a non-parametric ANOVA, comparing the medians of several univariate groups (given in columns). It can also be regarded as a multiple-group extension of the MannWhitney test (Zar 1996). It does not assume normal distribution, but does assume equal-shaped distribution for all groups. The null hypothesis is H0: The samples are taken from populations with equal medians. The test statistic H is computed as follows:  13 )1( 12 2             n n T nn H g g g , where ng is the number of elements in group g, n is the total number of elements, and Tg is the sum of ranks in group g. The test statistic Hc is adjusted for ties: nn ff H H i ii c      3 3 1 , 53 where fi is the number of elements in tie i. With G the number of groups, the p value is approximated from Hc using the chi-square distribution with G-1 degrees of freedom. This is less accurate if any ng<5. Mann-Whitney pairwise post-hoc tests Mann-Whitney pairwise test p values are given for all Np=G(G-1)/2 pairs of groups. The asymptotic approximation described under the Mann-Whitney module is used. If samples are very small, it may be useful to run the exact test available in that module instead. Four different views are available for the symmetric table: 1. Raw p values, uncorrected significance: The p values from each individual pairwise test, marked in pink if p<0.05, not corrected for multiple testing. 2. Raw p values, sequential Bonferroni significance: The p values from each individual pairwise test are shown uncorrected for multiple testing. Significance (pink marking) is assessed by first evaluating the smallest p value, with Bonferroni correction for Np pairs. If significant (pNp<0.05) the next smallest p value is significant if p(Np-1)<0.05, etc. 3. Bonferroni corrected p values: The values shown are p’ = pNp. Marked as significant if p’<0.05. 4. Mann-Whitney U: The test statistics. References Lund, R.E., Lund, J.R. 1983. Algorithm AS 190: Probabilities and upper quantiles for the studentized range. Journal of the Royal Statistical Society C 32:204-210. Zar, J.H. 1996. Biostatistical analysis. 3rd ed. Prentice Hall. 54 Several-samples repeated measures tests In repeated measures ANOVA, values in each row are observations on the same “subject”. Repeatedmeasures ANOVA is the extension of the paired t test to several samples. Each column (sample) must contain the same number of values. Missing values are not supported. The procedure begins like the independent-samples one-way ANOVA above. In short,    g Tg xxn 2 bgSS , where n is the sample size. The associated dfbg is the number of groups minus one.    g i gi xx 2 wgSS where the xi are those in group g. The associated dfwg is the total number of values minus the number of groups. The between-subjects sum of squares is    i Ti xxn 2 subSS , 55 where the ix are means of subject i across groups. The associated dfsub is the number of subjects minus one. The SSerror is simply SSwg – SSsub, with dferror = dfwg – dfsub. The mean squares are then the sum squares divided by their respective degrees of freedom: bg bg bg df SS MS  wg wg wg df SS MS  sub sub sub df SS MS  error error error df SS MS  . Finally, the F ratio is MSbg / MSerror, with dfbg and dferror degrees of freedom. Tukey’s pairwise post-hoc tests The "post-hoc" pairwise comparisons are based on the Tukey test. The Studentized Range Statistic Q is given in the lower left triangle of the array, and the probabilities p(equal) in the upper right. n XX Q SL errorMS   , where XL is the larger and XS the smaller mean of the two samples being compared. There are dferror degrees of freedom. Friedman test The Friedman test is a non-parametric test for equality of medians in several repeated-measures univariate groups. It can be regarded as the non-parametric version of repeated-measures ANOVA, or the repeated-measures version of the Kruskal-Wallis test. 56 The Friedman test follows Bortz et al. (2000). The basic test statistic is     k j j knT knk 1 22 )1(3 )1( 12  , where n are the number of rows, k the number of columns and Tj the column sums of the data table. The 2 value is then corrected for ties (if any):         m i ii tie tt knk 1 3 2 2 2 1 1 1   where m is the total number of tie groups and ti are the numbers of values in each tie group. For k=2, it is recommended to use one of the paired tests (e.g. sign or Wilcoxon test) instead. For small data sets where k=3 and n<10, or k=4 and n<8, the tie-corrected 2 value is looked up in a table of “exact” p values. When given, this is the preferred p value. The asymptotic p value (using the 2 distribution with k-1 degrees of freedom) is fairly accurate for larger data sets. It is computed from a continuity corrected version of 2 :   2 1 2 1          k j j kn TS      24 1112 32 2    kkn Skn  . 57 This 2 value is also corrected for ties using the equation above. The post hoc tests are by simple pairwise Wilcoxon, exact for n<20, asymptotic for n>=20. These tests have higher power than the Friedman test. Reference Bortz, J., Lienert, G.A. & Boehnke, K. 2000. Verteilungsfreie Methoden in der Biostatistik. 2nd ed. Springer. 58 Two-way ANOVA Two-way ANOVA (analysis of variance) tests the null hypotheses that several univariate samples have the same mean across each of two factors A and B, and that there are no dependencies (interactions) between factors. The samples are assumed to be close to normally distributed and have similar variances. If the sample sizes are equal, these two assumptions are not critical. The test assumes a fixed-factor design (the usual case). Three columns are needed: A group column (set data type to Group with ‘Column attributes’) with the levels for factor A, a group column with the levels for factor B, and a column of the corresponding measured values. The algorithm uses weighted means for unbalanced designs. Total sum of squares:    i i xx 2 TSS , taken over all points. The associated degrees of freedom dfT is the total number of values minus one. Within-group sum of squares:    1 2 21 2 wgSS g g i ggi xx , 59 where the xi are those in group (level) g1 for the first factor and g2 for the second factor, and the mean is taken within the same group combination. The associated dfwg is the total number of values minus the product of numbers of groups and columns. The between-groups sum of squares SSbg = SST – SSwg can be partitioned into three, namely the Factor A, Factor B and interaction terms.    i iA xxN 2 ASS , where the sum is over levels of Factor A, and the two means are the level mean and the total mean, respectively. NA is the number of levels of A. The degrees of freedom is dfA = NA – 1. Similarly, for Factor B:    j jB xxN 2 BSS where now the sum is over levels of Factor B. The degrees of freedom is dfB = NB – 1. The interaction sum of squares is SSAxB = SSbg – SSA – SSB, with dfAxB = (NA - 1)(NB - 1) degrees of freedom. Mean squares MS are the sum of squares divided by their respective degrees of freedom. Finally, the F ratios are FA = MSA / MSwg FB = MSB / MSwg FAxB = MSAxB / MSwg The graph of means is a simple graphical device, traditionally used to see the effects of factors and their interaction for a two-way ANOVA. The means are shown with either the A levels or the B levels on the x axis, and the other levels as separate lines: 60 Repeated measures (within-subjects) ANOVA Ticking the “Repeated measures” box selects another type of two-way ANOVA, where each of a number of “subjects” have received several treatments. NOT YET IN PAST 3 Missing values : Rows with missing values are deleted. 61 Two-way ANOVA without replication Two-way ANOVA for testing the null hypotheses that several univariate samples have the same mean across each of two factors. This module expects only one observation for each combination of levels for the two factors. The input data format is a table where the first factor levels enter in rows, and the second factor level in columns, e.g. a table of veterinary lab results: There is no interaction term. The equations are given by Ireland (2010), pp. 130-131. Reference Ireland, C.R. 2010. Experimental Statistics for Agriculture and Horticulture. CABI, 352 pp. 62 One-way ANCOVA ANCOVA (Analysis of Covariance) tests for equality of means for several univariate groups, adjusted for covariance with another variate. ANCOVA can be compared with ANOVA, but has the added feature that for each group, variance that can be explained by a specified "nuisance" covariate (x) is removed. This adjustment can increase the power of the test substantially. The program expects two or more pairs of columns, where each pair (group) is a set of correlated x-y data (means are compared for y, while x is the covariate). The example below uses three pairs (groups) a, b and c. The Plot tab presents a scatter plot and linear regression lines for all the groups. The ANOVA-like summary table contains sum-of-squares etc. for the adjusted means (between-groups effect) and adjusted error (within-groups), together with an F test for the adjusted means. An F test for the equality of regression slopes (as assumed by the ANCOVA) is also given. In the example, equal adjusted means in the three groups can be rejected at p<0.05. Equality of slopes can not be rejected (p=0.74). The Groups tab gives the summary statistics for each group (mean, adjusted mean and regression slope). Assumptions include similar linear regression slopes for all groups, normal distributions, similar variance and sample sizes. Missing data: x-y pairs with either x or y missing are disregarded. 63 Correlation table Two or more columns are required. A matrix is presented with the correlations between all pairs of columns. In the ‘Statistic \ p(uncorr)’ table format, correlation values are given in the lower triangle of the matrix, and the two-tailed probabilities that the columns are uncorrelated are given in the upper. Both parametric and non-parametric coefficients and tests are available. Missing data: Supported by pairwise deletion, except for partial correlation which uses mean value imputation. Linear r (Pearson) Pearson’s r is the most commonly used parametric correlation coefficient. It is given by       22 yyxx yyxx r i i i i ii      . The significance is computed using a two-tailed t test with n-2 degrees of freedom and 2 1 2 r n rt    . Spearman’s D and rs Spearman’s (non-parametric) rank-order correlation coefficient is the linear correlation coefficient (Pearson’s r) of the ranks. Following Press et al. (1992) it is computed as 64                                                nn gg nn ff ggffD nn r m mm k kk k m mmkk s 3 3 3 3 33 3 11 12 1 12 16 1 . Here, D is the sum squared difference of ranks (midranks for ties):    n i ii SRD 1 2 . The fk are the numbers of ties in the kth group of ties among the Ri’s, and the gm are the numbers of ties in the mth group of ties among the Si’s. For n>9, the probability of non-zero rs (two-tailed) is computed using a t test with n-2 degrees of freedom: 2 1 2 s s r n rt    . For small n this approximation is inaccurate, and for n<=9 the program therefore switches automatically to an exact test. This test compares the observed rs to the values obtained from all possible permutations of the first column. The asymptotic test on D is closely related to the test on rs (see Press et al. 1992). It is computed for all n (no exact test for small n). Kendall’s tau This non-parametric correlation coefficient is not in very common use. It is computed according to Press et al. (1992). All possible N(N-1)/2 pairs of bivariate data points are considered. If two pairs have the same direction in x as in y (x and y both decrease, or both increase), they are called concordant. If not, they are discordant. A tie in the x’s is called an extra-x, and a tie in the y’s is called an extra-y. Pairs with ties in both variables are discarded. The number of pairs in the four categories are counted. Then, extraxdiscordantconcordantextraydiscordantconcordant discordantconcordant    . The asymptotic test is based on Kendall’s tau being approximately normal, with zero mean and    19 104 var    NN N  . 65 Polyserial correlation This correlation is only carried out if the second column consists of integers with a range less than 100. It is designed for correlating a normally distributed continuous/interval variable (first column) with an ordinal variable (second column) that bins a normally distributed variable. For example, the second column could contain the numbers 1-3 coding for “small”, “medium” and “large”. There would typically be more “medium” than “small” or “large” values because of the underlying normal distribution of sizes. Past uses the two-step algorithm of Olsson et al. (1982). This is more accurate than their “ad hoc” estimator, and nearly as accurate as the full multivariate ML algorithm. The two-step algorithm was chosen because of speed, allowing a permutation test (but only for N<100; not yet in Past 3). For larger N the given asymptotic test (log-ratio test) is accurate. Partial linear correlation Using this option, for each pair of columns, the linear correlation is computed while controlling for all the remaining columns. For example, with three columns A, B, C the correlation AB is controlled for C; AC is controlled for B; BC is controlled for A. The partial linear correlation can be defined as the correlation of the residuals after regression on the controlling variable(s). The significance is estimated with a t test with n-2-k degrees of freedom, where k is the number of controlling variables: 2 1 2 r kn rt    References Olsson, U., F. Drasgow & N.J. Dorans. 1982. The polyserial correlation coefficient. Psychometrika 47:337-347. Press, W.H., S.A. Teukolsky, W.T. Vetterling & B.P. Flannery. 1992. Numerical Recipes in C. Cambridge University Press. 66 Intraclass correlation A typical use of the intraclass correlation coefficient (ICC) is to quantify rater reliability, i.e. level of agreement between several ‘raters’ measuring the same objects. It is a standard tool to assess measurement error. ICC=1 would indicate perfect reliability. Raters (or ‘judges’) go in columns, while the objects measured go in rows. In the example below there are four raters A-B, who have measured 6 objects. Past follows the standard reference, Shrout and Fleiss (1979), which provides a number of different coefficients, referred to as ICC(m,k) where m is the model type. If k=1, the coefficient evaluates individual measurements (by a single rater); otherwise it evaluates the average measurement across raters. The models are Model 1: the raters rating different objects are different, and randomly sampled from a larger set of raters Model 2: the same raters rate all objects, and the raters are a subset of a larger set of raters. Model 3: no assumptions about the raters. The most commonly used ICC is ICC(2,1), which is therefore marked in red in Past. 67 The analysis is based on a two-way ANOVA without replication, as described elsewhere in this manual. Confidence intervals are parametric, following the equations of Shrout and Fleiss (1979). The data in the example above are from the Shrout and Fleiss paper, the output from Past reproducing their results. Reference Shrout, P.E., Fleiss, J.L. 1979. Intraclass correlations: Uses in assessing rater reliability. Psychological Bulletin 86:420-428. 68 Normality tests Three statistical tests for normal distribution of one or several samples of univariate data, given in columns. The data used below were generated by the normal and uniform random number generators in Past (‘Evaluate expression’ module). For all three tests, the null hypothesis is H0: The sample was taken from a population with normal distribution. If the given p(normal) is less than 0.05, normal distribution can be rejected (marked in pink). Of the given tests, the Shapiro-Wilk and Anderson-Darling are considered to be the more exact, and the Jarque-Bera is given for reference. And even poorer test (four-bin chi-squared) was included in previous versions of Past. There is a maximum sample size of n=5000, while the minimum sample size is 3 (the tests will of course have extremely small power for such small n). Remember the multiple testing issue if you run these tests on several samples – a Bonferroni or other correction may be appropriate. Shapiro-Wilk test The Shapiro-Wilk test (Shapiro & Wilk 1965) returns a test statistic W, which is small for non-normal samples, and a p value. The implementation is based on the standard code “AS R94” (Royston 1995), correcting an inaccuracy in the previous algorithm “AS 181” for large sample sizes. Jarque-Bera test The Jarque-Bera test (Jarque & Bera 1987) is based on skewness S and kurtosis K. The test statistic is 69             4 3 6 2 2 K S n JB . In this context, the skewness and kurtosis used are     3 2 3 1 1              xx n xx n S i i ,     4 2 4 1 1              xx n xx n K i i . Note that these equations contain simpler estimators than the G1 and G2 given in the Univariate summary statistics module, and that the kurtosis here will be 3, not zero, for a normal distribution. Asymptotically (for large sample sizes), the test statistic has a chi-square distribution with two degrees of freedom, and this forms the basis for the p value given by Past. It is known that this approach works well only for large sample sizes, and Past therefore also includes a significance test based on Monte Carlo simulation, with 10,000 random values taken from a normal distribution. Anderson-Darling test The data Xi are sorted in ascending sequence, and normalized for mean and standard deviation:   ˆ ˆ  i i X Y . With F the normal cumulativedistribution function (CDF), the test statistic is          n i kni YFYFi n nA 1 1 2 1lnln12 1 . Significance is estimated according to Stephens (1986). First, a correction for small sample size is applied:        2 22* 25.275.0 1 nn AA . The p value is estimated as 70                         2.073.22314.101436.13exp1 6.02.0938.59796.42318.8exp1 6.034.038.1279.49177.0exp 6.00186.0709.52937.1exp 2*22*2* 2*22*2* 2*22*2* 2*22*2* Aaa AAa AAA AAA p Missing data: Supported by deletion. References Jarque, C. M. & Bera, A. K. 1987. A test for normality of observations and regression residuals. International Statistical Review 55:163–172. Royston, P. 1995. A remark on AS 181: The W-test for normality. Applied Statistics 44:547-551. Shapiro, S. S. & Wilk, M. B. 1965. An analysis of variance test for normality (complete samples). Biometrika 52:591–611. Stephens, M.A. 1986. Tests based on edf statistics. Pp. 97-194 in D'Agostino, R.B. & Stephens, M.A. (eds.), Goodness-of-Fit Techniques. New York: Marcel Dekker. 71 Contingency table (chi2 etc.) These tests expect a frequency table with numbers of elements in different categories (rows and columns). Rows represent the different states of one nominal variable, columns represent the states of another nominal variable, and cells contain the integer counts of occurrences of that specific state (row, column) of the two variables. The contingency table analysis then gives information on whether the two variables of taxon and locality are associated. For example, this test can be used to compare two samples (columns) with the number of individuals in each taxon organized in the rows. You should be cautious about this test if any of the cells contain less than five individuals (see Fisher’s exact test below). The significance of association between the two variables is given, with p values from the chi-squared distribution and from a permutation test with 9999 replicates. The "Sample vs. expected" box should be ticked if you have two columns, and your second column consists of counts from a theoretical distribution (expected values) with zero sampling error, possibly non-integer. This is not a small-sample correction. In this case, only the chi-squared test is available. The Monte Carlo permutation test uses the given number of random replicates. For "Sample vs. expected" these replicates are generated by keeping the expected values fixed, while the values in the first column are random with relative probabilities as specified by the expected values, and with constant sum. For two samples, all cells are random but with constant row and column sums. See e.g. Brown & Rothery (1993) or Davis (1986) for details. 72 The Fisher's exact test is also given (two-tailed). When available, the Fisher's exact test may be far superior to the chi-square. For large tables or large counts, the computation time can be prohibitive and will time out after one minute. In such cases the parametric test is probably acceptable in any case. The procedure is complex, and based on the network algorithm of Mehta & Patel (1986). Two further measures of association are given. Both are transformations of chi-squared (Press et al. 1992). With n the total sum of counts, M the number of rows and N the number of columns: Cramer’s V:  1N,1Mminn V 2    Contingency coefficient C: n C 2 2    Note that for nx2 tables, the Fisher’s exact test is available in the Chi^2 module. Missing data not supported. References Brown, D. & P. Rothery. 1993. Models in biology: mathematics, statistics and computing. John Wiley & Sons. Davis, J.C. 1986. Statistics and Data Analysis in Geology. John Wiley & Sons. Mehta, C.R. & N.R. Patel. 1986. Algorithm 643: FEXACT: a FORTRAN subroutine for Fisher's exact test on unordered r×c contingency tables. ACM Transactions on Mathematical Software 12:154-161. Press, W.H., S.A. Teukolsky, W.T. Vetterling & B.P. Flannery. 1992. Numerical Recipes in C. Cambridge University Press. 73 Risk/odds This module compares the counts of a binary outcome under two different treatments, with statistics that are in common use in medicine. The data are entered in a 2x2 table, with treatments in rows and counts of the two different outcomes in columns. The following example shows the results of a vaccination trial on 460 patients: Got influenza Did not get influenza Vaccine 20 220 Placebo 80 140 In general, the data take the form Outcome 1 Outcome 2 Treatment 1 d1 h1 Treatment 2 d0 h0 74 Let n1=d1+h1, n0=d0+h0 and p1=d1/n1, p0=d0/n0. The statistics are then computed as follows: Risk difference: RD = p1-p0 95% confidence interval on risk difference (Pearson’s chi-squared):     0 |00 1 11| 11 n pp n pp se     Interval: RD - 1.96 se to RD + 1.96 se Z test on risk difference (two-tailed): es RD z  Risk ratio: RR = p1/p0 95% confidence interval on risk ratio (“delta method”):   0011 1111 ln ndnd RRse  es eEF 96.1  Interval: RR / EF to RR x EF Z test on risk ratio (two-tailed): es RR z ln  Odds ratio: 00 11 hd hd OR  95% confidence interval on odds ratio (“Woolfs’s formula”):   0011 1111 ln hdhd ORse  es eEF 96.1  Interval: OR / EF to OR x EF Note there is currently no continuity correction. Missing data are not allowed and will give an error message. 75 Survival analysis (Kaplan-Meier curves, log-rank test etc.) Survival analysis for two groups (treatments) with provision for right censoring. The module draws Kaplan-Meier survival curves for the two groups and computes three different tests for equivalence. The program expects four columns. The first column contains times to failure (death) or censoring (failure not observed up to and including the given time) for the first group, the second column indicates failure (1) or censoring (0) for the corresponding individuals. The last two columns contain data for the second group. Failure times must be larger than zero. The program also accepts only one treatment (given in two columns), or more than two treatments in consecutive pairs of columns, plotting one or multiple Kaplan-Meier curves. The statistical tests are only comparing the first two groups, however. The Kaplan-Meier curves and the log-rank, Wilcoxon and Tarone-Ware tests are computed according to Kleinbaum & Klein (2005). Average time to failure includes the censored data. Average hazard is number of failures divided by sum of times to failure or censorship. The log-rank test is by chi-squared on the second group:                              j jjjj jjjjjjjj j jj nnnn mmnnmmnn em EO EO 1 var 21 2 21 21212121 2 22 22 2 222  . 76 Here, nij is the number of individuals at risk, and mij the number of failures, in group i at distinct failure time j. The expected number of failures in group 2 at failure time j is   jj jjj j nn mmn e 21 212 2    . The chi-squared has one degree of freedom. The Wilcoxon and Tarone-Ware tests are weighted versions of the log-rank test, where the terms in the summation formulas for O2-E2 and var(O2-E2) receive weights of nj and nj, respectively. These tests therefore give more weight to early failure times. They are not in common use compared with the log-rank test. This module is not strictly necessary for survival analysis without right censoring – the Mann-Whitney test may be sufficient for this simpler case. Missing data: Data points with missing value in one or both columns are disregarded. Reference Kleinbaum, D.G. & Klein, M. 2005. Survival analysis: a self-learning text. Springer. 77 Combine errors A simple module for producing a weighted mean and its standard deviation from a collection of measurements with errors (one sigma). Expects two columns: the data x and their one-sigma errors σ. The sum of the individual gaussian distributions is also plotted. The weighted mean and its standard deviation are computed as    i i i iix 2 2 1    ,   i i 2 1 1   . This is the maximum-likelihood estimator for the mean, assuming all the individual distributions are normal with the same mean. Missing data: Rows with missing data in one or both columns are deleted. 78 Multivar menu Principal components Principal components analysis (PCA) finds hypothetical variables (components) accounting for as much as possible of the variance in your multivariate data (Davis 1986, Harper 1999). These new variables are linear combinations of the original variables. PCA may be used for reduction of the data set to only two variables (the two first components), for plotting purposes. One might also hypothesize that the most important components are correlated with other underlying variables. For morphometric data, this might be size, while for ecological data it might be a physical gradient (e.g. temperature or depth). Bruton & Owen (1988) describe a typical morphometric application of PCA. The input data is a matrix of multivariate data, with items in rows and variates in columns. The PCA routine finds the eigenvalues and eigenvectors of the variance-covariance matrix or the correlation matrix, with the SVD algorithm. Use variance-covariance if all variables are measured in the same units (e.g. centimetres). Use correlation (normalized var-covar) if the variables are measured in different units; this implies normalizing all variables using division by their standard deviations. The eigenvalues give a measure of the variance accounted for by the corresponding eigenvectors (components). The percentages of variance accounted for by these components are also given. If most of the variance is accounted for by the first one or two components, you have scored a success, but if the variance is spread more or less evenly among the components, the PCA has in a sense not been very successful. In the example below (landmarks from gorilla skulls), component 1 is strong, explaining 45.9% of variance. The bootstrapped confidence intervals are not shown un less the ‘Bootstrap N’ value is non-zero. 79 Groups: If groups are specified with a group column, the PCA can optionally be carried out withingroup or between-group. In within-group PCA, the average within each group is subtracted prior to eigenanalysis, essentially removing the differences between groups. In between-group PCA, the eigenanalysis is carried out on the group means (i.e. the items analysed are the groups, not the rows). For both within-group and between-group PCA, the PCA scores are computed using vector products with the original data. Row-wise bootstrapping is carried out if a positive number of bootstrap replicates (e.g. 1000) is given in the 'Bootstrap N' box. The bootstrapped components are re-ordered and reversed according to Peres-Neto et al. (2003) to increase correspondence with the original axes. 95% bootstrapped confidence intervals are given for the eigenvalues. The 'Scree plot' (simple plot of eigenvalues) may also indicate the number of significant components. After this curve starts to flatten out, the components may be regarded as insignificant. 95% confidence intervals are shown if bootstrapping has been carried out. The eigenvalues expected under a random model (Broken Stick) are optionally plotted - eigenvalues under this curve may represent non-significant components (Jackson 1993). In the gorilla example above, the eigenvalues for the 16 components (blue line) lie above the broken stick values (red dashed line) for the first two components, although the broken stick is inside the 95% confidence interval for the second component. The 'View scatter' option shows all data points (rows) plotted in the coordinate system given by two of the components. If you have groups, they will be shown with different symbols and colours. The Minimal Spanning Tree is the shortest possible set of lines connecting all points. This may be used as a visual aid in grouping close points. The MST is based on an Euclidean distance measure of the original data points, and is most meaningful when all variables use the same unit. The 'Biplot' option shows a projection of the original axes (variables) onto the scattergram. This is another visualisation of the PCA loadings (coefficients) - see below. 80 If the "Eigenval scale" is ticked, the data points will be scaled by kd1 , and the biplot eigenvectors by kd - this is the correlation biplot of Legendre & Legendre (1998). If not ticked, the data points are not scaled, while the biplot eigenvectors are normalized to equal length (but not to unity, for graphical reasons) - this is the distance biplot. The 'View loadings' option shows to what degree your different original variables (given in the original order along the x axis) enter into the different components (as chosen in the radio button panel). These component loadings are important when you try to interpret the 'meaning' of the components. The 'Coefficients' option gives the PC coefficients, while 'Correlation' gives the correlation between a variable and the PC scores. If bootstrapping has been carried out, 95% confidence intervals are shown (only for the Coefficients option). Missing data can be handled by one of two methods: 1. Mean value imputation: Missing values are replaced by their column average. Not recommended. 2. Iterative imputation: Missing values are inititally replaced by their column average. An initial PCA run is then used to compute regression values for the missing data. The procedure is iterated until convergence. This is usually the preferred method, but can cause some overestimation of the strength of the components (see Ilin & Raiko 2010). 81 References Bruton, D.L. & A.W. Owen. 1988. The Norwegian Upper Ordovician illaenid trilobites. Norsk Geologisk Tidsskrift 68:241-258. Davis, J.C. 1986. Statistics and Data Analysis in Geology. John Wiley & Sons. Harper, D.A.T. (ed.). 1999. Numerical Palaeobiology. John Wiley & Sons. Ilin, A. & T. Raiko. 2010. Practical approaches to Principal Component Analysis in the presence of missing values. Journal of Machine Learning Research 11:1957-2000. Jackson, D.A. 1993. Stopping rules in principal components analysis: a comparison of heuristical and statistical approaches. Ecology 74:2204-2214. Peres-Neto, P.R., D.A. Jackson & K.M. Somers. 2003. Giving meaningful interpretation to ordination axes: assessing loading significance in principal component analysis. Ecology 84:2347-2363. 82 Principal coordinates Principal coordinates analysis (PCO) is another ordination method, also known as Metric Multidimensional Scaling. The algorithm is from Davis (1986). The PCO routine finds the eigenvalues and eigenvectors of a matrix containing the distances or similarities between all data points. The Gower measure will normally be used instead of Euclidean distance, which gives results similar to PCA. An additional eleven distance measures are available these are explained under Cluster Analysis. The eigenvalues, giving a measure of the variance accounted for by the corresponding eigenvectors (coordinates) are given for the first four most important coordinates (or fewer if there are fewer than four data points). The percentages of variance accounted for by these components are also given. The similarity/distance values are raised to the power of c (the "Transformation exponent") before eigenanalysis. The standard value is c=2. Higher values (4 or 6) may decrease the "horseshoe" effect (Podani & Miklos 2002). The 'View scatter' option allows you to see all your data points (rows) plotted in the coordinate system given by the PCO. If you have colored (grouped) rows, the different groups will be shown using different symbols and colours. The "Eigenvalue scaling" option scales each axis using the square root of the eigenvalue (recommended). The minimal spanning tree option is based on the selected similarity or distance index in the original space. Missing data is supported by pairwise deletion (not for the Raup-Crick, Rho or user-defined indices). References Davis, J.C. 1986. Statistics and Data Analysis in Geology. John Wiley & Sons. Podani, J. & I. Miklos. 2002. Resemblance coefficients and the horseshoe effect in principal coordinates analysis. Ecology 83:3331-3343. 83 Non-metric MDS Non-metric multidimensional scaling is based on a distance matrix computed with any of 21 supported distance measures, as explained under Similarity and Distance Indices above. The algorithm then attempts to place the data points in a two- or three-dimensional coordinate system such that the ranked differences are preserved. For example, if the original distance between points 4 and 7 is the ninth largest of all distances between any two points, points 4 and 7 will ideally be placed such that their euclidean distance in the 2D plane or 3D space is still the ninth largest. Non-metric multidimensional scaling intentionally does not take absolute distances into account. The program may converge on a different solution in each run, depending upon the initial conditions. Each run is actually a sequence of 11 trials, from which the one with smallest stress is chosen. One of these trials uses PCO as the initial condition, the others are random. The solution is automatically rotated to the major axes (2D and 3D). The algorithm implemented in PAST, which seems to work very well, is based on a new approach developed by Taguchi and Oono (2005). The minimal spanning tree option is based on the selected similarity or distance index in the original space. Environmental variables: It is possible to include one or more initial columns containing additional “environmental” variables for the analysis. These variables are not included in the ordination. The correlation coefficients between each environmental variable and the NMDS scores are presented as vectors from the origin. The length of the vectors are arbitrarily scaled to make a readable biplot, so only their directions and relative lengths should be considered. Shepard plot: This plot of obtained versus observed (target) ranks indicates the quality of the result. Ideally, all points should be placed on a straight ascending line (x=y). The R2 values are the coefficients of determination between distances along each ordination axis and the original distances 84 (perhaps not a very meaningful value, but is reported by other NMDS programs so is included for completeness). Missing data is supported by pairwise deletion (not for the Raup-Crick, Rho and user-defined indices). For environmental variables, missing values are not included in the computation of correlations. Reference Taguchi, Y.-H., Oono, Y. 2005. Relational patterns of gene expression via non-metric multidimensional scaling analysis. Bioinformatics 21:730-40. 85 Correspondence analysis Correspondence analysis (CA) is yet another ordination method, somewhat similar to PCA but for counted data. For comparing associations (columns) containing counts of taxa, or counted taxa (rows) across associations, CA is the more appropriate algorithm. Also, CA is more suitable if you expect that species have unimodal responses to the underlying parameters, that is they favour a certain range of the parameter, becoming rare for lower and higher values (this is in contrast to PCA, which assumes a linear response). The CA routine finds the eigenvalues and eigenvectors of a matrix containing the Chi-squared distances between all rows (or columns, if that is more efficient – the result is the same). The algorithm follows Greenacre (2010), with SVD. The eigenvalue, giving a measure of the similarity accounted for by the corresponding eigenvector, is given for each eigenvector. The percentages of similarity accounted for by these components are also given. The 'View scatter' option allows you to see all your data points (rows) plotted in the coordinate system given by the CA. If you have grouped rows, the different groups can be shown using separate convex hulls and concentration ellipses. In addition, the variables (columns, associations) can be plotted in the same coordinate system (Q mode), optionally including the column labels. If your data are 'well behaved', taxa typical for an association should plot in the vicinity of that association. Relay plot (NOT YET IN PAST 3): This is a composite diagram with one plot per column. The plots are ordered according to CA column scores. Each data point is plotted with CA first-axis row scores on the vertical axis, and the original data point value (abundance) in the given column on the horizontal axis. This may be most useful when samples are in rows and taxa in columns. The relay plot will then 86 show the taxa ordered according to their positions along the gradients, and for each taxon the corresponding plot should ideally show a unimodal peak, partly overlapping with the peak of the next taxon along the gradient (see Hennebert & Lees 1991 for an example from sedimentology). Missing data is supported by column average substitution. Reference Greenacre, M. 2010. Biplots in practice. Fundación BBVA, 237 pp. Hennebert, M. & A. Lees. 1991. Environmental gradients in carbonate sediments and rocks detected by correspondence analysis: examples from the Recent of Norway and the Dinantian of southwest England. Sedimentology 38:623-642. 87 Detrended correspondence analysis The Detrended Correspondence (DCA) module uses the same algorithm as Decorana (Hill & Gauch 1980), with modifications according to Oxanen & Minchin (1997). It is specialized for use on 'ecological' data sets with abundance data; samples in rows, taxa in columns. Eigenvalues for the four ordination axes are given as in CA, indicating their relative importance in explaining the spread in the data. Detrending is a sort of normalization procedure in two steps. The first step involves an attempt to 'straighten out' points lying in an arch, which is a common occurrence. The second step involves 'spreading out' the points to avoid clustering of the points at the edges of the plot. Detrending may seem an arbitrary procedure, but can be a useful aid in interpretation. Missing data is supported by column average substitution. References Hill, M.O. & H.G. Gauch Jr. 1980. Detrended Correspondence analysis: an improved ordination technique. Vegetatio 42:47-58. Oxanen, J. & P.R. Minchin. 1997. Instability of ordination results under changes in input data order: explanations and remedies. Journal of Vegetation Science 8:447-454. 88 Canonical correspondence Canonical Correspondence Analysis (Legendre & Legendre 1998) is correspondence analysis of a site/species matrix where each site has given values for one or more environmental variables (temperature, depth, grain size etc.). The ordination axes are linear combinations of the environmental variables. CCA is thus an example of direct gradient analysis, where the gradient in environmental variables is known a priori and the species abundances (or presence/absences) are considered to be a response to this gradient. Each site should occupy one row in the spreadsheet. The environmental variables should enter in the first columns, followed by the abundance data (the program will ask for the number of environmental variables). The implementation in PAST follows the eigenanalysis algorithm given in Legendre & Legendre (1998). The ordinations are given as site scores - fitted site scores are presently not available. Environmental variables are plotted as correlations with site scores. Both scalings (type 1 and 2) of Legendre & Legendre (1998) are available. Scaling 2 emphasizes relationships between species. Missing values are supported by column average substitution. Reference Legendre, P. & L. Legendre. 1998. Numerical Ecology, 2nd English ed. Elsevier, 853 pp. 89 Seriation Seriation of an absence-presence (0/1) matrix using the algorithm described by Brower & Kile (1988). This method is typically applied to an association matrix with taxa (species) in the rows and samples in the columns. For constrained seriation (see below), columns should be ordered according to some criterion, normally stratigraphic level or position along a presumed faunal gradient. The seriation routines attempt to reorganize the data matrix such that the presences are concentrated along the diagonal. There are two algorithms: Constrained and unconstrained optimization. In constrained optimization, only the rows (taxa) are free to move. Given an ordering of the columns, this procedure finds the 'optimal' ordering of rows, that is, the ordering of taxa which gives the prettiest range plot. Also, in the constrained mode, the program runs a 'Monte Carlo' simulation, generating and seriating 30 random matrices with the same number of occurences within each taxon, and compares these to the original matrix to see if it is more informative than a random one (this procedure is time-consuming for large data sets). In the unconstrained mode, both rows and columns are free to move. Missing data are treated as absences. Reference Brower, J.C. & K.M. Kile. 1988. Seriation of an original data matrix as applied to palaeoecology. Lethaia 21:79-93. 90 CABFAC factor analysis, NOT YET IN PAST 3 This module implements the classical Imbrie & Kipp (1971) method of factor analysis and environmental regression (CABFAC and REGRESS, see also Klovan & Imbrie 1971). The program asks whether the first column contains environmental data. If not, a simple factor analysis with Varimax rotation will be computed on row-normalized data. If environmental data are included, the factors will be regressed onto the environmental variable using the second-order (parabolic) method of Imbrie & Kipp, with cross terms. PAST then reports the RMA regression of original environmental values against values reconstructed from the transfer function. Different methods for cross-validation (leave-one-out and k-fold) are available. You can also save the transfer function as a text file that can later be used for reconstruction of palaeoenvironment (see below). This file contains:  Number of taxa  Number of factors  Factor scores for each taxon  Number of regression coefficients  Regression coefficients (second- and first-order terms, and intercept) Missing values are supported by column average substitution. References Imbrie, J. & N.G. Kipp. 1971. A new micropaleontological method for quantitative paleoclimatology: Application to a late Pleistocene Caribbean core. In: The Late Cenozoic Glacial Ages, edited by K.K. Turekian, pp. 71-181, Yale Univ. Press, New Haven, CT. Klovan, J.E. & J. Imbrie. 1971. An algorithm and FORTRAN-IV program for large scale Q-mode factor analysis and calculation of factor scores. Mathematical Geology 3:61-77. 91 Discriminant analysis This module provides discriminant analysis for two or more groups (the latter is sometimes called Canonical Variates Analysis). The groups must be specified with a group column. A scatter plot of specimens along the first two canonical axes produces maximal and second to maximal separation between all groups. The axes are linear combinations of the original variables as in PCA, and eigenvalues indicate amount of variation explained by these axes. If only two groups are given, a histogram is plotted instead. Missing data supported by column average substitution. Classifier Classifies the data, assigning each point to the group that gives minimal Mahalanobis distance to the group mean. The Mahalanobis distance is calculated from the pooled within-group covariance matrix, giving a linear discriminant classifier. The given and estimated group assignments are listed for each point. In addition, group assignment is cross-validated by a leave-one-out cross-validation (jackknifing) procedure. Mystery specimens: Rows with unknown group, i.e. ‘?’ In the group column, are not included in the discriminant analysis itself, but will be classified. In this way, it is possible to classify new specimens that are not part of the training set. 92 Confusion matrix A table with the numbers of points in each given group (rows) that are assigned to the different groups (columns) by the classifier. Ideally each point should be assigned to its respective given group, giving a diagonal confusion matrix. Off-diagonal counts indicate the degree of failure of classification. Computational details Different softwares use different versions of CVA. The computations used by Past are given below. Let B be the given data, with n items in rows and k variates in columns, centered on the grand means of columns (column averages subtracted). Let g be the number of groups, ni the number of items in group i. Compute the gxk matrix X of weighted means of within group residuals, for group i and variate j ijiij n BX  , where ijB is a column average within group i. Compute B2 from B by centering within groups. Now compute W and the normalized, pooled, within-group covariance matrix Wcov: 22BBW  WWcov gn   1 . e and U are the eigenvalues and eigenvectors of W; ec and Uc are the eigenvalues and eigenvectors of Wcov. Then, )1diag()1diag( eXUXUeZZ  . a and A are the eigenvalues and eigenvectors of Z’Z. We take only the first g-1 eigenvectors (columns of A), as the rest will be zero. The canonical variates are now  AeUC c1diag . The CVA scores are then BC. Reification of variables can be done along vectors WcovC. 93 Two-block PLS Two-block Partial Least squares can be seen as an ordination method comparable with PCA, but with the objective of maximizing covariance between two sets of variates on the same rows (specimens, sites). For example, morphometric and environmental data collected on the same specimens can be ordinated in order to study covariation between the two. The program will ask for the number of columns belonging to the first block. The remaining columns will be assigned to the second block. There are options for plotting PLS scores both within and across blocks, and PLS loadings. The algorithm follows Rohlf & Corti (2000). Permutation tests and biplots are not yet implemented. Partition the nxp data matrix Y into Y1 and Y2 (the two blocks), with p1 and p2 columns. The correlation or covariance matrix R of Y can then be partitioned as        2221 1211 RR RR R . The algorithm proceeds by singular value decomposition of the matrix R12 of correlations across blocks: t 2112 DFFR  . The matrix D contains the singular values i along the diagonal. F1 contains the loadings for block 1, and F2 the loadings for block 2 (cf. PCA). The "Squared covar %" is a measure of the overall squared covariance between the two sets of variables, in percent relative to the maximum possible (all correlations equal to 1) (Rohlf & Corti p. 741). The "% covar” of axes are the amounts of covariance explained by each PLS axis, in percents of the total covariance. They are calculated as  2 2 100 i i   . Missing data supported by column average substitution. Reference Rohlf, F.J. & M. Corti. 2000. Use of two-block partial least squares to study covariation in shape. Systematic Biology 49:740-753. 94 Cluster analysis The hierarchical clustering routine produces a 'dendrogram' showing how data points (rows) can be clustered. For 'R' mode clustering, putting weight on groupings of taxa, taxa should go in rows. It is also possible to find groupings of variables or associations (Q mode), by entering taxa in columns. Switching between the two is done by transposing the matrix (in the Edit menu). Three different algorithms are available:  Unweighted pair-group average (UPGMA). Clusters are joined based on the average distance between all members in the two groups.  Single linkage (nearest neighbour). Clusters are joined based on the smallest distance between the two groups.  Ward's method. Clusters are joined such that increase in within-group variance is minimized, One method is not necessarily better than the other, though single linkage is not recommended by some. It can be useful to compare the dendrograms given by the different algorithms, to informally assess the robustness of the clusters. For Ward's method, a Euclidean distance measure is inherent to the algorithm. For UPGMA and single linkage, the distance matrix can be computed using 24 different indices, as described under the ‘Similarity and distance indices’ section. Missing data: The cluster analysis algorithm can handle missing data, coded with question marks (?). This is done using pairwise deletion, meaning that when distance is calculated between two points, any variables that are missing are ignored in the calculation. For Raup-Crick, missing values are treated as absence. Missing data are not supported for Ward's method, nor for the Rho similarity measure. Two-way clustering: The two-way option allows simultaneous clustering in R-mode and Q-mode. 95 Stratigraphically constrained clustering: This option will allow only adjacent rows or groups of rows to be joined during the agglomerative clustering procedure. May produce strange-looking (but correct) dendrograms. Bootstrapping: If a number of bootstrap replicates is given (e.g. 100), the columns are subjected to resampling. Press Enter after typing to update the value in the “Boot N” number box. The percentage of replicates where each node is still supported is given on the dendrogram. Note on Ward’s method: PAST produces Ward’s dendrograms identical to those made by Stata, but somewhat different from those produced by Statistica. The reason for the discrepancy is unknown. 96 Neighbour joining Neigbour joining clustering (Saitou & Nei 1987) is an alternative method for hierarchical cluster analysis. The method was originally developed for phylogenetic analysis, but may be superior to UPGMA also for ecological data. In contrast with UPGMA, two branches from the same internal node do not need to have equal branch lengths. A phylogram (unrooted dendrogram with proportional branch lengths) is given. Distance indices and bootstrapping are as for other cluster analysis (above). To run the bootstrap analysis, type in the number of required bootstratp replicates (e.g. 1000, 10000) in the “Boot N” box and press Enter to update the value. Negative branch lengths are forced to zero, and transferred to the adjacent branch according to Kuhner & Felsenstein (1994). The tree is by default rooted on the last branch added during tree construction (this is not midpoint rooting). Optionally, the tree can be rooted on the first row in the data matrix (outgroup). Missing data supported by pairwise deletion. References Saitou, N. & M. Nei. 1987. The neighbor-joining method: a new method for reconstructing phylogenetic trees. Molecular Biology and Evolution 4:406-425. 97 K-means clustering K-means clustering (e.g. Bow 1984) is a non-hierarchical clustering method. The number of clusters to use is specified by the user, usually according to some hypothesis such as there being two sexes, four geographical regions or three species in the data set The cluster assignments are initially random. In an iterative procedure, items are then moved to the cluster which has the closest cluster mean, and the cluster means are updated accordingly. This continues until items are no longer "jumping" to other clusters. The result of the clustering is to some extent dependent upon the initial, random ordering, and cluster assignments may therefore differ from run to run. This is not a bug, but normal behaviour in k-means clustering. The cluster assignments may be copied and pasted back into the main spreadsheet, and corresponding colors (symbols) assigned to the items using the 'Numbers to colors' option in the Edit menu. Missing data supported by column average substitution. Reference Bow, S.-T. 1984. Pattern recognition. Marcel Dekker, New York. 98 Multivariate normality Multivariate normality is assumed by a number of multivariate tests. PAST computes Mardia's multivariate skewness and kurtosis, with tests based on chi-squared (skewness) and normal (kurtosis) distributions. A powerful omnibus (overall) test due to Doornik & Hansen (1994) is also given. If at least one of these tests show departure from normality (small p value), the distribution is significantly non-normal. Sample size should be reasonably large (>50), although a small-sample correction is also attempted for the skewness test. Missing data supported by column average substitution. References Doornik, J.A. & H. Hansen. 1994. An omnibus test for univariate and multivariate normality. W4&91 in Nuffield Economics Working Papers. Mardia, K.V. 1970. Measures of multivariate skewness and kurtosis with applications. Biometrika 36:519-530. 99 Box’s M Test for the equivalence of the covariance matrices for two or more multivariate samples marked with a group column. This is a test for homoscedasticity, as assumed by MANOVA. The Box's M statistic is given together with a significance value based on an F test. Note that this test is supposedly very sensitive. This means that a high p value will be a good, although informal, indicator of equality, while a highly significant result (low p value) may in practical terms be a somewhat too sensitive indicator of inequality. The statistic is computed as follows – note this equals the “-2 ln M” of some texts (Rencher 2002).    g i ingnM 1 ln1ln)( iSS where Si are the within-group covariance matrices, S is the pooled covariance matrix, g the number of groups, n the total number of rows, ni the number of rows in group i, and |•| denotes the determinant. For significance, with r the number of variates (columns), calculate                                         2 1 2 1 2 1 1 1 16 21 1 1 1 116 132 1 gnng rr gnngr rr g i i g i i   The degrees of freedom for the F test are then 100      2 1 2 1 1 2 211      df df rrgdf Finally, M df dfdf F 1 21   . The Monte Carlo test is based on 999 random permutations. Missing data supported by column average substitution. Reference Rencher, A.C. 2002. Methods of multivariate analysis, 2nd ed. Wiley. 101 MANOVA One-way MANOVA (Multivariate ANalysis Of VAriance) is the multivariate version of the univariate ANOVA, testing whether two or more groups (specified with a group column) have the same multivariate mean. Two statistics are provided: Wilk's lambda with it's associated Rao's F and the Pillai trace with it's approximated F. Wilk's lambda is probably more commonly used, but the Pillai trace may be more robust. Number of constraints: For correct calculation of the p values, the number of dependent variables (constraints) must be specified. It should normally be left at 0, but for Procrustes fitted landmark data use 4 (for 2D) or 6 (for 3D). Pairwise comparisons (post-hoc): If the MANOVA shows significant overall difference between groups, the analysis can proceed by pairwise comparisons. In PAST, the post-hoc analysis is simple, by pairwise Hotelling's tests. The following values can be displayed in the table:  Hotelling's p values, not corrected for multiple testing. Marked in pink if significant (p<0.05).  The same p values, but significance (pink) assessed using the sequential Bonferroni scheme.  Bonferroni corrected p values (multiplied by the number of pairwise comparisons). The Bonferroni correction gives very little power.  Squared Mahalanobis distances. Note: These pairwise comparisons use the within-group covariance matrix pooled over all groups participating in the MANOVA. They may therefore give slightly other results than if only two of the groups are selected for analysis. Missing data supported by column average substitution. 102 One-way ANOSIM ANOSIM (ANalysis Of Similarities) is a non-parametric test of significant difference between two or more groups, based on any distance measure (Clarke 1993). The distances are converted to ranks. ANOSIM is normally used for taxa-in-samples data, where groups of samples are to be compared. Items go in rows, variates in columns, and groups should be specified with a group column as usual. In a rough analogy with ANOVA, the test is based on comparing distances between groups with distances within groups. Let rb be the mean rank of all distances between groups, and rw the mean rank of all distances within groups. The test statistic R is then defined as   41   NN rr R wb . Large positive R (up to 1) signifies dissimilarity between groups. The one-tailed significance is computed by permutation of group membership, with 9,999 replicates (can be changed). Pairwise ANOSIMs between all pairs of groups are provided as a post-hoc test. Significant comparisons (at p<0.05) are shown in pink. The optional Bonferroni correction multiplies the p values with the number of comparisons. This correction is very conservative (produces large p values). The sequential Bonferroni option does not output corrected p values, but significance is decided based on step-down sequential Bonferroni, which is slightly more powerful than simple Bonferroni. Missing data supported by pairwise deletion (not for the Raup-Crick, Rho and user-defined indices). Reference Clarke, K.R. 1993. Non-parametric multivariate analysis of changes in community structure. Australian Journal of Ecology 18:117-143. 103 One-way PERMANOVA NPMANOVA (Non-Parametric MANOVA, also known as PERMANOVA) is a non-parametric test of significant difference between two or more groups, based on any distance measure (Anderson 2001). NPMANOVA is normally used for ecological taxa-in-samples data, where groups of samples are to be compared, but may also be used as a general non-parametric MANOVA. Items go in rows, variates in columns, and groups should be specified with a group column. NPMANOVA calculates an F value in analogy with ANOVA. In fact, for univariate data sets and the Euclidean distance measure, NPMANOVA is equivalent to ANOVA and gives the same F value. The significance is computed by permutation of group membership, with 9,999 replicates (can be changed by the user). Pairwise NPMANOVAs between all pairs of groups are provided as a post-hoc test. Significant comparisons (at p<0.05) are shown in pink. The Bonferroni correction shown in the upper triangle of the matrix multiplies the p values with the number of comparisons. This correction is very conservative (produces large p values). Missing data supported by pairwise deletion. Reference Anderson, M.J. 2001. A new method for non-parametric multivariate analysis of variance. Austral Ecology 26:32-46. 104 Two-way ANOSIM The two-way ANOSIM in PAST uses the crossed design (Clarke 1993). For more information see oneway ANOSIM, but two group columns are required. Reference Clarke, K.R. 1993. Non-parametric multivariate analysis of changes in community structure. Australian Journal of Ecology 18:117-143. 105 Two-way PERMANOVA The two-way NPMANOVA (Anderson, 2001) in PAST uses the crossed design. The design must be balanced, i.e. each combination of levels must contain the same number of rows. For more information see one-way NPMANOVA, but two group columns are required(as for two-way ANOSIM). Reference Anderson, M.J. 2001. A new method for non-parametric multivariate analysis of variance. Austral Ecology 26:32-46. 106 Mantel test and partial Mantel test The Mantel test (Mantel 1967, Mantel & Valand 1970) is a permutation test for correlation between two distance or similarity matrices. In PAST, these matrices can also be computed automatically from two sets of original data. The first matrix must be above the second matrix in the spreadsheet, and the rows be specified as two groups (with a group column). The two matrices must have the same number of rows. If they are distance or similarity matrices, they must also have the same number of columns. The R value is simply the Pearson’s correlation coefficient between all the entries in the two matrices (because the matrices are symmetric it is only necessary to correlate the lower triangles). It ranges from -1 to +1. The permutation test compares the original R to R computed in e.g. 9999 random permutations. The reported p value is one-tailed. In the example below, the first matrix (gpa) consists of Procrustes-fitted landmark coordinates from primate skulls, while the second matrix (seq) contains sequence data from the same primates. The user has selected the Euclidean measure for the first matrix, and Jukes-Cantor for the second. The two data sets seem to be negatively correlated (R=-0.19), and there is no significant positive correlation (the test is one-tailed). In other words, there is no correlation between morphology and genetics. 107 Partial Mantel test It is possible to add a third matrix C below the two matrices A and B as described above. This matrix must be marked as above, and contain the same number of rows as A and B. A separate similarity measure can then be selected for this matrix. If such a third matrix is included, the program will carry out a partial Mantel test for the correlation of A and B, controlling for similarities given in C (Legendre & Legendre 1998). Only matrix A is permutated, and the R value is computed as            22 11 BCAC BCACAB CAB RR RRR R    where R(AB) is the correlation coefficient between A and B. References Legendre, P. & L. Legendre. 1998. Numerical Ecology, 2nd English ed. Elsevier, 853 pp. Mantel, N. 1967. The detection of disease clustering and a generalized regression approach. Cancer Research 27:209-220. Mantel, N. & R.S. Valand 1970. A technique of nonparametric multivariate analysis. Biometrics 26:547-558. 108 SIMPER SIMPER (Similarity Percentage) is a simple method for assessing which taxa are primarily responsible for an observed difference between groups of samples (Clarke 1993). The overall significance of the difference is often assessed by ANOSIM. The Bray-Curtis similarity measure (multiplied with 100) is most commonly used with SIMPER, but the Euclidean, cosine and chord measures can also be used. If more than two groups are selected, you can either compare two groups (pairwise) by choosing from the lists of groups, or you can pool all samples to perform one overall multi-group SIMPER. In the latter case, all possible pairs of samples are compared using the Bray-Curtis measure. The overall average dissimilarity is computed using all the taxa, while the taxon-specific dissimilarities are computed for each taxon individually. Samples go in rows, grouped with a group column, and taxa in columns. In the output table, taxa are sorted in descending order of contribution to group difference. The last three columns show the mean abundance in each of the groups. Missing data supported by column average substitution. Reference Clarke, K.R. 1993. Non-parametric multivariate analysis of changes in community structure. Australian Journal of Ecology 18:117-143. 109 Paired Hotelling The paired Hotelling's test expects two groups of multivariate data, marked with a group column. Rows within each group must be consecutive. The first row of the first group is paired with the first row of the second group, the second row is paired with the second, etc. With n the number of pairs and p the number of variables:    2 2 2 )1( 1 1 1 T np pn F nT n n i iiy i i iii             ySy yYyYS Yy XXY 1 y T T 1 The F has p and n-p degrees of freedom. For n<=16, the program also calculates an exact p value based on the T2 statistic evaluated for all possible permutations. Missing data supported by column average substitution. 110 Modern Analog Technique The Modern Analog Technique is a calibration method for reconstructing a past environmental parameter (e.g. temperature) from faunal assosications. It works by finding modern sites with faunal associations close to those in downcore samples. Environmental data from the modern sites are then used to estimate the environment downcore. The (single) environmental variable, usually temperature, enters in the first column, and taxa in consecutive columns. All the modern sites, with known values for the environmental variable, go in the first rows, followed by all the downcore samples (these should have question marks in the environmental column). The plot on the first tab shows all the modern samples, with the observed temperature (for example) versus the MAT reconstructed temperature using leave-one-out cross-validation (jackknifing). Parameters to set:  Weighting: When several modern analogs are linked to one downcore sample, their environmental values can be weighted equally, inversely proportional to faunal distance, or inversely proportional to ranked faunal distance.  Distance measure: Several distance measures commonly used in MAT are available. "Squared chord" has become the standard choice in the literature.  Distance threshold: Only modern analogs closer than this threshold are used. A default value is given, which is the tenth percentile of distances between all sample pairs in the modern data. The "Dissimilarity distribution" histogram may be useful when selecting this threshold.  N analogs: This is the maximum number of modern analogs used for each downcore sample.  Jump method (on/off): For each downcore sample, modern samples are sorted by ascending distance. When the distance increases by more than the selected percentage, the subsequent modern analogs are discarded. 111 Note that one or more of these options can be disabled by entering a large value. For example, a very large distance threshold will never apply, so the number of analogs is decided only by the "N analogs" value and optionally the jump method. Cross validation The scatter plot and R2 value show the results of a leave-one-out (jackknifing) cross-validation within the modern data. The y=x line is shown in red. This only partly reflects the "quality" of the method, as it gives little information about the accuracy of downcore estimation. Dissimilarity distribution A histogram of all distances in the core-top (modern) data. Semivariogram Shows a semivariogram of variance in the environmental variable as a function of faunal difference. Several semivariogram models can be fitted. This type of plot is familiar from spatial geostatistics, but is also useful for MAT because it gives a good impression of the degree of “noise” in the faunal data with respect to environmental prediction. Reconstructions Reconstruction of the paleoenvironmental values using MAT. 112 Similarity and distance indices Computes a number of similarity or distance measures between all pairs of rows. The data can be univariate or (more commonly) multivariate, with variables in columns. The results are given as a symmetric similarity/distance matrix. This module is rarely used, because similarity/distance matrices are usually computed automatically from primary data in modules such as PCO, NMDS, cluster analysis and ANOSIM in Past. Euclidean Basic Euclidean distance (the value is adjusted for missing data).    i kijijk xxd 2 . Gower A distance measure that averages the difference over all variables, each term normalized for the range of that variable:     i sisi s kiji jk xx xx n d s minmax 1 . The Gower measure is similar to Manhattan distance (see below) but with range normalization. When using mixed data types (see below), this is the default measure for continuous and ordinal data. Chord Euclidean distance between normalized vectors. Commonly used for abundance data. Can be written as    i ki i ji i kiji jk xx xx d 22 22 . Manhattan The sum of differences in each variable:   i kijijk xxd . 113 Bray-Curtis Bray-Curtis is a popular similarity index for abundance data. Past calculates Bray-Curtis similarity as follows:       i kiji i kiji jk xx xx d 1 . This is algebraically equivalent to the form given originally by Bray and Curtis (1957):        i kiji i kiji jk xx xx d ,min 2 . Many authors operate with a Bray-Curtis distance, which is simply 1-d. Cosine The inner product of abundances each normalised to unit norm, i.e. the cosine of the angle between the vectors.    i ki i ji i kiji jk xx xx d 22 . Morisita For abundance data.              1 1 1 i ji i ji i jiji xx xx               1 1 2 i ki i ki i kiki xx xx        i ki i ji i kiji jk xx xx d 21 2  . 114 Horn Horn’s overlap index for abundance data (Horn 1966).  i jij xN  i kik xN          kkjjkjkj i i i kikijijikijikiji jk NNNNNNNN xxxxxxxx d lnlnln lnlnln       . Mahalanobis A distance measure taking into account the covariance structure of the data. With S the variancecovariance matrix:    kjkj xxSxx  1T jkd . Correlation The complement 1-r of Pearson’s r correlation across the variables:       22 1 kki i jji i kkijji jk xxxx xxxx d      . Taking the complement makes this a distance measure. See also the Correlation module, where Pearson’s r is given directly and with significance tests. Rho The complement 1-rs of Spearman’s rho, which is the correlation coefficient of ranks. See also the Correlation module, where rho is given directly and with significance tests. Dice Also known as the Sorensen coefficient. For binary (absence-presence) data, coded as 0 or 1 (any positive number is treated as 1). The Dice similarity puts more weight on joint occurences than on mismatches. 115 When comparing two rows, a match is counted for all columns with presences in both rows. Using M for the number of matches and N for the the total number of columns with presence in just one row, we have djk = 2M / (2M+N). Jaccard A similarity index for binary data. With the same notation as given for Dice similarity above, we have djk = M / (M+N). Geographical MISSING IN PAST 3 Distance in meters along a great circle between two points on the Earth’s surface. Exactly two variables (columns) are required, with latitudes and longitudes in decimal degrees (e.g. 58 degrees 30 minutes North is 58.5). Coordinates are expected in the WGS84 datum, and distance is calculated with respect to the WGS84 ellipsoid. Use of other datums will result in very slight errors. The accuracy of the algorithm used (Vincenty 1975) is on the order of 1 mm with respect to WGS84. Kulczynski A similarity index for binary data. With the same notation as given for Dice similarity above (with N1 and N2 referring to the two rows), we have 2 21 NM M NM M d jk     . Ochiai A similarity index for binary data, comparable to the cosine similarity for other data types: 21 NM M NM M d jk   . Simpson The Simpson index (Simpson 1943) is defined simply as M / Nmin, where Nmin is the smaller of the numbers of presences in the two rows. This index treats two rows as identical if one is a subset of the other, making it useful for fragmentary data. 116 Raup-Crick Raup-Crick index for absence-presence data. This index (Raup & Crick 1979) uses a randomization (Monte Carlo) procedure, comparing the observed number of species ocurring in both associations with the distribution of co-occurrences from 1000 random replicates from the pool of samples. Hamming Hamming distance for categorical data as coded with integers (or sequence data coded as CAGT). The Hamming distance is the number of differences (mismatches), so that the distance between (3,5,1,2) and (3,7,0,2) equals 2. In PAST, this is normalised to the range [0,1], which is known to geneticists as "p-distance". Jukes-Cantor Distance measure for genetic sequence data (CAGT). Similar to p (or Hamming) distance, but takes into account probability of reversals:        pd 3 4 1ln 4 3 Kimura The Kimura 2-parameter distance measure for genetic sequence data (CAGT). Similar to Jukes-Cantor distance, but takes into account different probabilities of nucleotide transitions vs. transversions (Kimura 1980). With P the observed proportion of transitions and Q the observed number of transversions, we have    QQPd 21ln 4 1 21ln 2 1  . Tajima-Nei Distance measure for genetic sequence data (CAGT). Similar to Jukes-Cantor distance, but does not assume equal nucleotide frequencies. User-defined similarity Expects a symmetric similarity matrix rather than original data. No error checking! User-defined distance Expects a symmetric distance matrix rather than original data. No error checking! 117 Mixed This option requires that data types have been assigned to columns (see Entering and manipulating data). A pop-up window will ask for the similarity/distance measure to use for each datatype. These will be combined using an average weighted by the number of variates of each type. The default choices correspond to those suggested by Gower, but other combinations may well work better. The "Gower" option is a range-normalised Manhattan distance. All-zeros rows: Some similarity measures (Dice, Jaccard, Simpson etc.) are undefined when comparing two all-zero rows. To avoid errors, especially when bootstrapping sparse data sets, the similarity is set to zero in such cases. Missing data: Most of these measures treat missing data (coded as ‘?’) by pairwise deletion, meaning that if a value is missing in one of the variables in a pair of rows, that variable is omitted from the computation of the distance between those two rows. The exceptions are rho distance, using column average substitution, and Raup-Crick, which treats missing data as zero. References Bray, J.R. & J.T. Curtis. 1957. An ordination of the upland forest communities of Southern Wisconsin. Ecological Monographs 27:325-349. Horn, H.S. 1966. Measurement of overlap in comparative ecological studies. American Naturalist 100:419-424. Kimura, M. 1980. A simple model for estimating evolutionary rates of base substitutions through comparative studies of nucleotide sequences. Journal of Molecular Evolution 16:111-120. Raup, D. & R.E. Crick. 1979. Measurement of faunal similarity in paleontology. Journal of Paleontology 53:1213-1227. Simpson, G.G. 1943. Mammals and the nature of continents. American Journal of Science 241:1-31. Vincenty, T. 1975. Direct and inverse solutions of geodesics on the ellipsoid with application of nested equations. Survey Review 176:88-93. 118 Genetic sequence stats A number of simple statistics on genetic sequence (DNA or RNA) data. The module expects a number of rows, each with a sequence. The sequences are expected to be aligned and of equal length including gaps (coded as ‘?’). Some of these statistics are useful for selecting appropriate distance measures elsewhere in Past. Total length: The total sequence length, including gaps, of one sequence Average gap: The number of gap positions, averaged over all sequences Average A, T/U, C, G: The average number of positions containing each nucleotide Average p distance: The p distance between two sequences, averaged over all pairs of sequences. The p (or Hamming) distance is defined as the proportion of unequal positions Average Jukes-Cantor d: The Jukes-Cantor d distance between two sequences, averaged over all pairs of sequences. d = -3ln(1 - 4p/3)/4, where p is the p distance Maximal Jukes-Cantor d: Maximal Jukes-Cantor distance between any two sequences Average transitions (P): Average number of transitions (a↔g, c↔t, i.e. within purines, pyrimidines) Average transversions (Q): Average number of transversions (a↔t, a↔c, c↔g, t↔g, i.e. across purines, pyrimidines) R=P/Q: The transition/transversion ratio Missing data: Treated as gaps. 119 Model menu Linear, bivariate If two columns are selected, they represent x and y values, respectively. If one column is selected, it represents y values, and x values are taken to be the sequence of positive integers (1,2,...). A straight line y=ax+b is fitted to the data. There are four different algorithms available: Ordinary Least Squares (OLS), Reduced Major Axis (RMA), Major Axis (MA), and Robust. OLS regression assumes the x values are fixed, and finds the line which minimizes the squared errors in the y values. Use this if your x values have very little error associated with them. RMA and MA try to minimize both the x and the y errors. RMA/MA fitting and standard error estimation is according to Warton et al. (2006). The “Robust” method is an advanced Model I (fixed x values) regression which is robust to outliers. It sometimes gives strange results, but can be very successful in the case of “almost” normally distributed errors but with some far-off values. The algorithm is “Least Trimmed Squares” based on the “FastLTS” code of Rousseeuw & Driessen (1999). Parametric error estimates are not available, but Past gives bootstrapped confidence intervals on slope and intercept (beware – this is extremely slow for large data sets). Both x and y values can be log-transformed (base 10), in effect fitting your data to the 'allometric' function y=10b xa . An a value around 1 indicates that a straight-line ('isometric') fit may be more applicable. The values for a and b, their errors, Pearson's r correlation, and the probability that the columns are not correlated are given. Note the r2 is simply the Pearson’s coefficient squared – it does not adjust for regression method. 120 The calculation of standard errors for slope and intercept assumes normal distribution of residuals and independence between the variables and the variance of residuals. If these assumptions are strongly violated, it is preferable to use the bootstrapped 95 percent confidence intervals (1999 replicates). The permutation test on correlation (r2 ) uses 9,999 replicates. Confidence band In OLS regression (not RMA/MA/Robust), a 95 percent "Working-Hotelling" confidence band for the fitted line (not for the data points!) is available. The confidence band is calculated as                  2 2 2 2,2/05.0 1 CI xx xx n SEtaxb i regn where the squared sum of residuals    22 iireg axbySE . When the intercept is forced to zero, the confidence band is calculated as   2 2 2 1,2/05.0CI i regn x x SEtax . Zero intercept Forces the regression line through zero. This has implications also for the calculation of slope and the standard error of the slope. All four methods handle this option. 121 Residuals The Residuals window reports the distances from each data point to the regression line, in the x and y directions. Only the latter is of interest when using ordinary linear regression rather than RMA or MA. The residuals can be copied back to the spreadsheet and inspected for normal distribution and independence between independent variable and residual variance (homoskedasticity). Durbin-Watson test The Durbin-Watson test for positive autocorrelation of residuals in y (violating an assumption of OLS regression) is given in the Residuals window. The test statistic varies from zero (total positive autocorrelation) through 2 (zero autocorrelation) to 4 (negative autocorrelation). For n<=400, an exact p value for no positive autocorrelation is calculated using the PAN algorithm (Farebrother 1980, with later corrections). The test is not accurate when using the Zero intercept option. Breusch-Pagan test The Breusch-Pagan test for heteroskedasticity, i.e. nonstationary variance of residuals (violating an assumption of OLS regression) is given in the Residuals window. The test statistic is LM = nr2 where r is the correlation coefficient between the x values and the squared residuals. It is asymptotically distributed as 2 with one degree of freedom. The null hypothesis of the test is homoskedasticity. Exponential functions Your data can be fitted to an exponential function y=eb eax by first log-transforming just your y column (in the Transform menu) and then performing a straight-line fit. RMA equations Slope      2 2 sign      xx yy ra . Standard error on   2 1 abs 2    n r aa . Intercept xayb  . Standard error on 22 2 a r sx n s b  , where sr is the estimate of standard deviation of residuals and sa is the standard error on slope. For zero intercept (b=0), set 0x and 0y for the calculation of slope and its standard error (including the calculation of r therein), and use n-1 instead of n-2 for the calculation of standard error. Missing data: Supported by row deletion. 122 References Farebrother, R.W. 1980. Pan's procedure for the tail probabilities of the Durbin-Watson statistic. Applied Statistics 29:224–227. Rousseeuw, P.J. & van Driessen, K. 1999. Computing LTS regression for large data sets. Institute of Mathematical Statistics Bulletin. Warton, D.I., Wright, I.J., Falster, D.S. & Westoby, M. 2006. Bivariate line-fitting methods for allometry. Biological Review 81:259-291. 123 Linear, multivariate (one independent, n dependent) When you have one independent variate and several dependent variates, you can fit each dependent variate separately to the independent variate using simple linear regression. This module makes the process more convenient by having a scroll button going through each dependent variate. The module expects two or more columns of measured data, with the independent in the first column and the dependents in consecutive columns. In addition, an overall MANOVA test of multivariate regression significance is provided. The Wilks' lambda test statistic is computed as the ratio of determinants HE E   , where E is the error (residuals) sum of squares and crossproducts, and H is the hypothesis (predictions) sum of squares and crossproducts. The Rao’s F statistic is computed from the Wilks’ lambda and subjected to a one-tailed F test (see ‘Linear, n independent, n dependent’ below). Missing data supported by column average substitution. Regression for geometric morhpometrics For Procrustes-fitted landmarks or Elliptic Fourier coefficients as the dependent variables, see the Geometry menu for regression with visualization of shape change. 124 Linear, multiple (one dependent, n independent) Two or more columns of measured data, with the dependent in the first column and the independents in consecutive columns. The program will present the multiple correlation coefficient R and R2 , together with the "adjusted" R2 and an overall ANOVA-type significance test. With SSR the regression sum of squares, SSE the error (residuals) sum of squares, n the number of points and k the number of independent variates, we have R2 =SSR/SST,    1 11 1 2 2    kn nR Radj , )1(   knSSE kSSR F . The coefficients (intercept, and slope for each independent variate) are presented with their estimated standard errors and t tests. Missing data supported by column average substitution. 125 Linear, multivariate multiple (m independent, n dependent) Requires two or more columns of measured data, with the dependent variables in the first column(s) and the independents in consecutive columns. The program will ask for the number of dependent variables. The output consists of four main parts. Overall MANOVA An overall test of multivariate regression significance. The Wilks' lambda test statistic is computed as the ratio of determinants HE E   , where E is the error (residuals) sum of squares and crossproducts, and H is the hypothesis (predictions) sum of squares and crossproducts. The Rao’s F statistic is computed from the Wilks’ lambda. With n the number of rows, p the number of dependent variables and q the number of independent variables, we have:   pq pqm F qp qp qp qpqnm 211 otherwise1 05if 5 4 1 2 1 1 1 1 22 22 22                    The F test has pq and m + 1- pq/2 degrees of freedom. Tests on independent variables The test for the overall effect of each independent variable (on all dependent variables) is based on a similar design as the overall MANOVA above, but comparing the residuals of regression with and without the independent variable in question. Tests on dependent variables See ‘Linear, n independent, one dependent’ above for details of the ANOVA tests for the overall effect of all independent variables on each dependent. Regression coefficients and statistics The complete set of coefficients and their significances for all combinations of independent and dependent variables. Missing data supported by column average substitution. 126 Generalized Linear Model This module computes a basic version of the Generalized Linear Model, for a single explanatory variable. It requires two columns of data (independent and dependent variables). GLM allows non-normal distributions, and also “transformation” of the model through a link function. Some particularly useful combinations of distribution and link function are: Normal distribution and the identity link: This is equivalent to ordinary least squares linear regression. Normal distribution and the reciprocal link: Fit to the function y=1/(ax+b). Normal or gamma distribution and the log link: Fit to the function y=exp(ax+b). Binomial (Bernoulli) distribution and the logit link: Logistic regression for a binary response variable (see figure above). Technical details The program uses the Iteratively Reweighted Least Squares (IRLS) algorithm for maximum likelihood estimation. The dispersion parameter φ, which is used only for the inference, not the parameter estimation, is fixed at φ =1 for the Poisson and binomial distributions. For the normal and gamma distributions, it is estimated using Pearson’s chi-square. The log-likelihood LL is computed from the deviance D by 2 D LL  . 127 The deviance is computed as follows: Normal:    i iiyD 2  Gamma:          i i ii i i yy D    ln2 Bernoulli:            i i i i i i i y y y yD  1 1 ln1ln2 (the first term defined as zero if yi=0) Poisson:          i ii i i i y y yD   ln2 The G statistic is the difference in D between the full model and an additional GLM run where only the intercept is fitted. G is approximately chi-squared with one degree of freedom, giving a significance for the slope. 128 Polynomial regression Two columns must be selected (x and y values). A polynomial of up to the fifth order is fitted to the data. The algorithm is based on a least-squares criterion and singular value decomposition (Press et al. 1992), with mean and variance standardization for improved numerical stability. The polynomial is given by 01 2 2 3 3 4 4 5 5 axaxaxaxaxay  . The chi-squared value is a measure of fitting error - larger values mean poorer fit. The Akaike Information Criterion has a penalty for the number of terms. The AIC should be as low as possible to maximize fit but avoid overfitting. R2 is the coefficient of determination, or proportion of variance explained by the model. Finally, a p value, based on an F test, gives the significance of the fit. Reference Press, W.H., S.A. Teukolsky, W.T. Vetterling & B.P. Flannery. 1992. Numerical Recipes in C. Cambridge University Press. 129 Nonlinear Attempts to fit two columns of x-y data to a number of nonlinear equations, using least squares. Select a function name from the list. To see more functions, grab a function name and drag up and down to scroll. The 95% confidence intervals are based on 1999 bootstrap replicates. Fitting to a nonlinear function can be a bit tricky. For most of the functions, Past uses an educated guess for the parameters, followed by Levenberg-Marquardt optimization. Note that the L-M algorithm has been improved from previous versions of Past, so results may differ slightly. The Akaike Information Criterion (AIC) may aid in the selection of model. Lower values for the AIC imply a better fit, adjusted for the number of parameters. Linear baxy  Included for comparison with the nonlinear functions. Fitting by ordinary least squares regression. Quadratic cbxaxy  2 Included for reference. Fitting by least-squares and SVD (the equation is linear in its coefficients). 130 Power caxy b  The usual power law equation. Initial guess by log-log transformation and linear regression (i.e. c = 0), followed by nonlinear optimization. Exponential caey bx  Initial guess by linearization (log-transforming y), followed by nonlinear optimization. Von Bertalanffy  cx beay   1 This equation is used for modelling growth of multi-celled animals (Brown & Rothery 1993). It is sometimes given in a slightly different form:    0 1 txK eLy    It is easy to see that aL  , K = c and   cbt ln0  . The value of a is first estimated by the maximal value of y, and b and c using a straight-line fit to a linearized model. Finally nonlinear optimization. Michaelis-Menten xb ax y   The Michaelis-Menten curve can make accurate fits to rarefaction curves, and may therefore (somewhat controversially) be used for extrapolating these curves to estimate biodiversity (Colwell & Coddington 1994). It is also an important model equation for chemical kinetics. The algorithm uses maximum-likelihood estimators for the so-called Eadie-Hofstee transformation (Raaijmakers 1987; Colwell & Coddington 1994), followed by nonlinear optimization. Logistic cx be a y    1 A sigmoidal (S-shaped) curve. The logistic equation can model growth with saturation (Brown & Rothery 1993), and was used by Sepkoski (1984) to describe the proposed stabilization of marine diversity in the late Palaeozoic. 131 The value of a is first estimated by the maximal value of y, and b and c using a straight-line fit to a linearized model. Finally nonlinear optimization. Gompertz cx be aey  Initial estimate is computed using regression on a linearized model , followed by nonlinear optimization. Gaussian   2 2 2c bx aey    The ‘bell curve’ with mean b and standard deviation c. Initial guess by fitting a parabola to ln y, followed by nonlinear optimization. References Brown, D. & P. Rothery. 1993. Models in biology: mathematics, statistics and computing. John Wiley & Sons. Colwell, R.K. & J.A. Coddington. 1994. Estimating terrestrial biodiversity through extrapolation. Philosophical Transactions of the Royal Society of London B 345:101-118. Press, W.H., S.A. Teukolsky, W.T. Vetterling & B.P. Flannery. 1992. Numerical Recipes in C. Cambridge University Press. Raaijmakers, J.G.W. 1987. Statistical analysis of the Michaelis-Menten equation. Biometrics 43:793- 803. Sepkoski, J.J. 1984. A kinetic model of Phanerozoic taxonomic diversity. Paleobiology 10:246-267. 132 Sinusoidal regression Two columns must be selected (x and y values). A sum of up to eight sinusoids with periods specified by the user, but with unknown amplitudes and phases, is fitted to the data. This can be useful for modeling periodicities in time series, such as annual growth cycles or climatic cycles, usually in combination with spectral analysis. The algorithm is based on a least-squares criterion and singular value decomposition. By default, the periods are set to the range of the x values, and harmonics (1/2, 1/3, 1/4, 1/5, 1/6, 1/7 and 1/8 of the fundamental period). These values can be changed, and need not be in harmonic proportion. The “Fit periods” option will sequentially optimize the period of each sinusoid (over the full meaningful range from one period to the Nyquist frequency), after subtracting all previously fitted sinusoids. This is a simple example of the “Matching pursuit” algorithm. The algorithm is slow but robust and will fairly reliably find the global optimum. The chi-squared value is a measure of fitting error - larger values mean poorer fit. The Akaike Information Criterion has a penalty for the number of sinusoids (the equation used assumes that the periods are estimated from the data). The AIC should be as low as possible to maximize fit but avoid overfitting. R2 is the coefficient of determination, or proportion of variance explained by the model. Finally, a p value, based on an F test, gives the significance of the fit. 133 It is not meaningful to specify periodicities that are smaller than two times the typical spacing of data points. Each sinusoid is given by y=a*cos(2*pi*(x-x0) / T - p), where a is the amplitude, T is the period and p is the phase. x0 is the first (smallest) x value. An overall constant offset (mean) is also given. There are also options to enforce a pure sine or cosine series, i.e. with fixed phases. 134 Smoothing spline Two columns must be selected (X and Y values). The data are fitted to a smoothing spline, which is a sequence of third-order polynomials continuous up to the second derivative. A typical application is the construction of a smooth curve going through a noisy data set. The algorithm follows de Boor (2001). Sharp jumps in your data can give rise to oscillations in the curve, and you can also get large excursions in regions with few data points. Multiple data points at the same X value are collapsed to a single point by weighted averaging and calculation of a combined standard deviation. An optional third columns specifies standard deviations on the data points. These are used for weighting the data. If unspecified, they are all set to 10% of the standard deviation of the Y values. The smoothing value set by the user is a normalized version of the smoothing factor of de Boor (default 1). Larger values give smoother curves. A value of 0 will start a spline segment at every point. Clicking "Optimize smoothing" will calculate an "optimal" smoothing by a crossvalidation procedure. "View given points" gives a table of the given data points X, Y and stdev(Y), the corresponding Y values on the spline curve (ys) and the residuals. The chi-squared test for each point may be used to identify outliers. The final column suggests an stdev(Y) value to use if forcing the p value to 0.5. An optional fourth input column (if used then the third column must also be filled with stdev values) may contain a different number of values from the previous columns. It contains X values to be used for interpolation between the data points. Optional columns 5-7 contain lower and upper limits for X values (rectangular distribution) and standard deviation for Y values (normal distribution), to be used by bootstrapping (Monte Carlo) simulation providing error bars for the interpolated values. These functions are included mainly for computing boundary ages for the geological time scale. Reference de Boor, Carl. 2001. A practical guide to splines. Springer. 135 LOESS smoothing Two columns must be selected (x and y values). The algorithm used is “LOWESS” (LOcally WEighted Scatterplot Smoothing; Cleveland 1979, 1981), with its recommended default parameters (including two robustness iterations). Given a number of points n and a smoothing parameter q specified by the user, the program fits the nq points around each given point to a straight line, with a weighting function decreasing with distance. The new smoothed point is the value of the fitted linear function at the original x position. The Bootstrap option will estimate a 95% confidence band for the curve based on 999 random replicates. In order to retain the structure of the interpolation, the procedure uses resampling of residuals rather than resampling of original data points. LOESS or smoothing spline? This is almost a matter of taste. Compare the curves above, for the same dataset. The spline often gives a more aesthetically pleasing curve because of its continuous derivatives, but can suffer from overshooting near sharp bends in the data. References Cleveland, W.S. 1979. Robust locally weighted fitting and smoothing scatterplots. Journal of the American Statistical Association 74:829-836. Cleveland, W.S. 1981. A program for smoothing scatterplots by robust locally weighted fitting. The American Statistician 35:54. 136 Mixture analysis Mixture analysis is a maximum-likelihood method for estimating the parameters (mean, standard deviation and proportion) of two or more univariate normal distributions, based on a pooled univariate sample. The program can also estimate mean and proportion of exponential and Poisson distributions. For example, the method can be used to study differences between sexes (two groups), or several species, or size classes, when no independent information about group membership is available. The program expects one column of univariate data, assumed to be taken from a mixture of normally distributed populations (or exponential or Poisson). In the example below, sizes of two brachiopod samples have been pooled into one sample. The means, standard deviations and proportions of the two original samples have been almost perfectly recovered. PAST uses the EM algorithm (Dempster et al. 1977), which can get stuck on a local optimum. The procedure is therefore automatically run 20 times, each time with new, random starting positions for the means. The starting values for standard deviation are set to s/G, where s is the pooled standard deviation and G is the number of groups. The starting values for proportions are set to 1/G. The user is still recommended to run the program a few times to check for stability of the solution ("better" solutions have less negative log likelihood values). The Akaike Information Criterion (AIC; Akaike 1974) is calculated with a small-sample correction: 1 )1(2 ln22AICc    kn kk Lk where k is the number of parameters, n the number of data points and L the likelihood of the model given the data. A minimal value for AIC indicates that you have chosen the number of groups that produces the best fit without overfitting. 137 It is possible to assign each of the data points to one of the groups with a maximum likelihood approach. This can be used as a non-hierarchical clustering method for univariate data. The “Assignments” button will open a window where the value of each probability density function is given for each data point. The data point can be assigned to the group that shows the largest value. Missing data: Supported by deletion. References Akaike, H. 1974. A new look at the statistical model identification. IEEE Transactions on Automatic Control 19: 716-723. Dempster, A.P., Laird, N.M. & Rubin, D.B. 1977. Maximum likelihood from incomplete data via the EM algorithm". Journal of the Royal Statistical Society, Series B 39:1-38. 138 Abundance models This module can be used for plotting taxon abundances in descending rank order on a linear or logarithmic (Whittaker plot) scale, or number of species in abundance octave classes (as shown when fitting to log-normal distribution). Taxa go in rows. It can also fit the data to one of four different standard abundance models:  Geometric, where the 2nd most abundant species should have a taxon count of k<1 times the most abundant, the 3rd most abundant a taxon count of k times the 2nd most abundant etc. for a constant k. With ni the count of the ith most abundant taxon, we have 1 1   i i knn . This will give a straight descending line in the Whittaker plot. Fitting is by simple linear regression of the log abundances.  Log-series, with two parameters alpha and x. The fitting algorithm is from Krebs (1989). The number of species with n individuals (this equation does not translate directly to the Whittaker plot representation): n x S n n    Broken stick (MacArthur 1957). There are no free parameters to be fitted in this model. With Stot the total number of species and ntot the total number of individuals:      iS j tottot tot i tot jSS n n 0 1 . 139  Log-normal. The fitting algorithm is from Krebs (1989). The logarithm (base 10) of the fitted mean and variance are given. The octaves refer to power-of-2 abundance classes: Octave Abundance 1 1 2 2-3 3 4-7 4 8-15 5 16-31 6 32-63 7 64-127 ... ... A significance value based on chi-squared is given for each of these models, but the power of the test is not the same for the four models and the significance values should therefore not be compared. It is important, as always, to remember that a high p value can not be taken to imply a good fit. A low value does however imply a bad fit. Also note that the chi-squared tests in Past do not seem to correspond with some other software, possibly because Past use counts rather than the logtransformed values in the Whittaker plots. References Krebs, C.J. 1989. Ecological Methodology. Harper & Row, New York. MacArthur, R.H. 1957. On the relative abundance of bird species. Proceedings of the National Academy of Sciences, USA 43:293-295. 140 Species packing (Gaussian) NOT YET IN PAST 3 This module fits Gaussian response models to species abundances along a gradient, for one or more species. The fitted parameters are optimum (average), tolerance (standard deviation) and maximum. One column of environmental measurements in samples (e.g. temperature), and one or more columns of abundance data (taxa in columns). The algorithm is based on weighted averaging according to ter Braak & van Dam (1989). Reference ter Braak, C.J.F & H. van Dam. 1989. Inferring pH from diatoms: a comparison of old and new calibration methods. Hydrobiologia 178:209-223. 141 Logarithmic spiral Fits a set of points in the plane to a logarithmic spiral. Useful for characterizing e.g. mollusc shells, teeth, claws and horns. Requires two columns of coordinates (x and y). The points must be given in sequence, either inwards or outwards. Left-handed and right-handed spirals are both acceptable. The fitted spiral in polar coordinates: b aer  . The scale a and the exponent b are given, together with the estimated center point, marked with a red cross. The whorl expansion rate W (factor increase in radius per whorl) is calculated from b as b eW 2  . The center position is estimated by nonlinear optimization and the spiral itself by linearization and regression. 142 Diversity menu Diversity indices These statistics apply to association data, where number of individuals are tabulated in rows (taxa) and possibly several columns (samples). The available statistics are as follows, for each sample:  Number of taxa (S)  Total number of individuals (n)  Dominance = 1-Simpson index. Ranges from 0 (all taxa are equally present) to 1 (one taxon dominates the community completely).         i i n n D 2 where ni is number of individuals of taxon i.  Simpson index 1-D. Measures 'evenness' of the community from 0 to 1. Note the confusion in the literature: Dominance and Simpson indices are often interchanged!  Shannon index (entropy). A diversity index, taking into account the number of individuals as well as number of taxa. Varies from 0 for communities with only a single taxon to high values for communities with many taxa, each with few individuals.  i ii n n n n H ln  Buzas and Gibson's evenness: eH /S 143  Brillouin’s index:     n nn HB i i  !ln!ln  Menhinick's richness index: n S  Margalef's richness index: (S-1) / ln(n)  Equitability. Shannon diversity divided by the logarithm of number of taxa. This measures the evenness with which individuals are divided among the taxa present.  Fisher's alpha - a diversity index, defined implicitly by the formula S=a*ln(1+n/a) where S is number of taxa, n is number of individuals and a is the Fisher's alpha.  Berger-Parker dominance: simply the number of individuals in the dominant taxon relative to n.  Chao1, bias corrected: An estimate of total species richness. Chao1 = S + F1(F1 - 1) / (2 (F2 + 1)), where F1 is the number of singleton species and F2 the number of doubleton species. Many of these indices are explained in Harper (1999). Approximate confidence intervals for all these indices can be computed with a bootstrap procedure. The given number of random samples (default 9999) are produced, each with the same total number of individuals as in the original sample. For each individual in the random sample, the taxon is chosen with probabilities proportional to the original abundances. A 95 percent confidence interval is then calculated. Note that the diversity in the replicates will often be less than, and never larger than, the pooled diversity in the total data set – this bias can optionally be “fixed” by centering the confidence interval on the original value. Bootstrapped comparison of diversity indices in two samples is provided in the Compare diversities module. Reference Harper, D.A.T. (ed.). 1999. Numerical Palaeobiology. John Wiley & Sons. 144 Quadrat richness Requires two or more columns, each containing presence/absence (1/0) of different taxa down the rows (positive abundance is treated as presence). Four non-parametric species richness estimators are included in PAST: Chao 2, first- and secondorder jackknife, and bootstrap. All of these require presence-absence data in two or more sampled quadrats of equal size. Colwell & Coddington (1994) reviewed these estimators, and found that the Chao2 and the second-order jackknife performed best. The output from Past is divided into two panels. First, the richness estimators and their analytical standard deviations (only for Chao2 and Jackknife1) are computed from the given set of samples. Then the estimators are computed from 1000 random resamplings of the samples with replacement (bootstrapping), and their means and standard deviations are reported. In other words, the standard deviations reported here are bootstrap estimates, and not based on the analytical equations. Chao2 The Chao2 estimator (Chao 1987) is calculated as in EstimateS version 8.2.0 (Colwell 2009), with bias correction:    12 11ˆ 2 11 2           Q QQ m m SS obsChao where Sobs is the total observed number of species, m the number of samples, Q1 the number of uniques (species that occur in precisely one sample) and Q2 the number of duplicates (species that occur in precisely two samples). If Q1>0 and Q2>0, variance is estimated as 145              4 2 2 12 2 1 2 2 2 2 11 2 2 11 2 14 11 14 121 12 11ˆraˆv                               Q QQQ m m Q QQ m m Q QQ m m SChao . If Q1>0 but Q2=0:       2 4 1 22 11 2 11 2 ˆ4 1 4 121 2 11ˆraˆv Chao Chao S Q m mQQ m mQQ m m S                           . If Q1=0:    obsobs SMSM obsChao eeSS   1ˆraˆv 2 , where M is the total number of occurrences of all species in all samples. Jackknife 1 First-order jackknife (Burnham & Overton 1978, 1979; Heltshe & Forrester 1983): 11 1ˆ Q m m SS obsjack         .                    S j jjack m Q fj m m S 0 2 12 1 1ˆraˆv , where fj is the number of samples containing j unique species. Jackknife 2 Second-order jackknife (Smith & van Belle 1984):      1 232 ˆ 2 21 2      mm mQ m mQ SS obsjack . No analytical estimate of variance is available. Bootstrap Bootstrap estimator (Smith & van Belle 1984):    obsS k m kobsboot pSS 1 1ˆ , where pk is the proportion of samples containing species k. No analytical estimate of variance is available. 146 References Burnham, K.P. & W.S. Overton. 1978. Estimation of the size of a closed population when capture probabilities vary among animals. Biometrika 65:623-633. Burnham, K.P. & W.S. Overton. 1979. Robust estimation of population size when capture probabilities vary among animals. Ecology 60:927-936. Chao, A. 1987. Estimating the population size for capture-recapture data with unequal catchability. Biometrics 43, 783-791. Colwell, R.K. & J.A. Coddington. 1994. Estimating terrestrial biodiversity through extrapolation. Philosophical Transactions of the Royal Society (Series B) 345:101-118. Heltshe, J. & N.E. Forrester. 1983. Estimating species richness using the jackknife procedure. Biometrics 39:1-11. Smith, E.P. & G. van Belle. 1984. Nonparametric estimation of species richness. Biometrics 40:119- 129. 147 Beta diversity Two or more rows (samples) of presence-absence (0/1) data, with taxa in columns. The beta diversity module in Past can be used for any number of samples (not limited to only two samples). The eight measures available are described in Koleff et al. (2003): Past Koleff et al. Equation Ref. Whittaker bw 1  S Whittaker (1960) Harrison b-1 1 1   N S  Harrison et al. (1992) Cody bc     2 HlHg  Cody (1975) Routledge bI                    i ii i ii T ee T T  101010 log 1 log 1 log Routledge (1977) Wilson-Shmida bt     2 HlHg  Wilson & Shmida (1984) Mourelle bme      12   N HlHg  Mourelle & Ezcurra (1997) Harrison 2 b-2 1 1 max   N S  Harrison et al. (1992) Williams b-3 S max 1   Williams (1996) 148 S: total number of species;  : average number of species; N: number of samples; g(H): total gain of species along gradient (samples ordered along columns); l(H): total loss of species; ei: number of samples containing species i; T: total number of occurrences. References Harrison, S., S.J. Ross & J.H. Lawton. 1992. Beta diversity on geographic gradients in Britain. Journal of Animal Ecology 61:151-158. Koleff, P., K.J. Gaston & J.J. Lennon. 2003. Measuring beta diversity for presence-absence data. Journal of Animal Ecology 72:367-382. Routledge, R.D. 1977. On Whittaker’s components of diversity. Ecology 58:1120-1127. Whittaker, R.H. 1960. Vegetation of the Siskiyou mountains, Oregon and California. Ecological Monographs 30:279-338. 149 Taxonomic distinctness Requires one or more columns (samples), each containing counts of individuals of different taxa down the rows. In addition, one or more group columns with names of genera/families etc. (see below). Taxonomic diversity and taxonomic distinctness as defined by Clarke & Warwick (1998), including confidence intervals computed from 1000 random replicates taken from the pooled data set (all samples). Note that the "global list" of Clarke & Warwick is not entered directly, but is calculated internally by pooling (summing) the given samples. These indices depend on taxonomic information also above the species level, which has to be entered for each species as follows. Species names go in the name column (leftmost; in the Row attributes), genus names in the first group column, family in second group column etc., up to four group columns. Of course you can substitute for other taxonomic levels as long as they are in ascending order. Species counts for the samples follow in the columns thereafter. Taxonomic distinctness in one sample is given by (note other, equivalent forms exist):         ji i iiji ji jiij xxxx xxw 21 , 150 where the wij are weights such that wij = 0 if i and j are the same species, wij = 1 if they are the same genus, etc. The x are the abundances. Taxonomic distinctness:      ji ji ji jiij xx xxw * . For presence-absence data, taxonomic diversity and distinctness will be valid but equal to each other. Reference Clarke, K.R. & Warwick, R.M. 1998. A taxonomic distinctness index and its statistical properties. Journal of Applied Ecology 35:523-531. 151 Individual rarefaction For comparing taxonomical diversity in samples of different sizes. Requires one or more columns of counts of individuals of different taxa (each column must have the same number of values). When comparing samples: Samples should be taxonomically similar, obtained using standardised sampling and taken from similar 'habitat'. Given one or more columns of abundance data for a number of taxa, this module estimates how many taxa you would expect to find in a sample with a smaller total number of individuals. With this method, you can compare the number of taxa in samples of different size. Using rarefaction analysis on a large sample, you can read out the number of expected taxa for any smaller sample size (including that of the smallest sample). The algorithm is from Krebs (1989), using a log Gamma function for computing combinatorial terms. An example application in paleontology can be found in Adrain et al. (2000). Let N be the total number of individuals in the sample, s the total number of species, and Ni the number of individuals of species number i. The expected number of species E(Sn) in a sample of size n and the variance V(Sn) are then given by                                s i i n n N n NN SE 1 1 152                                                                                                                        s j j i jiji s i ii n n N n N n NN n NN n N n NNN n N n NN n N n NN SV 2 1 1 1 2 1 Standard errors (square roots of resampling variances) are given by the program. In the graphical plot, these standard errors are converted to 95 percent confidence intervals. References Adrain, J.M., S.R. Westrop & D.E. Chatterton. 2000. Silurian trilobite alpha diversity and the endOrdovician mass extinction. Paleobiology 26:625-646. Krebs, C.J. 1989. Ecological Methodology. Harper & Row, New York. 153 Sample rarefaction (Mao’s tau) Sample rarefaction requires a matrix of presence-absence data (abundances treated as presences), with taxa in rows and samples in columns. Sample-based rarefaction (also known as the species accumulation curve) is applicable when a number of samples are available, from which species richness is to be estimated as a function of number of samples. PAST implements the analytical solution known as "Mao’s tau", with standard deviation. In the graphical plot, the standard errors are converted to 95 percent confidence intervals. See Colwell et al. (2004) for details. With H samples and Sobs the total number of observed species, let sj be the number of species found in j samples, such that s1 is the number of species found in exactly one sample, etc. The total number of species expected in h≤H samples is then     H j jjhobs sSh 1 ~  . The combinatorial coefficients α are               Hhj Hhj HjhH jHhH jh for0 for !! !!  . These coefficients are computed via a log Gamma function. The variance estimator is       H j jjh S h s 1 2 22 ~ ~ 1~   , where S ~ is an estimator for the unknown total species richness. Following Colwell et al. (2004), a Chao2-type estimator is used. For s2>0,   2 2 1 2 1~ Hs sH SS obs   . For s2=0,      12 11~ 2 11    sH ssH SS obs . For modeling and extrapolating the curve using the Michaelis-Menten equation, use the Copy Data button, paste to a new Past spreadsheet, and use the nonlinear fitting module in the Model menu. Reference Colwell, R.K., C.X. Mao & J. Chang. 2004. Interpolating, extrapolating, and comparing incidence-based species accumulation curves. Ecology 85:2717-2727. 154 SHE analysis SHE analysis (Hayek & Buzas 1997, Buzas & Hayek 1998) requires a matrix of integer abundance data (counts), with taxa in rows and samples in columns. The program calculates log species abundance (ln S), Shannon index (H) and log evenness (ln E = H – ln S) for the first sample. Then the second sample is added to the first, and the process continues. The resulting cumulative SHE profiles can be interpreted ecologically. If the samples are taken not from one homogenous population but across a gradient or up a stratigraphic section, breaks in the curve may be used to infer discontinuities (e.g. biozone boundaries). References Buzas, M.A. & L.-A. C. Hayek. 1998. SHE analysis for biofacies identification. The Journal of Foraminiferal Research 28:233-239. Hayek, L.-A. C. & M.A. Buzas. 1997. Surveying natural populations. Columbia University Press. 155 Diversity permutation test Expects two columns of abundance data with taxa down the rows.This module computes a number of diversity indices for two samples, and then compares the diversities using random permutations. 9999 random matrices with two columns (samples) are generated, each with the same row and column totals as in the original data matrix. 156 Diversity t test Comparison of the Shannon and Simpson diversities in two samples. Then Shannon t test is described by e.g. Hutcheson (1970), Poole (1974), Magurran (1988). This is an alternative to the randomization test available in the Diversity permutation test module. Requires two columns of abundance data with taxa down the rows. The Shannon index here include a bias correction and may diverge slightly from the uncorrected estimates calculated elsewhere in PAST, at least for small samples. With pi the proportion (0-1) of taxon i, S the number of taxa and N the number of individuals, the estimator of the index is    S i ii N S ppH 1 2 1 ln' (note that the second term is incorrect in Magurran 1988). The variance of the estimator is      2 22 2 1lnln 'Var N S N pppp H iiii       . The t test statistic is given by 21 21 VarVar HH HH t    . The degrees of freedom for the t test is       2 2 2 1 2 1 2 21 VarVar VarVar N H N H HH df      . The Simpson index (dominance) has estimated variance (Brower et al. 1998):           22 2223 1 321212214 Var      NN pNNNpNNpNNN D iii . References Brower, J.E., Zar, J.H., von Ende, C.N. 1998. Field and Laboratory Methods for General Ecology. McGraw-Hill, Boston. Hutcheson, K. 1970. A test for comparing diversities based on the Shannon formula. Journal of Theoretical Biology 29:151-154. Magurran, A. 1988. Ecological Diversity and Its Measurement. Princeton University Press. Poole, R.W. 1974. An introduction to quantitative ecology. McGraw-Hill, New York. 157 Diversity profiles This module requires one or more columns of abundance data with taxa down the rows. The main purpose is to compare diversities in several samples. The validity of comparing diversities across samples can be criticized because of arbitrary choice of diversity index. One sample may for example contain a larger number of taxa, while the other has a larger Shannon index. A number of diversity indices may be compared to make sure that the diversity ordering is robust. A formal way of doing this is to define a family of diversity indices, dependent upon a single continuous parameter (Tothmeresz 1995). PAST uses the exponential of the so-called Renyi index, which depends upon a parameter . For =0, this function gives the total species number. =1 (in the limit) gives an index proportional to the Shannon index, while =2 gives an index which behaves like the Simpson index.            S i ipH 1 ln 1 1 expexp    The program can plot several such diversity profiles together. If the profiles cross, the diversities are non-comparable. The bootstrapping option (giving a 95% confidence interval) is based on 2000 replicates. Reference Tothmeresz, B. 1995. Comparison of different methods for diversity ordering. Journal of Vegetation Science 6:283-290. 158 Time series menu Simple periodogram Since palaeontological data are often unevenly sampled, Fourier-based methods can be difficult to use. PAST therefore includes the Lomb periodogram algorithm for unevenly sampled data (Press et al. 1992), with time values given in the first column and dependent values in the second column. If only one column is selected, an even spacing of one unit between data points is assumed. The Lomb periodogram should then give similar results as the FFT. The data are automatically detrended prior to analysis. The frequency axis is in units of 1/(x unit). If for example, your x values are given in millions of years, a frequency of 0.1 corresponds to a period of 10 million years. The power axis is in units proportional to the square of the amplitudes of the sinusoids present in the data. Also note that the frequency axis extends to very high values. If your data are evenly sampled, the upper half of the spectrum is a mirror image of the lower half, and is of little use. If some of your regions are closely sampled, the algorithm may be able to find useful information even above the half-point (Nyquist frequency). The highest peak in the spectrum is presented with its frequency and power value, together with a probability that the peak could occur from random data. The 0.01 and 0.05 significance levels ('white noise lines') are shown as red dashed lines. The example above shows a spectral analysis of a foram oxygen isotope record from 1 Ma to Recent, with an even spacing of 0.003 Ma (3 ka). There are periodicities at frequencies of about 9 (split peak), 25 and 43 Ma-1 , corresponding to periods of 111 ka, 40 ka and 23 ka – clearly orbital forcing. Reference Press, W.H., S.A. Teukolsky, W.T. Vetterling & B.P. Flannery. 1992. Numerical Recipes in C. Cambridge University Press. 159 REDFIT spectral analysis This module is an implementation of the REDFIT procedure of Schulz and Mudelsee (2002). It is a more advanced version of the simple Lomb periodogram described above. REDFIT includes an option for “Welch overlapped segment averaging”, which implies splitting the time series into a number of segments, overlapping by 50%, and averaging their spectra. This reduces noise but also reduces spectral resolution. In addition, the time series is fitted to an AR(1) red noise model which is usually a more appropriate null hypothesis than the white noise model described above. The given “falsealarm lines” are based on both parametric approximations (chi2) and Monte Carlo (using 1000 random realizations of an AR(1) process). The input must be in the form of two columns with time and data values, or one column of equallyspaced data values. The data are automatically detrended. The fitting to AR(1) implies that the data must have the correct time direction (in contrast with the simple spectrogram above where the time direction is arbitrary). The time values are expected to be ages before present. If not, it will be necessary to give them negative signs. The frequency oversampling value controls the number of points along the frequency axis (but having many points does not increase frequency resolution!). Increasing the number of segments will reduce noise, but also decrease the resolution. The window function influences the trade-off between spectral resolution and attenuation of side lobes. The (average) tau value is the characteristic time scale (the parameter of the AR model). The bandwidth is the spectral resolution given as the width between the -6dB points. The fit to an AR(1) model can be assessed using the runs value and its 5% acceptance interval. This test is only available with Monte Carlo on, oversampling=1, segments=1, window=rectangular. 160 In addition to a fixed set of false-alarm levels (80%, 90%, 95% and 99%), the program also reports a “critical” false-alarm level (False-al) that depends on the segment length (Thomson 1990). Important: Because of long computation time, the Monte Carlo simulation is not run by default, and the Monte Carlo false-alarm levels are therefore not available. When the Monte Carlo option is enabled, the given spectrum may change slightly because the Monte Carlo results are then used to compute a “bias-corrected” version (see Schulz and Mudelsee 2002). References Schulz, M. & M. Mudelsee. 2002. REDFIT: estimating red-noise spectra directly from unevenly spaced paleoclimatic time series. Computers & Geosciences 28:421-426. Thomson, D.J. 1990. Time series analysis of Holocene climate data. Philosophical Transactions of the Royal Society of London, Series A 330:601-616. 161 Multitaper spectral analysis In traditional spectral estimation, the data are often “windowed” (multiplied with a bell-shaped function) in order to reduce spectral leakage. In the multitaper method, several different (orthogonal) window functions are applied, and the results combined. The resulting spectrum has low leakage, low variance, and retains information contained in the beginning and end of the time series. In addition, statistical testing can take advantage of the multiple spectral estimates. One possible disadvantage is reduced spectral resolution. The multitaper method requires evenly spaced data, given in one column. The implementation in Past is based on the code of Lees & Park (1995). The multitaper spectrum can be compared with a simple periodogram (FFT with a 10% cosine window) and a smoothed periodogram. The number of tapers (NWIN) can be set to 3, 4 or 5, for different tradeoffs between variance reduction and resolution. The “time-bandwidth product” p is fixed at 3.0. The F test for significance of periodicity follows Lees & Park (1995). The 0.05 and 0.01 significance levels are shown as horizontal lines, based on 2 and 2*NWIN-2 degrees of freedom. The data are zero-padded to the second lowest power of 2 above the length of the sequence. This is required to reproduce the test results given by Lees & Park (1995). Reference Lees, J.M. & J. Park. 1995. Multiple-taper spectral analysis: a stand-alone C-subroutine. Computers & Geosciences 21:199-236. 162 Walsh transform The Walsh transform is a type of spectral analysis (for finding periodicities) of binary or ordinal data. It assumes even spacing of data points, and expects one column of binary (0/1) or ordinal (integer) data. The normal methods for spectral analysis are perhaps not optimal for binary data, because they decompose the time series into sinusoids rather than "square waves". The Walsh transform may then be a better choice, using basis functions that flip between -1 and +1. These basis functions have varying "frequencies" (number of transitions divided by two), known as sequencies. In PAST, each pair of even ("cal") and odd ("sal") basis functions is combined into a power value using cal2 +sal2 , producing a "power spectrum" that is comparable to the Lomb periodogram. In the example above, compare the Walsh periodogram (top) to the Lomb periodogram (bottom). The data set has 0.125 periods per sample. Both analyses show harmonics. The Walsh transform is slightly exotic compared with the Fourier transform, and the results must be interpretated cautiously. For example, the effects of the duty cycle (percentage of ones versus zeros) are somewhat difficult to understand. In PAST, the data values are pre-processed by multiplying with two and subtracting one, bringing 0/1 binary values into the -1/+1 range optimal for the Walsh transform. The data are zero-padded to the next power of 2 if necessary, as required by the method. 163 Short-time Fourier transform Spectral analysis using the Fourier transform (FFT), but dividing the signal into a sequence of overlapping windows, which are analysed individually. This allows development of the spectrum in time, in contrast with the global analysis provided by the other spectral analysis modules. Sample position is shown on the x axis, frequency (in periods per sample) on the y axis, and power on a logarithmic scale as colour or grey scale. The Short-time Fourier Transform (STFT) can be compared with wavelet analysis, but with a linear frequency scale and with constant time resolution independent of frequency. The window size controls the trade-off between resolution in time and frequency; small windows give good time resolution but poor frequency resolution. Windows are zero-padded by a factor eight to give a smoother appearance of the diagram along the frequency axis. The window functions (Rectangle, Welch, Hanning, Hamming, Blackman-Harris, multitaper with 3, 4 or 5 tapers) give different trade-offs between frequency resolution and sideband rejection. 164 Wavelet transform Inspection of time series at different scales. Requires one column of ordinal or continuous data with even spacing of points. The continuous wavelet transform (CWT) is an analysis method where a data set can be inspected at small, intermediate and large scales simultaneously. It can be useful for detecting periodicities at different wavelengths, self-similarity and other features. The vertical axis in the plot is a logarithmic size scale (base 2), with the signal observed at a scale of only two consecutive data points at the top, and at a scale of one fourth of the whole sequence at the bottom. One unit on this axis corresponds to a doubling of the size scale. The top of the figure thus represents a detailed, fine-grained view, while the bottom represents a smoothed overview of longer trends. Signal power (or more correctly squared correlation strength with the scaled mother wavelet) is shown with a grayscale or in colour. The shape of the mother wavelet can be set to Morlet (wavenumber 6), Paul (4th order) or DOG (Derivative Of Gaussian, 2nd or 6th derivative). The Morlet wavelet usually performs best. The example above is based on a foram oxygen isotope record from 1 Ma to Recent, with an even spacing of 0.003 Ma (3 ka). A band can be seen at a scale of about 25 =32 samples, or about 100 ka. A weaker band around 23.7 =13 samples corresponds to a scale of about 40 ka. These are orbital periodicities. In contrast with the “bulk” spectral analysis, the scalogram makes visible changes in strength and frequency over time. The so-called “cone of influence” can be plotted to show the region where boundary effects are present. The algorithm is based on fast convolution of the signal with the wavelet at different scales, using the FFT. 165 Significance test: The significance level corresponding to p=0.05 can be plotted as a contour (chisquared test according to Torrence & Compo 1998). The “Lag” value, as given by the user, specifies the null hypothesis. Lag=0 specifies a white-noise model. Values 030, and the option is not available for n>40. Missing data supported. Reference Davis, J.C. 1986. Statistics and Data Analysis in Geology. John Wiley & Sons. 169 Cross-correlation Cross-correlation (Davis 1986)is carried out on two column(s) of evenly sampled temporal/stratigraphic data. The x axis shows the displacement of the second column with respect to the first, the y axis the correlation between the two time series for a given displacement. The "p values" option will draw the significance of the correlation, after Davis (1986). For two time series x and y, the cross-correlation value at lag time m is               22 yyxx yyxx r mii mii m . The summations and the mean values are only taken over the parts where the sequences overlap for a given lag time. The equation shows that for positive lags, x is compared with a y that has been delayed by m samples. A high correlation value at positive lags thus means that features in y are leading, while x lags behind. For negative lags, features in x are leading. A reminder of this is given by the program. The p value for a given m is given by a t test with n-2 degrees of freedom, with n the number of samples that overlap: 2 1 2 m m r n rt    . It is important to note that this test concerns one particular m. Plotting p as a function of all m raises the issue of multiple testing – p values smaller than 0.05 are expected for 5% of lag times even for completely random (uncorrelated) data sets. In the example above, the “earthquakes” data seem to lag behind the “injection” data with a delay of 0-2 samples (months in this case), where the correlation values are highest. The p values (red curve) indicates significance at these lags. Curiously, there also seems to be significance for negative correlation at large positive and negative lags. Missing data supported. Reference Davis, J.C. 1986. Statistics and Data Analysis in Geology. John Wiley & Sons. 170 Mantel correlogram (and periodogram) This module expects several rows of multivariate data, one row for each sample. Samples are assumed to be evenly spaced in time. The Mantel correlogram (e.g. Legendre & Legendre 1998) is a multivariate extension to autocorrelation, based on any similarity or distance measure. The Mantel correlogram in PAST shows the average similarity between the time series and a time lagged copy, for different lags. The Mantel periodogram is a power spectrum of the multivariate time series, computed from the Mantel correlogram (Hammer 2007). The Mantel scalogram is an experimental plotting of similarities between all pairs of points along the time series. The apex of the triangle is the similarity between the first and last point. The base of the triangle shows similarities between pairs of consecutive points. 171 References Hammer, Ø. 2007. Spectral analysis of a Plio-Pleistocene multispecies time series using the Mantel periodogram. Palaeogeography, Palaeoclimatology, Palaeoecology 243:373-377. Legendre, P. & L. Legendre. 1998. Numerical Ecology, 2nd English ed. Elsevier, 853 pp. 172 Runs test The runs test is a non-parametric test for randomness in a sequence of values such as a time series. Non-randomness may include such effects as autocorrelation, trend and periodicity. The module requires one column of data, which are internally converted to 0 (x≤0) or 1 (x>0). The test is based on a dichotomy between two values (x≤0 or x>0). It counts the number of runs (groups of consecutive equal values) and compares this to a theoretical value. The runs test can therefore be used directly for sequences of binary data. There are also options for “runs about the mean” (the mean value subtracted from the data prior to testing), or “runs up and down” (the differences from one value to the next taken before testing). With n the total number of data points, n1 the number of points ≤0 and n2 the number of points >0, the expected number of runs in a random sequence, and the variance, are   n nnn RE 212  .      1 22 2 2121    nn nnnnn RVar . With the observed number of runs R, a z statistic can be written as   )(RVar RER z   . The resulting two-tailed p value is not accurate for n<20. A Monte Carlo procedure is therefore also included, based on 10,000 random replicates using the observed n, n1 and n2. 173 ARMA (and intervention analysis) NOT YET IN PAST 3 Analysis and removal of serial correlations in time series, and analysis of the impact of an external disturbance ("intervention") at a particular point in time. Stationary time series, except for a single intervention. One column of equally spaced data. This powerful but somewhat complicated module implements maximum-likelihood ARMA analysis, and a minimal version of Box-Jenkins intervention analysis (e.g. for investigating how a climate change might impact biodiversity). By default, a simple ARMA analysis without interventions is computed. The user selects the number of AR (autoregressive) and MA (moving-average) terms to include in the ARMA difference equation. The log-likelihood and Akaike information criterion are given. Select the numbers of terms that minimize the Akaike criterion, but be aware that AR terms are more "powerful" than MA terms. Two AR terms can model a periodicity, for example. The main aim of ARMA analysis is to remove serial correlations, which otherwise cause problems for model fitting and statistics. The residual should be inspected for signs of autocorrelation, e.g. by copying the residual from the numerical output window back to the spreadsheet and using the autocorrelation module. Note that for many paleontological data sets with sparse data and confounding effects, proper ARMA analysis (and therefore intervention analysis) will be impossible. The program is based on the likelihood algorithm of Melard (1984), combined with nonlinear multivariate optimization using simplex search. Intervention analysis proceeds as follows. First, carry out ARMA analysis on only the samples preceding the intervention, by typing the last pre-intervention sample number in the "last samp" box. It is also possible to run the ARMA analysis only on the samples following the intervention, by typing the first post-intervention sample in the "first samp" box, but this is not recommended because of the post-intervention disturbance. Also tick the "Intervention" box to see the optimized intervention model. The analysis follows Box and Tiao (1975) in assuming an "indicator function" u(i) that is either a unit step or a unit pulse, as selected by the user. The indicator function is transformed by an AR(1) process with a parameter delta, and then scaled by a magnitude (note that the magnitude given by PAST is the coefficient on the transformed indicator function: first do y(i)=delta*y(i-1)+u(i), then scale y by the magnitude). The algorithm is based on ARMA transformation of the complete sequence, then a corresponding ARMA transformation of y, and finally linear regression to find the magnitude. The parameter delta is optimized by exhaustive search over [0,1]. For small impacts in noisy data, delta may end up on a sub-optimum. Try both the step and pulse options, and see what gives smallest standard error on the magnitude. Also, inspect the "delta optimization" data, where standard error of the estimate is plotted as a function of delta, to see if the optimized value may be unstable. The Box-Jenkins model can model changes that are abrupt and permanent (step function with delta=0, or pulse with delta=1), abrupt and non-permanent (pulse with delta<1), or gradual and permanent (step with delta<0). 174 Be careful with the standard error on the magnitude - it will often be underestimated, especially if the ARMA model does not fit well. For this reason, a p value is deliberately not computed (Murtaugh 2002). The example data set (blue curve) is Sepkoski’s curve for percent extinction rate on genus level, interpolated to even spacing at ca. 5 million years. The largest peak is the Permian-Triassic boundary extinction. The user has specified an ARMA(2,0) model. The residual is plotted in red. The user has specified that the ARMA parameters should be computed for the points before the P-T extinction at time slot 37, and a pulse-type intervention. The analysis seems to indicate a large time constant (delta) for the intervention, with an effect lasting into the Jurassic. References Box, G.E.P. & G.C. Tiao. 1975. Intervention analysis with applications to economic and environental problems. Journal of the American Statistical Association 70:70-79. Melard, G. 1984. A fast algorithm for the exact likelihood of autoregressive-moving average models. Applied Statistics 33:104-114. Murtaugh, P.A. 2002. On rejection rates of paired intervention analysis. Ecology 83:1752-1761. 175 Insolation (solar forcing) model NOT YET IN PAST 3 This module computes solar insolation at any latitude and any time from 100 Ma to the Recent (the results are less accurate before 50 Ma). The calculation can be done for a "true" orbital longitude, "mean" orbital longitude (corresponding to a certain date in the year), averaged over a certain month in each year, or integrated over a whole year. The implementation in PAST is ported from the code by Laskar et al. (2004), by courtesy of these authors. Please reference Laskar et al. (2004) in any publications. It is necessary to specify a data file containing orbital parameters. Download the file INSOLN.LA2004.BTL.100.ASC from http://www.imcce.fr/Equipes/ASD/insola/earth/La2004 and put in anywhere on your computer. The first time you run the calculation, PAST will ask for the position of the file. The amount of data can become excessive for long time spans and short step sizes! Reference Laskar, J., P. Robutel, F. Joutel, M. Gastineau, A.C.M. Correia & B. Levrard. 2004. A long-term numerical solution for the insolation quantities of the Earth. Astronomy & Astrophysics 428:261-285. 176 Point events Expects one column containing times of events (e.g. earthquakes or clade divergences) or positions along a line (e.g. a transect). The times do not have to be in increasing order. Density trend (Laplace test) The “Laplace” test for a trend in density (intensity) is described by Cox & Lewis (1978). It is based on the test statistic n L L t U 12 1 2   where t is the mean event time, n the number of events and L the length of the interval. L is estimated as the time from the first to the last event, plus the mean waiting time. U is approximately normally distributed with zero mean and unit variance under the null hypothesis of constant intensity. This is the basis for the given p value. If p<0.05, a positive U indicates an increasing trend in intensity (decreasing waiting times), while a negative U indicates a decreasing trend. Note that if a trend is detected by this test, the sequence is not stationary and the assumptions of the exp test below are violated. Exp test for Poisson process The exp test (Prahl 1999) for a stationary Poisson process (random, independent events) is based on the set of n waiting times Δti between successive events in the sorted sequence . The test statistic is:          Tt i i T t n M 1 1 177 where T is the mean waiting time. M will tend to zero for a regularly spaced (overdispersed) sequence, and to 1 for a highly clustered sequence. For the null hypothesis of a Poisson process, M is asymptotically normally distributed with mean 1/e - /n and standard deviation β/n, where =0.189 and β=0.2427. This is the basis for the given z test. In summary, if p<0.05 the sequence is not Poisson. You can then inspect the M statistic; if smaller than the expected value this indicates regularity, if higher it indicates clustering. For both tests, p values are also estimated by Monte Carlo simulation with 9999 random data sets. References Cox, D. R. & P. A. W. Lewis. 1978. The Statistical Analysis of Series of Events. Chapman and Hall, London. Prahl, J. 1999. A fast unbinned test on event clustering in Poisson processes. Arxiv, Astronomy and Astrophysics September 1999. 178 Markov chain This module requires a single column containing a sequence of nominal data coded as integer numbers. For example, a stratigraphic sequence where 1 means limestone, 2 means shale and 3 means sand. A transition matrix containing counts or proportions (probabilities) of state transitions is displayed. The “from”-states are in rows, the “to”-states in columns. It is also possible to specify several columns, each containing one or more state transitions (two numbers for one transition, n numbers for a sequence giving n-1 transitions). The chi-squared test reports the probability that the data were taken from a system with random proportions of transitions (i.e. no preferred transitions). The transitions with anomalous frequencies can be identified by comparing the observed and expected transition matrices. The “Embedded (no repeats)” option should be selected if the data have been collected in such a way that no transitions to the same state are possible (data points are only collected when there is a change). The transition matrix will then have zeroes on the diagonal. The algorithms, including an iterative algorithm for embedded Markov chains, are according to Davis (1986). Reference Davis, J.C. 1986. Statistics and Data Analysis in Geology. John Wiley & Sons. 179 Simple smoothers A set of simple smoothers for a single column of evenly spaced data. Missing data are supported. Moving average Simple n-point, centered moving average (n must be odd). Commonly used, but has unfortunate properties such as a non-monotonic frequency response. Gaussian Weighted moving average using a Gaussian kernel with standard deviation set to 1/5 of the window size (of n points). This is probably the best overall method in the module. Moving median Similar to moving average but takes the median instead of the mean. This method is more robust to outliers. AR 1 (exponential) Recursive (autoregressive) filter, yi = yi-1 + (1-)xi with  a smoothing coefficient from 0 to 1. This corresponds to weighted averaging with exponentially decaying weights. Gives a phase delay and also a transient in the beginning of the series. Included for completeness. 180 FIR filter Filtering out certain frequency bands in a time series can be useful to smooth a curve, remove slow variation, or emphasize certain periodicities (e.g. Milankovitch cycles). One column of evenly spaced data is expected. For most applications in data analysis, it is crucial that the filter has linear phase response. Past therefore uses FIR (Finite Impulse Response) filters, which are designed using the Parks-McClellan algorithm. The following filter types are available: Lowpass, highpass, bandpass and bandstop. Filter parameters To design an optimal filter takes a little effort. Frequencies are specified in the range 0-0.5, i.e. T0/T where T0 is the sampling interval (not specified to the computer) and T is the required period. For example, if your real sampling interval is 1,000 years, a frequency corresponding to a period of 23,000 years is specified as 1,000/23,000=0.043. After setting the filter type, you should select a transition width (or leave the default of 0.02). Decreasing the transition width will make a sharper filter, at the cost of larger ripple (“waves” in the frequency response). Note that the values in text fields are not updated until you press Enter. Also, if an invalid combination is entered (e.g. a transition band crossing 0 or 0.5, or upper limit less than lower limit) the program will reset some value to avoid errors. It is therefore required to enter the numbers in an order so that the filter is always valid. 181 The filter types are as follows: 1. Lowpass. The From frequency is forced to zero. Frequencies up to the To frequency pass the filter. Frequencies from To+Transition to 0.5 are blocked. 2. Highpass. The To frequency is forced to 0.5. Frequencies above the From frequency pass the filter. Frequencies from 0 to From-Transition are blocked. 3. Bandpass. Frequencies from From to To pass the filter. Frequencies below From-Transition and above To+Transition are blocked. 4. Bandstop. Frequencies from From to To are blocked. Frequencies from 0 to From-Transition and from To+Transition to 0.5 pass the filter. Filter order The filter order should be large enough to give an acceptably sharp filter with low ripple. However, a filter of length n will give less accurate results in the first and last n/2 samples of the the time series, which puts a practical limit on filter order for short series. The Parks-McClellan algorithm will not always converge. This gives an obviously incorrect frequency response, and attempting to apply such a filter to the data will give a warning message. Try to change the filter order (usually increase it) to fix the problem. 182 Date/time conversion Utility to convert dates and/or times to a continuous time unit for analysis. The program expects one or two columns, each containing dates or times. If both are given, then time is added to date to give the final time value. Dates can be given in the formats Year/Month/Day or Day/Month/Year. Years need all digits (a year given as 11 will mean 11 AD, not 2011). Only Gregorian calendar dates are supported. Leap years are taken into account. Time can be given as Hours:Minutes or Hours:Minutes:Seconds (seconds can include decimals). The output units can be years (using the Gregorian mean year of 365.2425 days), days (of 86400 seconds), hours, minutes or seconds. The starting time (time zero) can be the smallest given time, the beginning of the first day, the beginning of the first year, year 0 (note the “astronomical” convention where the year before year 1 is year 0), or the beginning of the first Julian day (noon, year -4712). The program operates with simple (UT) time, defined with respect to the Earth’s rotation and with a fixed number of seconds (86400) per day. If your input data consists of space-separated date-time values, such as “2011/12/24 18:00:00.00”, then you may have to use the “Import text file” function to read the data such that dates and times are split into separate columns. The calculation of Julian day (which is used to find number of days between two dates) follows Meeus (1991): if month <= 2 begin year := year - 1; month := month + 12; end; A = floor(year/100); B = 2 – A + floor(A/4); JD = floor(365.25(year + 4716)) + floor(30.6001(month+1)) + day + B – 1524.5; Reference Meeus, J. 1991. Astronomical algorithms. Willmann-Bell, Richmond. 183 Geometrical menu Directions (one sample) The module plots a rose diagram (polar histogram) of directions. Used for plotting current-oriented specimens, orientations of trackways, fault lines, etc. Also appropriate for time-of day data (0-24 hours). One column of directional (0-360) or orientational (0-180) data in degrees is expected. Directional or periodic data in other forms (radians, hours, etc.) must be converted to degrees using e.g. the Evaluate Expression module (Transform menu). By default, the 'mathematical' angle convention of anticlockwise from east is chosen. If you use the 'geographical' convention of clockwise from north, tick the box. You can also choose whether to have the abundances proportional to radius in the rose diagram, or proportional to area (equal area). The "Kernel density" option plots a circular kernel density estimate. Descriptive statistics The mean angle takes circularity into account: i i       cos sin tan 1 (taken to the correct quadrant). 184 The 95 percent confidence interval on the mean is estimated according to Fisher (1983). It assumes circular normal distribution, and is not accurate for very large variances (confidence interval larger than 45 degrees) or small sample sizes. The bootstrapped 95% confidence interval on the mean uses 5000 bootstrap replicates. The graphic uses the bootstrapped confidence interval. The concentration parameter κ is estimated by iterative approximation to the solution to the equation     R I I    0 1 where I0 and I1 are imaginary Bessel functions of orders 0 and 1, estimated according to Press et al. (1992), and R defined below (see e.g. Mardia 1972). Rayleigh’s test for uniform distribution The R value (mean resultant length) is given by: nR n i i n i i 2 1 2 1 sincos                 . R is further tested against a random distribution using Rayleigh's test for directional data (Davis 1986). Note that this procedure assumes evenly or unimodally (von Mises) distributed data - the test is not appropriate for e.g. bimodal data. The p values are computed using an approximation given by Mardia (1972): 2 RnK             2 4322 288 97613224 4 2 1 n KKKK n KK ep K Rao’s spacing test for uniform distribution The Rao's spacing test (Batschelet 1981) for uniform distribution has test statistic   n i iTU 12 1  , where no 360 . iiiT   1 for i=20. Axial data The 'Orientations' option allows analysis of linear (axial) orientations (0-180 degrees). The Rayleigh and Watson tests are then carried out on doubled angles (this trick is described by Davis 1986); the Chi-square uses four bins from 0-180 degrees; the rose diagram mirrors the histogram around the origin. References Batschelet, E. 1981. Circular statistics in biology. Academic Press. Davis, J.C. 1986. Statistics and Data Analysis in Geology. John Wiley & Sons. Fisher, N.I. 1983. Comment on "A Method for Estimating the Standard Deviation of Wind Directions". Journal of Applied Meteorology 22:1971. Lockhart, R.A. & M.A. Stephens 1985. Tests of fit for the von Mises distribution. Biometrika 72:647- 652. Mardia, K.V. 1972. Statistics of directional data. Academic Press, London. Russell, G. S. & D.J. Levitin 1995. An expanded table of probability values for Rao's spacing test. Communications in Statistics: Simulation and Computation 24:879-888. 186 Directions (two samples) The module expects two columns of directional (0-360) or orientational (0-180) data in degrees. Watson-Williams test The Watson-Williams test for equal mean angle in two samples is a parametric test, assuming von Mises distribution, but is fairly robust. The concentration parameter κ should be larger than 1.0 for accurate testing. In addition, the test assumes similar angular variances (R values). The two samples φ and θ have n1 and n2 values. Rayleigh’s spread R is calculated for each sample and for the combined sample: 2 1 2 1 1 11 sincos                n i i n i iR  2 1 2 1 2 22 sincos                n i i n i iR  2 11 2 11 2121 sinsincoscos                n i i n i i n i i n i iR  . 187 The test statistic U is computed as    21 21 2 RRn RRR nU    . The significance is computed by first correcting U according to Mardia (1972a):                  95.0 8 3 1 45.0 1 8 1 2 2 nRU nR n U U    , where n=n1+n2. The p value is then given by the F distribution with 1 and n-2 degrees of freedom. The combined concentration parameter κ is maximum-likelihood, computed as described under “Directions (one sample)” above. Mardia-Watson-Wheeler test This non-parametric test for equal distribution is computed according to Mardia (1972b).           2 2 2 2 2 1 2 1 2 1 2 n SC n SC W where, for the first sample,    1 1 11 2cos n i i NrC  ,    1 1 11 2sin n i i NrS  and similarly for the second sample (N=n1+n2). The r1i are the ranks of the values of the first sample within the pooled sample. For N>14, W is approximately chi-squared with 2 degrees of freedom. References Mardia, K.V. 1972a. Statistics of directional data. Academic Press, London. Mardia, K.V. 1972b. A multi-sample uniform scores test on a circle and its parametric competitor. Journal of the Royal Statistical Society Series B 34:102-113. 188 Circular correlation Testing for correlation between two directional or orientational variates. Assumes “large” number of observations. Requires two columns of directional (0-360) or orientational (0-180) data in degrees. This module uses the circular correlation procedure and parametric significance test of Jammalamadaka & Sengupta (2001). The circular correlation coefficient r between vectors of angles α and β is               n i ii n i ii r 1 22 1 sinsin sinsin   , where the angular means are calculated as described previously. The test statistic T is computed as                 n k kk n k n k kk rT 1 22 1 1 22 sinsin sinsin   . For large n, this statistic has asymptotically normal distribution with mean 0 and variance 1 under the null hypothesis of zero correlation, which is the basis for the calculation of p. Reference Jammalamadaka, S.R. & A. Sengupta. 2001. Topics in circular statistics. World Scientific. 189 Spherical (one sample) This module makes stereo plots of axial, spherical data (e.g. strike-dip measurements in structural geology). Spherical statistics may be added in future versions. Three data formats can be used, all using the geographic angle convention (degrees, clockwise from north):  Trend (azimuth) and plunge (angle down from the horizontal) for axial data  Dip azimuth and dip angle (down from the horizontal) for planes. The pole (normal vector) of the plane is plotted.  Strike and dip for planes, using the right-hand rule convention with the dip down to the right from the strike. The pole to the plane is plotted. Density contouring is based on a modified Kamb method algorithm by Vollmer (1995). Both equal area (Schmidt) and equal angle (Wulff) projections are available. Projections are to the lower hemisphere. Density estimates can use an inverse area, inverse area squared or exponential law, giving progressively higher smoothing. Reference Vollmer, F.W. 1995. C program for automatic contouring of spherical orientation data using a modified Kamb method. Computers & Geosciences 21:31-49. 190 Point pattern analysis - nearest neighbours This module tests for clustering or overdispersion of points given as two-dimensional coordinate values. The procedure assumes that elements are small compared to their distances, that the domain is predominantly convex, and n>50. Two columns of x/y positions are required. Applications of this module include spatial ecology (are in-situ brachiopods clustered), morphology (are trilobite tubercles overdispersed), and geology (distribution of e.g. volcanoes, earthquakes, springs). The calculation of point distribution statistics using nearest neighbour analysis follows Davis (1986) with modifications. The area is estimated either by the smallest enclosing rectangle or using the convex hull, which is the smallest convex polygon enclosing the points. Both are inappropriate for points in very concave domains. Two different edge effect adjustment methods are available: wraparound ("torus") and Donnelly's correction. Wrap-around edge detection is only appropriate for rectangular domains. The null hypothesis is a random Poisson process, giving a modified exponential nearest neighbour distribution (see below) with mean 2 nA  where A is the area and n the number of points. The probability that the distribution is Poisson is presented, together with the R value: 191 nA dd R 2   where d is the observed mean distance between nearest neighbours. Clustered points give R<1, Poisson patterns give R~1, while overdispersed points give R>1. The expected (theoretical) distribution under the null hypothesis is plotted as a continuous curve together with the histogram of observed distances. The expected probability density function as a function of distance r is    2 πexpπ2 rrrg   where  = n/A is the point density (Clark & Evans 1954). The orientations (0-180 degrees) and lengths of lines between nearest neighbours, are also included. The orientations can be subjected to directional analysis to test whether the points are organised along lineaments (see Hammer 2009 for more advanced methods). References Clark, P.J. & Evans, F.C. 1954. Distance to nearest neighbor as a measure of spatial relationships in populations. Ecology 35:445-453. Davis, J.C. 1986. Statistics and Data Analysis in Geology. John Wiley & Sons. Hammer, Ø. 2009. New methods for the statistical analysis of point alignments. Computers & Geosciences 35:659-666. 192 Ripley’s K point pattern analysis Ripley's K (Ripley 1979) is the average point density as a function of distance from every point. It is useful when point pattern characteristics change with scale, e.g. overdispersion over small distances but clustering over large distances. Two columns of x/y coordinates in a rectangular domain are expected. Define the estimated intensity of the point pattern, with n points in an area A, as An . The distance between points i and j is dij The estimate of Ripley’s K, as a function of distance, is then computed as       n i ij ij ddI n dK 1 1  , where the indicator function I is one if the argument is true, zero otherwise. The normalization of K is such that for complete spatial randomness (CSR), K(d) is expected to increase as the area of circles, i.e.   2 ddK  . The L(d) function is a corresponding transformation of K(d):      dK dL  For CSR, L(d)=d, and L(d)-d=0. A 95% confidence interval for CSR is estimated using 1000 Monte Carlo simulations within the bounding rectangle (previous versions used the approximation nA42.1 ). 193 Ripley's edge correction is included, giving weights to counts depending on the proportion of the test circle that is inside the rectangular domain. The example above shows locations of volcanic pipes. L(d)-d is below the 95% confidence interval of CSR, indicating lateral inhibition, up to a distance of ca. 70 m. For larger distances, the curve flattens in the manner expected from CSR. Area For the correct calculation of Ripley's K, the area must be known. In the first run, the area is computed using the smallest bounding rectangle, but this can both over- and underestimate the real area. The area can therefore be adjusted by the user. An overestimated area will typically show up as a strong overall linear trend with positive slope for L(d)-d. Fractal dimension The fractal dimension (if any) can be estimated as the asymptotic linear slope in a log-log plot of K(d). For CSR, the log-log slope should be 2.0. Fractals should have slopes less than 2. References Ripley, B.D. 1979. Tests of 'randomness' for spatial point patterns. Journal of the Royal Statistical Society, ser. B 41:368-374. 194 Kernel density Makes a smooth map of point density in 2D. Two columns of x/y coordinates in a rectangular domain are expected. The user can specify the size of the grid (number of rows and columns). The “Radius” value sets the scale r of the kernel. There is currently no automatic selection of “optimal” radius, so this value must be set by the user depending on the scale of interest. The density estimate is based on one of four kernel functions, with radius parameter r. With    22 iii yyxxd  : Gaussian (default):             i i r d r yx 2 2 2 2 exp 1 ,f  Paraboloid:          i i i i rd rd r d r yx 0 1 2 3 ,f 2 2 2  Triangular:           i i i i rd rd r d r yx 0 12 ,f 2  Uniform:          i i i rd rd r yx 0 11 ,f 2  195 The scaling gives an estimate of the number of points per area, not a probability density. The gaussian and paraboloid (quadratic) kernels usually perform best. The uniform kernel gives very low smoothness. 196 Point alignments Detection of linear alignments in a 2D point pattern, using the continuous sector method (Hammer 2009). Typical applications are in geology and geography, to study the distribution of earthquakes, volcanoes, springs etc. associated with faults and other linear structures. The Radius parameter sets the scale of analysis. In the example above, lineaments of length 1200 m (twice the radius) are detected. Alpha sets the significance level for the Rayleigh test used by the procedure. Note that this is a pointwise significance, not corrected for the multiple testing of all the points. The Dispersion filter disables alignments with uneven distribution of points along the lineament. View numbers lists the alignment positions and their orientations, which can be subjected to circular statistics if required (Directions module). Reference Hammer, Ø. 2009. New methods for the statistical detection of point alignments. Computers & Geosciences 35:659-666. 197 Spatial autocorrelation (Moran’s I) Spatial autocorrelation in Past requires three columns, containing x and y coordinates and corresponding data values z for a number of points. The Moran’s I correlation statistic is then computed within each of a number of distance classes (bins), ranging from small to large distances. The one-tailed critical value for p<0.05 can be plotted for each bin. Moran’s I values exceeding the critical value may be considered significant, but Bonferroni or other adjustment for multiple testing should be considered because of the several bins. The calculation follows Legendre & Legendre (1998). For each distance class d, compute               n i i n h n i ihhi zz n zzzzw W dI 1 2 1 1 1 1 . Here, n is the total number of points, W is the number of pairs of points having distances within the distance class, and whi a weight function such that whi=1 if points h and i are within the distance class and whi=0 otherwise (Kronecker delta). Note that this equation is incorrect in some publications. 198 For the one-tailed critical level I0.05, compute                           1 05.005.0 22 2 21 2 2 2 21 2 2 1 2 1 4 2 1 2 2 1 1 2 1 1Var6452.1 1 1 321 62333 Var 2 1                             nkII nWnnn WnSSnnbWnSSnnn I zz zzn b wwS wwS n i i n i i n i ii n h n i ihhi Here the wi+ and w+i are the row and column sums. The correction factor k0.05 is set to 707.005.010  if    13244  nnWnn , otherwise k0.05=1. Reference Legendre, P. & Legendre, L. 1998. Numerical Ecology, 2nd English ed. Elsevier, 853 pp. 199 Gridding (spatial interpolation) “Gridding” is the operation of spatial interpolation of scattered 2D data points onto a regular grid. Three columns with position (x,y) and corresponding data values are required. Gridding allows the production of a map showing a continuous spatial estimate of some variate such as fossil abundance or thickness of a rock unit, based on scattered data points. The user can specify the size of the grid (number of rows and columns). The spatial coverage of the map is generated automatically as a square covering the data points. When plotting, this can be reduced to the convex hull of the points. A least-squares linear surface (trend) is automatically fitted to the data, removed prior to gridding and finally added back in. This is primarily useful for the semivariogram modelling and the kriging method. Cross validation: This option will remove each data point in turn and re-compute the surface based on the remaining points (“jackknife”). The differences between the original data values and the crossvalidated values indicate the prediction accuracy of the surface model. These differences are reported for each point, together with the mean squared error (MSE) over all points. Four interpolation algorithms are available: Inverse distance weighting The value at a grid node is simply the average of the N closest data points, as specified by the user (the default is to use all data points). The points are weighted in inverse proportion to distance. This algorithm is fast but will not always give good (smooth) results. A typical artefact is “bull’s eyes” around data points. One advantage is that the interpolated values will never exceed the range of the 200 data points. By setting N=1, this algorithm reduces to the nearest-neighbour method, which sets the value at a grid node to the value of the nearest data point. Thin-plate spline Maximally smooth interpolator. Can overshoot in the presence of sharp bends in the surface. This is a radial basis method with radial basis function φ = r ln r. Multiquadric Radial basis function φ = r. Popular for terrain modelling. Kriging The user is required to specify a model for the semivariogram, by choosing one of four common models and corresponding parameters to fit the empirical semivariances (the residual sum of squares should be as small as possible). The semivariogram is computed within each of a number of bins. Using the histogram option, choose a number of bins so that each bin (except possibly the rightmost ones) contains at least 30 distances. The nugget parameter is a constant added to the model. It implies a non-zero variance at zero distance, and will therefore allow the surface to not pass exactly through the given data points. The range controls the extent of the curve along the distance axis. In the equations below, the normalized distance value h represents distance/range. The scale controls the extent of the curve along the variance axis. Spherical:                1 1 2 1 2 3 3 hscalenugget hh h scalenugget h Exponential:    h escalenuggeth   1 Gaussian:    2 1 h escalenuggeth   Cubic:           1 175.05.375.87 7532 hscalenugget hhhhhscalenugget h The “Optimize all” button will select the model and parameters giving the smallest residual sum of squares in the semivariogram. This may not be what you want: For example you may wish to use a specific model or to have zero nugget in order to ensure exact interpolation. This will require setting the values manually. The kriging procedure also provides an estimate of standard errors across the map (this depends on the semivariogram model being accurate). Kriging in PAST does not provide for anisotropic semivariance. Warning: Kriging is slow, do not attempt it for more than ca. 1000 data points on a 100x100 grid. 201 See e.g. Davis (1986) or de Smith et al. (2009) for more information on gridding. References Davis, J.C. 1986. Statistics and Data Analysis in Geology. John Wiley & Sons. de Smith, M.J., M.F. Goodchild & P.A. Longley. 2009. Geospatial Analysis, 3rd ed. Matador. 202 Multivariate allometry NOT YET IN PAST 3 This module is used for investigating allometry in a multivariate morphometric data set. It expects a multivariate data set with variables (distance measurements) in columns, specimens in rows. This method for investigating allometry in a multivariate data set is based on Jolicoeur (1963) with extensions by Kowalewski et al. (1997). The data are (automatically) log-transformed and subjected to PCA. The first principal component (PC1) is then regarded as a size axis (this is only valid if the variation accounted for by PC1 is large, say more than 80%). The allometric coefficient for each original variable is estimated by dividing the PC1 loading for that variable by the mean PC1 loading over all variables. 95% confidence intervals for the allometric coefficients are estimated by bootstrapping specimens. 2000 bootstrap replicates are made. Missing data is supported by column average substitution. References Jolicoeur, P. 1963. The multivariate generalization of the allometry equation. Biometrics 19:497-499. Kowalewski, M., E. Dyreson, J.D. Marcot, J.A. Vargas, K.W. Flessa & D.P. Hallmann. 1997. Phenetic discrimination of biometric simpletons: paleobiological implications of morphospecies in the lingulide brachiopod Glottidia. Paleobiology 23:444-469. 203 PCA of 2D landmarks (relative warps) This module is very similar to the standard PCA module, but with some added functionality for analyzing 2D landmark configurations. The expected data are specimens in rows, alternating x and y coordinates in columns. Procrustes standardization recommended. The relative warps (principal components) are ordered according to importance, and the first and second warps are usually the most informative. Note that this module does a straightforward PCA of the landmarks, meaning that the affine component is included in the analysis. The relative warps are visualized with vectors or thin-plate spline transformation grids. When you increase or decrease the score factor away from zero, the original landmark configuration and grid will be progressively deformed according to the selected relative warp. 204 Thin-plate splines for 2D landmarks This module shows a shape deformation from one landmark configuration to another. The expected data are specimens in rows, alternating x and y coordinates in columns. Procrustes standardization recommended. Any shape selected in the “From shape” menu, is taken as a reference, with an associated square grid. The warps from this to all other specimens can be viewed. You can also choose the mean shape as the reference. The 'Expansion factors' option will display the area expansion (or contraction) factor around each landmark in yellow numbers, indicating the degree of local growth. This is computed using the Jacobian of the warp. Also, the expansions are colour-coded for all grid elements, with green for expansion and purple for contraction. At each landmark, the principal strains can also be shown, with the major strain in black and minor strain in brown. These vectors indicate directional stretching. A description of thin-plate spline transformation grids is given by Dryden & Mardia (1998). Reference Dryden, I.L. & K.V. Mardia 1998. Statistical Shape Analysis. Wiley. 205 Size from landmarks (2D or 3D) NOT YET IN PAST 3 Digitized x/y or x/y/z landmark coordinates. Specimens in rows, coordinates with alternating x and y (and z for 3D) values in columns. Must not be Procrustes fitted or normalized for size! Calculates the centroid size for each specimen (Euclidean norm of the distances from all landmarks to the centroid). The values in the 'Normalized' column are centroid sizes divided by the square root of the number of landmarks - this might be useful for comparing specimens with different numbers of landmarks. Normalize size The 'Normalize size' option in the Transform menu allows you to remove size by dividing all coordinate values by the centroid size for each specimen. For 2D data you may instead use Procrustes coordinates, which are also normalized with respect to size. See Dryden & Mardia (1998), p. 23-26. Reference Dryden, I.L. & K.V. Mardia 1998. Statistical Shape Analysis. Wiley. 206 Distance from landmarks (2D or 3D) NOT YET IN PAST 3 Digitized x/y or x/y/z landmark coordinates. Specimens in rows, coordinates with alternating x and y (and z for 3D) values in columns. May or may not be Procrustes fitted or normalized for size. Calculates the Euclidean distances between two fixed landmarks for one or many specimens. You must choose two landmarks - these are named according to the name of the first column for the landmark (x value). 207 All distances from landmarks (EDMA) NOT YET IN PAST 3 Digitized x/y or x/y/z landmark coordinates. Specimens in rows, coordinates with alternating x and y (and z for 3D) values in columns. May or may not be Procrustes fitted or normalized for size. This function will replace the landmark data in the data matrix with a data set consisting of distances between all pairs of landmarks, with one specimen per row. The number of pairs is N(N-1)/2 for N landmarks. This transformation will allow multivariate analysis of distance data, which are not sensitive to rotation or translation of the original specimens, so a Procrustes fitting is not mandatory before such analysis. Using distance data also allows log-transformation, and analysis of fit to the allometric equation for pairs of distances. Missing data is supported by column average substitution. 208 Landmark linking NOT YET IN PAST 3 This function in the Geomet menu allows the selection of any pairs of landmarks to be linked with lines in the morphometric plots (thin-plate splines, partial and relative warps, etc.), to improve readability. The landmarks must be present in the main spreadsheet before links can be defined. Pairs of landmarks are selected or deselected by clicking in the symmetric matrix. The set of links can also be saved in a text file. Note that there is little error checking in this module. 209 Fourier shape (2D) NOT YET IN PAST 3 Analysis of fossil outline shape (2D). Shape expressible in polar coordinates, sufficient number of digitized points to capture features. Digitized x/y coordinates around an outline. Specimens in rows, coordinates of alternating x and y values in columns (see Procrustes fitting in the Transform menu). Accepts X-Y coordinates digitized around an outline. More than one shape (row) can be simultaneously analyzed. Points do not need to be totally evenly spaced. The shape must be expressible as a unique function in polar co-ordinates, that is, any straight line radiating from the centre of the shape must cross the outline only once. The algorithm follows Davis (1986). The origin for the polar coordinate system is found by numerical approximation to the centroid. 128 points are then produced at equal angular increments around the outline, through linear interpolation. The centroid is then re-computed, and the radii normalized (size is thus removed from the analysis). The cosine and sine components are given for the first twenty harmonics, but note that only N/2 harmonics are 'valid', where N is the number of digitized points. The coefficients can be copied to the main spreadsheet for further analysis (e.g. by PCA). The 'Shape view' window allows graphical viewing of the Fourier shape approximation(s). Reference Davis, J.C. 1986. Statistics and Data Analysis in Geology. John Wiley & Sons. 210 Elliptic Fourier shape analysis Requires digitized x/y coordinates around outlines. Specimens in rows, coordinates of alternating x and y values in columns. Elliptic Fourier shape analysis is in several respects superior to simple Fourier shape analysis. One advantage is that the algorithm can handle complicated shapes which may not be expressible as a unique function in polar co-ordinates. Elliptic Fourier shapes is now a standard method of outline analysis. The algorithm used in PAST is described by Ferson et al. (1985). EFA coefficients Cosine and sine components of x and y increments along the outline for the first 30 harmonics are given, but only the first N/2 harmonics should be used, where N is the number of digitized points. Size and positional translation are normalized away, and do not enter in the coefficients. The size (before normalization) is given in the first column. The optional standardization for rotation or starting point, following Ferson et al., sometimes flips shapes around . This should be checked with the ‘Shape view’ (see below) – it may be necessary to remove such specimens. The coefficients can be copied to the main spreadsheet for further analysis such as discriminant analysis. The 'Shape view' window allows graphical viewing of the elliptic Fourier shape approximation(s). EFA PCA Principal Components Analysis of the EFA coefficients of the given outlines, with visualization of the principal components as EFA deformations. For more details on PCA in Past, see the description of PCA. Reference Ferson, S.F., F.J. Rohlf & R.K. Koehn. 1985. Measuring shape variation of two-dimensional outlines. Systematic Zoology 34:59-68. 211 Hangle Fourier shape analysis Requires digitized x/y coordinates around outlines. Specimens in rows, coordinates of alternating x and y values in columns. The “Hangle” method for analysing closed outlines, proposed by Haines & Crampton (2000) is a competitor to Elliptic Fourier Analysis. Hangle has certain advantages over EFA, the most important being that fewer coefficients are needed to capture the outline to a given precision. This is of importance for statistical testing (e.g. MANOVA) and discriminant analysis. The implementation in Past is based on the Hangle/Hmatch/Htree/Hshape package of Haines & Crampton (thanks to the authors for providing the source code). The output consists of 46 Fourier coefficients, which are the cos and sin coefficients of the first 24 harmonics (modes), starting on harmonic number 2. Copy these numbers back to a Past spreadsheet for further multivariate shape analysis. Starting point normalization Usually leave at ‘Match all’, either with the ‘Hmatch’ or (perhaps preferably) the ‘Htree’ method to align all the outlines. Alternatively, select 2.-4. harmonic, which will phase shift each outline according to the selected mode (see Haines & Crampton 2000). Smoothing Increasing the smoothing parameter can reduce high-frequency noise, at the cost of dampening potentially informative high-frequency shape information. Shape view Use this function to inspect the shapes reconstructed from the Fourier coefficients. Check that the matching routine has not rotated any shape incorrectly. Also, use this function to select the minimum number of modes necessary for capturing the shape. In the example above, the number of modes has been set to 14, which captures 99.88% of the total integrated power (amplitude squared) of the selected shape. The number of modes is shown by the red line in the power spectrum – make sure that the main features of the spectrum are to the left of this line for all the shapes. Note: PCA visualization and regression (as for EFA) has not yet been implemented for Hangle. Reference Haines, A.J. & J.S. Crampton. 2000. Improvements to the method of Fourier shape analysis as applied in morphometric studies. Palaeontology 43:765-783 212 Eigenshape analysis NOT YET IN PAST 3 Digitized x/y coordinates around an outline. Specimens in rows, coordinates of alternating x and y values in columns (see Procrustes fitting in the Transform menu). Eigenshapes are principal components of outlines. The scatter plot of outlines in principal component space can be shown, and linear combinations of the eigenshapes themselves can be visualized. The implementation in PAST is partly based on MacLeod (1999). It finds the optimal number of equally spaced points around the outline using an iterative search, so the original points need not be equally spaced. The eigenanalysis is based on the covariance matrix of the non-normalized turning angle increments around the outlines. The algorithm does not assume a closed curve, and the endpoints are therefore not constrained to coincide in the reconstructed shapes. Landmarkregistered eigenshape analysis is not included. All outlines must start at the 'same' point. Reference MacLeod, N. 1999. Generalizing and extending the eigenshape method of shape space visualization and analysis. Paleobiology 25:107-138. 213 Coordinate transformation Conversion between geographical coordinates in different grids and datums. The number of input columns depends on the data type, as described below. Decimal degrees (WGS84) Two columns: Latitude and longitude, in decimal degrees (60.5 is 60 degrees, 30 minutes). Negative values for south of equator and west of Greenwich. Referenced to the WGS84 datum. Deg/decimal mins (WGS84) Four columns: Latitude degrees, deci mal minutes (40.5 is 40 minutes, 30 seconds), longitude degrees, decimal minutes. Referenced to the WGS84 datum. Deg/min/sec (WGS84) Six columns: Latitude degrees, minutes, seconds, longitude degrees, minutes, seconds. Referenced to the WGS84 datum. UTM-ED50 (Intl 1924) Three columns: Easting (meters), northing (meters), and zone. Use negative zone numbers for the southern hemisphere. The handling of UTM zones takes into account the special cases of Svalbard and western Norway. Referenced to the ED50 European datum at Potsdam. UTM-WGS84 (WGS84) Three columns: Easting (meters), northing (meters), and zone. Referenced to the WGS84 datum. UTM-NAD27 (Clarke 1866) Three columns: Easting (meters), northing (meters), and zone. Referenced to the NAD27 datum. Conversion to/from this format is slightly inaccurate (5-6 meters). UTM-NAD83 (GRS80) Three columns: Easting (meters), northing (meters), and zone. Referenced to the NAD83 datum (practically identical to WGS84). Sweden (RT90) Two columns: Easting (meters) and northing (meters). The transformations are based on code generously provided by I. Scollar. 214 Stratigraphy menu Unitary Associations Unitary Associations analysis (Guex 1991) is a method for biostratigraphical correlation (see Angiolini & Bucher 1999 for a typical application). The data input consists of a presence/absence matrix with samples in rows and taxa in columns. Samples belonging to the same section (locality) must be assigned to the same group, and ordered stratigraphically within each section such that the lowermost sample enters in the lowest row. Overview of the method The method of Unitary Associations is logical, but rather complicated, consisting of a number of steps. For details, see Guex (1991). The implementation in PAST includes most of the features found in the original program, called BioGraph (Savary & Guex 1999), and thanks to a fruitful co-operation with Jean Guex it also includes a number of additional options and improvements. The basic idea is to generate a number of assemblage zones (similar to 'Oppel zones') which are optimal in the sense that they give maximal stratigraphic resolution with a minimum of superpositional contradictions. One example of such a contradiction would be a section containing a species A above a species B, while assemblage 1 (containing species A) is placed below assemblage 2 (containing species B). PAST carries out the following steps: 215 1. Residual maximal horizons The method makes the range-through assumption, meaning that taxa are considered to have been present at all levels between the first and last appearance in any section. Then, any samples with a set of taxa that is contained in another sample are discarded. The remaining samples are called residual maximal horizons. The idea behind this throwing away of data is that the absent taxa in the discarded samples may simply not have been found even though they originally existed. Absences are therefore not as informative as presences. 2. Superposition and co-occurrence of taxa Next, all pairs (A,B) of taxa are inspected for their superpositional relationships: A below B, B below A, A together with B, or unknown. If A occurs below B in one locality and B below A in another, they are considered to be co-occurring although they have never actually been found together. The superpositions and co-occurrences of taxa can be viewed in the biostratigraphic graph. In this graph, taxa are coded as numbers. Co-occurrences between pairs of taxa are shown as solid blue lines. Superpositions are shown as dashed red lines, with long dashes from the above-occurring taxon and short dashes from the below-occurring taxon. Some taxa may occur in so-called forbidden sub-graphs, which indicate inconsistencies in their superpositional relationships. Two of the several types of such sub-graphs can be plotted in PAST: Cn cycles, which are superpositional cycles (A->B->C->A), and S3 circuits, which are inconsistencies of the type 'A co-occurring with B, C above A, and C below B'. Interpretations of such forbidden sub-graphs are suggested by Guex (1991). 216 3. Maximal cliques Maximal cliques are groups of co-occurring taxa not contained in any larger group of co-occurring taxa. The maximal cliques are candidates for the status of unitary associations, but will be further processed below. In PAST, maximal cliques receive a number and are also named after a maximal horizon in the original data set which is identical to, or contained in (marked with asterisk), the maximal clique. 4. Superposition of maximal cliques The superpositional relationships between maximal cliques are decided by inspecting the superpositional relationships between their constituent taxa, as computed in step 2. Contradictions (some taxa in clique A occur below some taxa in clique B, and vice versa) are resolved by a 'majority vote'. The contradictions between cliques can be viewed in PAST. The superpositions and co-occurrences of cliques can be viewed in the maximal clique graph. In this graph, cliques are coded as numbers. Co-occurrences between pairs of cliques are shown as solid blue lines. Superpositions are shown as dashed red lines, with long dashes from the above-occurring clique and short dashes from the below-occurring clique. Also, cycles between maximal cliques (see below) can be viewed as green lines. 5. Resolving cycles It will sometimes be the case that maximal cliques are now ordered in cycles: A is below B, which is below C, which is below A again. This is clearly contradictory. The 'weakest link' (superpositional relationship supported by fewest taxa) in such cycles is destroyed. 6. Reduction to unique path At this stage, we should ideally have a single path (chain) of superpositional relationships between maximal cliques, from bottom to top. This is however often not the case, for example if A and B are below C, which is below D, or if we have isolated paths without any relationships (A below B and C below D). To produce a single path, it is necessary to merge cliques according to special rules. 7. Post-processing of maximal cliques Finally, a number of minor manipulations are carried out to 'polish' the result: Generation of the 'consecutive ones' property, reinsertion of residual virtual co-occurrences and superpositions, and compaction to remove any generated non-maximal cliques. For details on these procedures, see Guex (1991). At last, we now have the Unitary Associations, which can be viewed in PAST. 217 The unitary associations have associated with them an index of similarity from one UA to the next, called D: Di=|UAi-UAi-1| / |UAi| +|UAi-1-UAi| / |UAi-1| 8. Correlation using the Unitary Associations The original samples are now correlated using the unitary associations. A sample may contain taxa which uniquely places it in a unitary association, or it may lack key taxa which could differentiate between two or more unitary associations, in which case only a range can be given. These correlations can be viewed in PAST. 9. Reproducibility matrix Some unitary associations may be identified in only one or a few sections, in which case one may consider to merge unitary associations to improve the geographical reproducibility (see below). The reproducibility matrix should be inspected to identify such unitary associations. A UA which is uniquely identified in a section is shown as a black square, while ranges of UAs (as given in the correlation list) are shown in gray. 10. Reproducibility graph and suggested UA merges (biozonation) The reproducibility graph (Gk' in Guex 1991) shows the superpositions of unitary associations that are actually observed in the sections. PAST will internally reduce this graph to a unique maximal path (Guex 1991, section 5.6.3), and in the process of doing so it may merge some UAs. These mergers are 218 shown as red lines in the reproducibility graph. The sequence of single and merged UAs can be viewed as a suggested biozonation. Special functionality The implementation of the Unitary Associations method in PAST includes a number of options and functions which have not yet been described in the literature. For questions about these, please contact us. References Angiolini, L. & H. Bucher. 1999. Taxonomy and quantitative biochronology of Guadalupian brachiopods from the Khuff Formation, Southeastern Oman. Geobios 32:665-699. Guex, J. 1991. Biochronological Correlations. Springer Verlag. Savary, J. & J. Guex. 1999. Discrete Biochronological Scales and Unitary Associations: Description of the BioGraph Computer Program. Meomoires de Geologie (Lausanne) 34. 219 Ranking-Scaling Ranking-Scaling (Agterberg & Gradstein 1999) is a method for quantitative biostratigraphy based on events in a number of wells or sections. The data input consists of wells in rows with one well per row, and events (e.g. FADs and/or LADs) in columns. The values in the matrix are depths of each event in each well, increasing upwards (you may want to use negative values to achieve this). Absences are coded as zero. If only the order of events is known, this can be coded as increasing whole numbers (ranks, with possible ties for co-occurring events) within each well. The implementation of ranking-scaling in PAST is not comprehensive, and advanced users are referred to the RASC and CASC programs of Agterberg and Gradstein. Overview of the method The method of Ranking-Scaling proceeds in two steps: 1. Ranking The first step of Ranking-Scaling is to produce a single, comprehensive stratigraphic ordering of events, even if the data contains contradictions (event A over B in one well, but B over A in another), or longer cycles (A over B over C over A). This is done by 'majority vote', counting the number of times each event occurs above, below or together with all others. Technically, this is achieved by Presorting followed by the Modified Hay Method (Agterberg & Gradstein 1999). 2. Scaling The biostratigraphic analysis may end with ranking, but additional insight may be gained by estimating stratigraphic distances between the consecutive events. This is done by counting the number of observed superpositional relationships (A above or below B) between each pair (A,B) of consecutive events. A low number of contradictions implies long distance. Some computed distances may turn out to be negative, indicating that the ordering given by the ranking step was not optimal. If this happens, the events are re-ordered and the distances recomputed in order to ensure only positive inter-event distances. RASC in PAST Parameters  Well threshold: The minimum number of wells in which an event must occur in order to be included in the analysis 220  Pair threshold: The minimum number of times a relationship between events A and B must be observed in order for the pair (A,B) to be included in the ranking step  Scaling threshold: Pair threshold for the scaling step  Tolerance: Used in the ranking step (see Agterberg & Gradstein) Ranking The ordering of events after the ranking step is given, with the first event at the bottom of the list. The "Range" column indicates uncertainty in the position. Scaling The ordering of the events after the scaling step is given, with the first event at the bottom of the list. For an explanation of all the columns, see Agterberg & Gradstein (1999). Event distribution A plot showing the number of events in each well, with the wells ordered according to number of events. Scattergrams For each well, the depth of each event in the well is plotted against the optimum sequence (after scaling). Ideally, the events should plot in an ascending sequence. Dendrogram Plot of the distances between events in the scaled sequence, including a dendrogram which may aid in zonation. Reference Agterberg, F.P. & F.M. Gradstein. 1999. The RASC method for Ranking and Scaling of Biostratigraphic Events. In: Proceedings Conference 75th Birthday C.W. Drooger, Utrecht, November 1997. Earth Science Review 46(1-4):1-25. 221 Constrained optimization (CONOP) Table of depths/levels, with wells/sections in rows and event pairs in columns: FADs in odd columns and LADs in even columns. Missing events coded with zeros. PAST includes a simple version of Constrained Optimization (Kemple et al. 1989). Both FAD and LAD of each taxon must be specified in alternate columns. Using so-called Simulated Annealing, the program searches for a global (composite) sequence of events that implies a minimal total amount of range extension (penalty) in the individual wells/sections. The parameters for the optimization procedure include an initial annealing temperature, the number of cooling steps, the cooling ratio (percentage lower than 100), and the number of trials per step. For explanation and recommendations, see Kemple et al. (1989). Output windows include the optimization history with the temperature and penalty as function of cooling step, the global composite solution and the implied ranges in each individual section. The implementation of CONOP in PAST is based on a FORTRAN optimization core provided by Sadler and Kemple. Reference Kemple, W.G., P.M. Sadler & D.J. Strauss. 1989. A prototype constrained optimization solution to the time correlation problem. In Agterberg, F.P. & G.F. Bonham-Carter (eds), Statistical Applications in the Earth Sciences. Geological Survey of Canada Paper 89-9:417-425. 222 Appearance Event Ordination NOT YET IN PAST 3 Diversity curve NOT YET IN PAST 3 Range confidence intervals NOT YET IN PAST 3 Distribution-free range confidence intervals NOT YET IN PAST 3 Spindle diagram NOT YET IN PAST 3 Filter events Cladistics NOT YET IN PAST 3 Scripting NOT YET IN PAST 3