Package {lande}

lande: Tools for Evolutionary Selection Analysis

Description

Estimates selection differentials and linear, quadratic and correlational selection gradients on relative fitness following Lande and Arnold (1983) doi:10.1111/j.1558-5646.1983.tb00236.x, with p-values from the logistic model for survival and from a Poisson or negative binomial model for counts. Fits cubic spline fitness functions as in Schluter (1988) doi:10.1111/j.1558-5646.1988.tb02507.x and two-trait fitness surfaces drawn only where there are data, and computes adaptive landscapes of mean fitness against the population mean for one or two traits, by period and for several groups on one surface. Includes bootstrap intervals, a combined results table and a Shiny app.

Data

Each dataset has its own help page, with its source: bumpus, crescent_pond_pupfish and little_lake_pupfish, finch_yearly and finch_community. The pupfish and finch data are also in inst/extdata as CSV files.

Author(s)

Maintainer: James T. Stroud stroud@gatech.edu (affiliation: School of Biological Sciences, Georgia Institute of Technology, Atlanta, GA, United States)

Authors:

• James T. Stroud stroud@gatech.edu (affiliation: School of Biological Sciences, Georgia Institute of Technology, Atlanta, GA, United States)

• Vanessa Zhou (affiliation: College of Computing, Georgia Institute of Technology, Atlanta, GA, United States)

• Ali Sean Khan (affiliation: College of Computing, Georgia Institute of Technology, Atlanta, GA, United States)

• Calvin Yu (affiliation: College of Computing, Georgia Institute of Technology, Atlanta, GA, United States)

• Taylor Papke (affiliation: College of Computing, Georgia Institute of Technology, Atlanta, GA, United States)

• Dima Al Nabhani (affiliation: College of Computing, Georgia Institute of Technology, Atlanta, GA, United States)

• Guillermo Garcia-Costoya (affiliation: School of Biological Sciences, Georgia Institute of Technology, Atlanta, GA, United States)

See Also

Useful links:

• https://github.com/Human-Augment-Analytics/lande

• https://human-augment-analytics.github.io/lande/

• Report bugs at https://github.com/Human-Augment-Analytics/lande/issues

Calculate Adaptive Landscape

Description

Computes the adaptive landscape (mean fitness as a function of population mean phenotype) using a fitted individual-level fitness model.

Usage

adaptive_landscape(
  data,
  fitness_model,
  trait_cols,
  group_col = NULL,
  population_variance = NULL,
  simulation_n = 1000,
  grid_n = 50,
  custom_range = NULL,
  clamp = TRUE,
  support_warn = 0.25
)

Arguments

Details

For a single trait the grid also carries .ind_fit, the individual fitness function evaluated at each population mean, for drawing the two curves together (see plot_adaptive_landscape()).

The simulated populations spread beyond the data, especially towards the edge of the grid, and the fitness of those individuals is extrapolated. .outside in the grid is the share of each simulated population falling outside the convex hull of the observed trait pairs, or outside the observed range for one trait, and the result's support gives that share over the whole grid and at the optimum.

Value

An object of class "adaptive_landscape". Its optimum_edge is TRUE when the highest mean fitness lies on the edge of the grid, where the landscape may keep rising beyond it.

Examples

prep <- prepare_selection_data(bumpus, "survival", c("total_length", "weight"))
surf <- correlated_fitness_surface(prep, "survival", c("total_length", "weight"), grid_n = 30)
land <- adaptive_landscape(prep, surf$model, c("total_length", "weight"),
                           simulation_n = 100, grid_n = 15)
land$optimum

# one trait, from the spline fitness function
uni <- univariate_spline(prep, "survival", "total_length")
land1 <- adaptive_landscape(prep, uni$model, "total_length", simulation_n = 100, grid_n = 30)
plot_adaptive_landscape(land1, "total_length")

Disruptive or stabilising selection on one trait

Description

Estimates the linear (beta) and quadratic (gamma) selection gradients of a single trait. Gamma is positive under disruptive selection and negative under stabilising selection.

Usage

analyze_disruptive_selection(
  data,
  fitness_col,
  trait_col,
  fitness_type = c("auto", "binary", "count", "continuous"),
  standardize = TRUE,
  group = NULL,
  return_grouped = FALSE
)

Arguments

Value

A data frame with one row per gradient (Term, Type, Beta_Coefficient, Standard_Error, P_Value, Variance).

Examples

analyze_disruptive_selection(bumpus, "survival", "total_length")
analyze_disruptive_selection(bumpus, "survival", "total_length",
                             group = "sex", return_grouped = TRUE)

Analyze linear selection gradients (beta)

Description

Estimates linear selection gradients by OLS on relative fitness. Binary fitness gets its p-values from a logistic GLM and count fitness from a Poisson GLM (negative binomial if overdispersed), since the OLS tests are not valid for either.

Usage

analyze_linear_selection(
  data,
  fitness_col,
  trait_cols,
  fitness_type,
  binary_response_col = NULL,
  group = NULL,
  se_type = c("ols", "hc3")
)

Arguments

Value

A list containing the fitted models, summaries, ANOVA tables, and VIFs.

Examples

prep <- prepare_selection_data(bumpus, "survival", c("total_length", "weight"))
fit <- analyze_linear_selection(prep, "survival", c("total_length", "weight"), "binary")
extract_linear_coefficients(c("total_length", "weight"), fit)

Analyze nonlinear selection gradients (gamma)

Description

Estimates quadratic and correlational selection gradients by OLS on relative fitness.

Usage

analyze_nonlinear_selection(
  data,
  fitness_col,
  trait_cols,
  fitness_type,
  binary_response_col = NULL,
  group = NULL,
  se_type = c("ols", "hc3")
)

Arguments

Value

A list containing the fitted nonlinear models, summaries, ANOVA tables, and VIFs.

Examples

prep <- prepare_selection_data(bumpus, "survival", c("total_length", "weight"))
fit <- analyze_nonlinear_selection(prep, "survival", c("total_length", "weight"), "binary")
extract_quadratic_coefficients(c("total_length", "weight"), fit)
extract_interaction_coefficients(c("total_length", "weight"), fit)

Bootstrap selection gradients

Description

Resamples individuals with replacement and re-estimates the selection gradients on each resample to obtain bootstrap standard errors and percentile confidence intervals for every coefficient (beta, gamma, gamma_ij).

Usage

bootstrap_selection(
  data,
  fitness_col,
  trait_cols,
  fitness_type = c("auto", "binary", "count", "continuous"),
  standardize = TRUE,
  group = NULL,
  use_relative_for_fit = TRUE,
  n_boot = 1000,
  conf = 0.95
)

Arguments

Details

Each resample goes through the whole procedure again: with standardize = TRUE the traits are restandardised, and fitness is made relative to the resample's own mean, before both models are refitted, so the intervals carry the uncertainty in the standardisation as well as in the fit. Resamples whose fit fails are dropped, and N_Boot counts the ones that were used.

When group is given, individuals are resampled within each group so every resample keeps the original group sizes. A resample in which a trait has no variance within a group (so it could only be centred) is discarded rather than fitted on a degenerate value. The point estimate is the ordinary selection_coefficients() fit on the full data and its warnings are reported as usual.

With continuous fitness the intervals for beta can be too narrow when the traits have heavy tails (skewed or not) and the linear model leaves out curvature, though less so than the parametric ones (see check_selection_assumptions()). For survival and counts in the package's simulation they covered a little less than the parametric ones, 92 to 93% for beta. Refitting on a transformed trait changes the scale of the gradients and is only a check.

Value

A data frame with one row per coefficient and columns Term, Type, Estimate (point estimate on the full data), Boot_SE, CI_lower, CI_upper, P_Value, and N_Boot (usable resamples for that coefficient). P_Value is the parametric p-value from the point-estimate fit (OLS t-test, or the logistic-GLM Wald test for binary fitness), not a bootstrap p-value.

Examples

set.seed(1)
bootstrap_selection(bumpus, "survival", c("total_length", "weight"), n_boot = 50)

Bumpus house sparrow survival data

Description

Morphological measurements and overwinter survival of house sparrows (Passer domesticus) collected by Hermon Bumpus in 1898 after a severe winter storm in Providence, Rhode Island.

Usage

bumpus

Format

A data frame with 136 rows and 11 variables:

sex

Sex of the bird (male or female).

survival

Survival indicator: 1 if the bird survived, 0 otherwise.

total_length

Total length (mm).

wingspread

Alar extent / wingspread (mm).

weight

Body weight (g).

head_length

Length of head and beak (mm).

humerus

Length of humerus (in).

femur

Length of femur (in).

tibiotarsus

Length of tibiotarsus (in).

skull_width

Width of skull (in).

sternum

Length of sternum / keel (in).

Source

Bumpus, H. C. (1899). The elimination of the unfit as illustrated by the introduced sparrow, Passer domesticus. Biological Lectures, Marine Biological Laboratory, Woods Hole, 209-226.

Canonical analysis of the gamma matrix

Description

Rotates the matrix of quadratic and correlational selection gradients to its canonical axes (Phillips and Arnold 1989; Blows and Brooks 2003). Each axis is a combination of the traits with a single curvature, its eigenvalue: negative for stabilising and positive for disruptive selection along that axis. Correlational selection that is spread over many \gamma_{ij} shows up as curvature on a few axes.

Usage

canonical_analysis(
  data,
  fitness_col,
  trait_cols,
  fitness_type = c("auto", "binary", "count", "continuous"),
  standardize = TRUE,
  group = NULL,
  bootstrap = FALSE,
  n_boot = 200,
  conf = 0.95,
  test = c("permutation", "double_regression"),
  n_perm = 999
)

Arguments

Details

The gamma matrix \gamma comes from selection_coefficients(), with the quadratic gradients doubled and the correlational gradients as they stand. Its eigenvectors are the columns of M and its eigenvalues the curvatures \lambda; directional selection along the axes is \theta = M^T \beta. Standard errors come from the double regression of Bisgaard and Ankenman (1996): fitness is refitted on the canonical scores and their squares, without cross-products, and twice the coefficient of a squared score is that axis's \lambda.

The double-regression tests treat the axes as known when they were estimated from the same data. With 136 individuals and fitness unrelated to five traits, the largest eigenvalue came out significant in about three samples in four, so by default the p-values come from the permutation test of Reynolds et al. (2010): fitness is shuffled among individuals, within groups when group is given, the canonical analysis and the double regression are repeated, and each axis's p-value is the share of shuffles, counting the data as one, whose F for the eigenvalue of the same rank is at least the observed one. It tests whether an axis curves more than it would with fitness unrelated to the traits; once one axis is strongly curved, the tests of the others are only rough. Every fitness type uses least squares on relative fitness, and the result depends on the random seed; call set.seed() first. With test = "double_regression" the p-values come from the double regression itself, for survival and counts from the logistic or count model on the same terms. The standard errors treat the axes as known either way.

Sampling error in \gamma also spreads its eigenvalues, so the largest curvatures are overestimated, especially with many traits and few individuals (Reynolds et al. 2010), and the estimated axes lean towards directions of phenotype with little variance (Morrissey 2014). Read the eigenvalues together with the fitness surface along the same axes, which plot_canonical_axes() draws, and with the bootstrap intervals. With the bootstrap, each resample's axes are matched to the original ones by their largest absolute cosine, since order and sign change between resamples.

Value

An object of class "canonical_analysis": gamma, the gamma matrix; beta; M, the loadings with one column per axis; axes, a data frame with lambda, its standard error and p-value, theta, and with the bootstrap CI_lower, CI_upper and N_Boot; scores, the prepared data with the canonical scores m1, m2, ... added; the fitness type, n, test and n_perm.

References

Bisgaard, S. and Ankenman, B. (1996) Standard errors for the eigenvalues in second-order response surface models. Technometrics 38, 238-246. Blows, M. W. and Brooks, R. (2003) Measuring nonlinear selection. The American Naturalist 162, 815-820. Morrissey, M. B. (2014) In search of the best methods for multivariate selection analysis. Methods in Ecology and Evolution 5, 1095-1109. Phillips, P. C. and Arnold, S. J. (1989) Visualizing multivariate selection. Evolution 43, 1209-1222. Reynolds, R. J., Childers, D. K. and Pajewski, N. M. (2010) The distribution and hypothesis testing of eigenvalues from the canonical analysis of the gamma matrix of quadratic and correlational selection gradients. Evolution 64, 1076-1085.

Examples

ca <- canonical_analysis(bumpus, "survival", c("total_length", "weight", "humerus"))
ca

Check the assumptions behind a selection analysis

Description

Runs the checks a Lande and Arnold analysis rests on and returns them in one table: multivariate normality of the traits, normality of each trait, collinearity, rows per model term, and for the gradient models a residual normality and heteroscedasticity test (continuous fitness), a separation check (binary fitness) or the dispersion ratio (count fitness). Nothing is enforced; the table is there to report alongside the gradients, as the protocol of Palacio et al. (2019) asks.

Usage

check_selection_assumptions(
  data,
  fitness_col,
  trait_cols,
  fitness_type = c("auto", "binary", "count", "continuous"),
  standardize = TRUE,
  group = NULL
)

Arguments

Details

The regression gradients equal the covariance-adjusted differentials, P^{-1} S, by least squares whatever the trait distribution. Normal traits let them be read as the average slope and curvature of the fitness surface (Lande and Arnold 1983; Morrissey and Sakrejda 2013) and tie gamma to the change in the phenotypic covariance. Fitness doesn't need to be normal, and normality doesn't change how the gradients are computed. Mardia's (1970) skewness and kurtosis test multivariate normality, and each trait is also tested on its own with Shapiro and Wilk's test when there are 5000 rows or fewer.

In the package's simulation (validation/simulation.R), on a curved fitness surface, nominal 95% intervals for beta covered 88% of the time with normal traits, 67% with log-normal ones and 77% with symmetric heavy-tailed ones (92, 80 and 91% by bootstrap). Most of that loss comes from the curvature, and HC3 standard errors (selection_coefficients(se_type = "hc3")) recover most of it. Estimating the SD of a heavy-tailed trait lowers coverage too, for gamma and, with skewed traits, for beta even on a straight surface (90%), and HC3 doesn't help there. Skew also biases beta a little. Refitting with a transformed trait changes the scale selection is measured on, so report it with the untransformed result. Collinearity is the largest variance inflation factor from the linear gradient model. The rows-per-term check counts the individuals, or for binary fitness the rarer outcome, per term of the quadratic model, with ten as the working minimum. With more than 2000 individuals Mardia's statistics use a random subsample of 2000, drawn the same way every run. Like the gradient models, the logistic and count models fit an intercept per group, with unlabelled rows as one more group. A group in which every individual survived, or none did, is fitted at 0 or 1 by its intercept and is left out of the separation check. For counts the note names the model behind each set of p-values; the linear and quadratic models are checked for overdispersion separately. If the performance package is installed its heteroscedasticity and overdispersion tests are added next to the package's own, along with an R squared for the gradient model.

Value

A data frame of class "selection_assumptions" with one row per check: check, statistic, p_value and note.

References

Lande, R. and Arnold, S. J. (1983) The measurement of selection on correlated characters. Evolution 37, 1210-1226. Mardia, K. V. (1970) Measures of multivariate skewness and kurtosis with applications. Biometrika 57, 519-530. Morrissey, M. B. and Sakrejda, K. (2013) Unification of regression-based methods for the analysis of natural selection. Evolution 67, 2094-2100. Palacio, F. X., Ordano, M. and Benitez-Vieyra, S. (2019) Measuring natural selection on multivariate phenotypic traits: a protocol for verifiable and reproducible analyses of natural selection. Israel Journal of Ecology and Evolution 65, 130-136.

Examples

check_selection_assumptions(bumpus, "survival", c("total_length", "weight", "humerus"))

Compare Correlated Fitness Surface vs Adaptive Landscape

Description

Compares individual-level fitness (correlated surface) to population-level mean fitness (adaptive landscape).

Usage

compare_fitness_surfaces_data(
  correlated_surface,
  adaptive_landscape,
  trait_cols,
  calculate_correlation = TRUE
)

Arguments

Value

A list containing combined grid data, summary stats, optima points, and correlation metrics.

Examples

prep <- prepare_selection_data(bumpus, "survival", c("total_length", "weight"))
surf <- correlated_fitness_surface(prep, "survival", c("total_length", "weight"), grid_n = 30)
land <- adaptive_landscape(prep, surf$model, c("total_length", "weight"),
                           simulation_n = 100, grid_n = 15)
comp <- compare_fitness_surfaces_data(surf, land, c("total_length", "weight"))
plot_fitness_surfaces_comparison(comp)

Calculate the Correlated Fitness Surface

Description

Fits a model for individual fitness based on multiple individual phenotypes (w ~ z1 + z2 + interactions). Traits MUST be standardized BEFORE calling this function.

Usage

correlated_fitness_surface(
  data,
  fitness_col,
  trait_cols,
  grid_n = 60,
  method = "auto",
  scale_traits = FALSE,
  group = NULL,
  group_effect = TRUE,
  k = NULL,
  mask = TRUE,
  too_far = NULL,
  bs = c("tp", "cr", "ps"),
  smoothing = c("REML", "GCV.Cp", "ML"),
  clamp = TRUE,
  level = 0.95,
  by_group = FALSE,
  count_family = c("poisson", "quasipoisson", "nb")
)

Arguments

Details

The family follows the fitness column: 0/1 fitness gets a binomial family, non-negative whole numbers with more than two values (lifespan in years, offspring) a Poisson family with a log link or the one set by count_family, and anything else a Gaussian family. With method = "auto" binary and count fitness use the GAM and continuous fitness the thin-plate spline. By default the GAM uses a thin-plate smooth with the smoothing parameter chosen by REML; bs and smoothing are there to match another study's smoother and do not apply to method = "tps". With "cr" or "ps", which are one-dimensional bases, the two traits enter as a tensor product smooth. Predicting a fitted surface over the full rectangle of the grid extrapolates into corners that no individual occupies; mask = TRUE leaves those blank rather than showing a fitted value there. The hull still fills gaps between separate clusters of individuals, such as several species on one surface; too_far blanks those too.

Value

A list containing the fitted model, grid predictions, and metadata. With method = "gam" the grid also carries .se, the standard error of the fitted fitness, and .fit_lo and .fit_hi, the band at level, worked out on the scale of the link and blanked where the fit is; the thin-plate spline gives none. peaks lists the local maxima of the kept surface, highest first, with interior TRUE when every neighbouring cell is kept and lower, and FALSE when the maximum sits against the edge of the data, which may only be where the data end. original_data holds the rows that were analysed. With a group the groups data frame has one row per group with its size, mean traits and the highest point of the surface within its own range, flagged peak_interior when that point is an interior maximum of the surface and peak_edge when it lies at the edge of the data; with both FALSE the surface keeps rising past the group's range. For count fitness count_family and dispersion record the family used and the Pearson dispersion of the fit.

Examples

prep <- prepare_selection_data(bumpus, "survival", c("total_length", "weight"))
surf <- correlated_fitness_surface(prep, "survival", c("total_length", "weight"), grid_n = 30)
surf$grid[which.max(surf$grid$.fit), ]

# two lakes of pupfish on one surface: cells far from any fish blank, each
# lake's mean and local peak marked, the lake kept out of the model
pup <- rbind(crescent_pond_pupfish, little_lake_pupfish)
pup <- pup[pup$density == "H", ]
prep2 <- prepare_selection_data(pup, "survival", c("nose", "noseangle"))
surf2 <- correlated_fitness_surface(prep2, "survival", c("nose", "noseangle"), grid_n = 30,
                                    too_far = 0.15, group = "lake", group_effect = FALSE)
surf2$groups

Pupfish hybrid survival in two San Salvador lakes

Description

Laboratory-reared F2 hybrids of the three pupfish species of San Salvador Island, Bahamas, photographed, tagged and released into field enclosures in Crescent Pond and Little Lake from March to June 2011, with laboratory-reared fish of the three species measured the same way. Martin (2016) analysed the high-density enclosures, density == "H", as do the examples and the app.

Usage

crescent_pond_pupfish

little_lake_pupfish

Format

Data frames with 23 variables, 993 rows for Crescent Pond and 1062 for Little Lake:

lake

"CP" (Crescent Pond) or "LL" (Little Lake).

density

The enclosure of an F2 hybrid, "H" (high density) or "L" (low density), or the species of a parental fish: "norm", the generalist Cyprinodon variegatus; "bozo", the molluscivore C. brontotheroides; "bull", the scale-eater C. desquamator.

survival

1 if the hybrid survived the three months in the enclosure, 0 otherwise; NA for parental fish.

ln.growth

Log growth rate of survivors over the three months. Hybrids that died are coded 0, so drop them before analysing growth; NA for parental fish.

color

Martin's plotting colour: orange for hybrids; blue, green and red for the generalist, molluscivore and scale-eater.

d13C, d15N

Carbon and nitrogen stable isotope ratios (per mil) of muscle from surviving hybrids; NA for the rest.

jaw

Lower jaw length, from the jaw joint to the tip of the dentary.

eye

Eye diameter, the mean of the major and minor axes of the iris (orbit diameter).

eye2

Eye roundness, the ratio of the major to the minor axis of the iris.

pmx

Craniofacial height, from the jaw joint to the tip of the premaxilla (upper jaw length).

snout

Lateral snout length, from the tip of the premaxilla to the front of the iris.

body

Dorsal to anal distance, between the first rays of the dorsal and anal fins (body depth).

caudal

Caudal peduncle height.

nasal

Dorsal snout length, from the front of the orbit to the tip of the maxilla in dorsal view.

mouth

Buccal width in dorsal view.

width

Head width across the opercula in dorsal view.

boteyeangle

Ventral orbit angle, at the tip of the premaxilla between the lower edge of the iris and the quadrate.

topeyeangle

Orbit angle, at the tip of the premaxilla between the upper and lower edges of the iris.

nose

Maxillary head protrusion (nasal protrusion).

noseangle

Maxillary head protrusion angle (nasal angle).

adduct

Preopercular height.

SL

Standard length.

An object of class data.frame with 1062 rows and 23 columns.

Details

The trait columns are Martin's scores for the 16 measurements of Martin and Wainwright (2013), taken from the photographs. He log-transformed the distances, regressed each trait on a size index (the first principal component of seven size-related distances) and kept the residuals, fitting the regression to the parental fish alone for jaw, mouth, nasal and width; eye2, nose and noseangle were not size-corrected. Every trait is standardised to mean 0 and standard deviation 1 within each lake, hybrids and parental fish together.

The column names are Martin's. The trait names below are the ones defined in the supplement to Martin and Wainwright (2013), matched to the columns through his size-correction script (in the Dryad archive) and their discriminant loadings, with the name Martin (2016) uses in brackets where it differs. The six functional traits of his Table 3 are jaw, pmx, nose, noseangle, body and eye. Their loadings table swaps the labels of snout and nasal; the names here follow the size-correction script.

Source

Martin, C. H. (2016) Context dependence in complex adaptive landscapes: frequency and trait-dependent selection surfaces within an adaptive radiation of Caribbean pupfishes. Evolution 70, 1265-1282. doi:10.1111/evo.12932. Data: doi:10.5061/dryad.n3mj3.

References

Martin, C. H. and Wainwright, P. C. (2013) Multiple fitness peaks on the adaptive landscape drive adaptive radiation in the wild. Science 339, 208-211. doi:10.1126/science.1227710.

Automatically detect fitness type and GLM family

Description

Heuristically determines whether a fitness vector represents binary, count, proportion, or continuous data.

Usage

detect_family(y)

Arguments

Value

A list containing type (string), family (GLM family object), and note (string).

Examples

detect_family(c(0, 1, 0, 0, 1, 1))
detect_family(bumpus$survival)$type
detect_family(rpois(50, 2))$type

Extract correlational selection coefficients

Description

Extract correlational selection coefficients

Usage

extract_interaction_coefficients(trait_cols, results)

Arguments

Value

A data frame with correlational (interaction) selection coefficients and statistics.

Examples

prep <- prepare_selection_data(bumpus, "survival", c("total_length", "weight"))
fit <- analyze_nonlinear_selection(prep, "survival", c("total_length", "weight"), "binary")
extract_interaction_coefficients(c("total_length", "weight"), fit)

Extract linear selection coefficients

Description

Extract linear selection coefficients

Usage

extract_linear_coefficients(trait_cols, results)

Arguments

Value

A data frame with linear selection coefficients and statistics.

Examples

prep <- prepare_selection_data(bumpus, "survival", c("total_length", "weight"))
fit <- analyze_linear_selection(prep, "survival", c("total_length", "weight"), "binary")
extract_linear_coefficients(c("total_length", "weight"), fit)

Extract quadratic selection coefficients

Description

Extract quadratic selection coefficients

Usage

extract_quadratic_coefficients(trait_cols, results)

Arguments

Value

A data frame with quadratic selection coefficients (doubled estimates) and statistics.

Examples

prep <- prepare_selection_data(bumpus, "survival", c("total_length", "weight"))
fit <- analyze_nonlinear_selection(prep, "survival", c("total_length", "weight"), "binary")
extract_quadratic_coefficients(c("total_length", "weight"), fit)

Darwin's finch community at El Garrapatero

Description

Four ground finch species at El Garrapatero, Santa Cruz Island, Galapagos, with Geospiza fortis split into its small and large beak morphs. One row per bird, with its mean beak measurements and its apparent lifespan, the fitness measure of Beausoleil et al. (2023). Birds first caught late in the study had fewer years in which to be seen again. Built by data-raw/finch_community.R from the file in the authors' code repository (GPL-3).

Usage

finch_community

Format

A data frame with 3428 rows and 7 variables:

band

Band code of the bird.

species

A factor: "fortis small" and "fortis large" (G. fortis), "fuliginosa" (G. fuliginosa), "magnirostris" (G. magnirostris) or "scandens" (G. scandens).

beak_length

Mean beak length (mm).

beak_depth

Mean beak depth (mm).

beak_width

Mean beak width (mm).

lifespan

Apparent lifespan in years: the last year the bird was seen minus the first.

first_year

The year the bird was first caught.

Source

Beausoleil, M.-O. et al. (2023) The fitness landscape of a community of Darwin's finches. Evolution 77, 2533-2546. doi:10.1093/evolut/qpad160. Data: doi:10.5683/SP3/0YIWSE.

Medium ground finch survival from year to year

Description

Medium ground finches (Geospiza fortis) marked and recaptured at El Garrapatero, Santa Cruz Island, Galapagos. Each row is a bird in a year it was seen, from 2004 to 2010, with whether it was seen again in any later year (apparent survival). Built by data-raw/finch_yearly.R.

Usage

finch_yearly

Format

A data frame with 1087 rows and 7 variables:

band

Band code of the bird.

year

Year the bird was seen.

survived

1 if the bird was seen again in any later year up to 2018, the last in the file, 0 otherwise. Death cannot be told from emigration.

beak_pc1

Beak size: the first principal component of the three beak measurements over all the birds, signed so that larger beaks score higher.

beak_length

Median beak length (mm).

beak_width

Median beak width (mm).

beak_depth

Median beak depth (mm).

Source

Beausoleil, M.-O. et al. (2019) Temporally varying disruptive selection in the medium ground finch (Geospiza fortis). Proceedings of the Royal Society B 286, 20192290. doi:10.1098/rspb.2019.2290. Data: doi:10.5061/dryad.zcrjdfn6q.

Compare the fitted fitness at two points of a surface

Description

Takes two points on a GAM surface, given as group names or trait values, and returns the difference in fitted fitness on the scale of the link, with its standard error. Two peaks are separate only if each is higher than the pass between them, the lowest point on the highest route from one to the other, which valley = TRUE compares.

Usage

peak_difference(
  surface,
  from,
  to,
  valley = FALSE,
  route = c("pass", "line"),
  n_path = 50
)

Arguments

Details

Differences are on the scale of the link, log fitness for counts and log odds for survival. The standard error comes from the covariance of the model's coefficients, corrected for smoothing parameter uncertainty when the fit was by REML or ML. The points are taken as fixed, although the peaks and the valley were found on the same fitted surface, so the comparisons describe the fitted surface and are not planned tests. The valley is the lowest point on the route, so its differences are biased upwards and their z values too large, which is why the valley rows have no p-value. The pass is found cell by cell over the grid, joining each kept cell to its eight neighbours, from the two kept cells nearest the points in grid steps; its precision follows grid_n, and a point far from any kept cell gets a message. If the route never drops below the cell a point starts from, or the lower point is no higher than the pass, there is no dip to compare. The straight line ignores the mask. For overdispersed counts fit the surface with count_family = "quasipoisson" first, since Poisson standard errors are then too small.

Value

A data frame with one row per comparison: the fitted fitness at the two points (fit_a, fit_b), their difference on the link scale, its standard error, z and the two-sided p_value (NA for the valley rows). With valley = TRUE the attribute "valley" holds the trait values of the pass or of the lowest point on the line.

Examples

prep <- prepare_selection_data(bumpus, "survival", c("total_length", "weight"))
surf <- correlated_fitness_surface(prep, "survival", c("total_length", "weight"), grid_n = 30)
peak_difference(surf, c(-1, -1), c(1, 1))

Plot Adaptive Landscape

Description

Two traits give a contour map of mean fitness against the two population means. One trait gives a curve of mean fitness against the population mean, with the individual fitness function drawn alongside it for comparison.

Usage

plot_adaptive_landscape(
  landscape,
  trait_cols,
  original_data = NULL,
  group_col = NULL,
  bins = 12,
  show_optimum = TRUE,
  show_actual_means = TRUE,
  show_individual = TRUE,
  point_alpha = 0.8,
  show_support = FALSE,
  support_level = 0.5,
  ...
)

Arguments

Value

A ggplot object representing the adaptive landscape.

Examples

prep <- prepare_selection_data(bumpus, "survival", c("total_length", "weight"))
surf <- correlated_fitness_surface(prep, "survival", c("total_length", "weight"), grid_n = 30)
land <- adaptive_landscape(prep, surf$model, c("total_length", "weight"),
                           simulation_n = 100, grid_n = 15)
plot_adaptive_landscape(land, c("total_length", "weight"))

Plot Adaptive Landscape (3D Perspective)

Description

Plot Adaptive Landscape (3D Perspective)

Usage

plot_adaptive_landscape_3d(
  landscape,
  trait_cols,
  theta = -30,
  phi = 30,
  grid_n = 200,
  color_palette = NULL,
  ...
)

Arguments

Value

A 3D plot produced by fields::drape.plot().

Examples

if (requireNamespace("fields", quietly = TRUE)) {
  prep <- prepare_selection_data(bumpus, "survival", c("total_length", "weight"))
  surf <- correlated_fitness_surface(prep, "survival", c("total_length", "weight"), grid_n = 30)
  land <- adaptive_landscape(prep, surf$model, c("total_length", "weight"),
                             simulation_n = 100, grid_n = 15)
  plot_adaptive_landscape_3d(land, c("total_length", "weight"))
}

Fitness along the canonical axes

Description

Draws the fitness surface in the space of two canonical axes, or the fitness function along one, from the scores of canonical_analysis(). A negative \lambda is stabilising selection along an axis only when fitness along that axis has a peak inside the data, otherwise it is only curvature.

Usage

plot_canonical_axes(ca, which = 1:2, grid_n = 40, ...)

Arguments

Value

A ggplot object.

Examples

ca <- canonical_analysis(bumpus, "survival", c("total_length", "weight", "humerus"))
plot_canonical_axes(ca, which = c(1, 3), grid_n = 25)

Plot Correlated Fitness Surface

Description

Plot Correlated Fitness Surface

Usage

plot_correlated_fitness(
  tps,
  trait_cols,
  bins = 12,
  point_alpha = 0.7,
  show_points = FALSE,
  show_optimum = TRUE,
  show_groups = TRUE,
  group_lines = TRUE,
  uncertainty = c("none", "se", "band"),
  ...
)

Arguments

Value

A ggplot object representing the correlated fitness surface.

Examples

prep <- prepare_selection_data(bumpus, "survival", c("total_length", "weight"))
surf <- correlated_fitness_surface(prep, "survival", c("total_length", "weight"), grid_n = 30)
plot_correlated_fitness(surf, c("total_length", "weight"), show_points = TRUE)

Plot Enhanced Correlated Fitness Surface

Description

Plot Enhanced Correlated Fitness Surface

Usage

plot_correlated_fitness_enhanced(
  tps,
  trait_cols = NULL,
  original_data = NULL,
  fitness_col = NULL,
  bins = 12,
  point_alpha = 0.7,
  show_optimum = TRUE,
  show_groups = TRUE,
  group_lines = TRUE,
  uncertainty = c("none", "se", "band"),
  ...
)

Arguments

Value

A ggplot object with enhanced visualizations.

Examples

prep <- prepare_selection_data(bumpus, "survival", c("total_length", "weight"))
surf <- correlated_fitness_surface(prep, "survival", c("total_length", "weight"), grid_n = 30)
plot_correlated_fitness_enhanced(surf, c("total_length", "weight"),
                                 original_data = prep, fitness_col = "survival")

Plot Fitness Surfaces Comparison

Description

Visualizes the comparison between individual-level correlated fitness and population-level adaptive landscape.

Usage

plot_fitness_surfaces_comparison(
  comparison_data,
  bins = 10,
  title = NULL,
  point_alpha = 0.7,
  linewidth = 0.8,
  show_optima = TRUE
)

Arguments

Value

A list containing side_by_side (a patchwork object) and overlay (a ggplot object).

Examples

prep <- prepare_selection_data(bumpus, "survival", c("total_length", "weight"))
surf <- correlated_fitness_surface(prep, "survival", c("total_length", "weight"), grid_n = 30)
land <- adaptive_landscape(prep, surf$model, c("total_length", "weight"),
                           simulation_n = 100, grid_n = 15)
comp <- compare_fitness_surfaces_data(surf, land, c("total_length", "weight"))
plot_fitness_surfaces_comparison(comp)

Plot fitness functions and landscapes over time

Description

Plot fitness functions and landscapes over time

Usage

plot_temporal_landscape(
  tl,
  type = c("panels", "heatmap"),
  show_points = TRUE,
  show_landscape = TRUE,
  show_optimum = TRUE,
  connect = FALSE,
  bins = 10,
  ncol = NULL,
  ...
)

Arguments

Value

A ggplot object.

Examples

prep <- prepare_selection_data(finch_yearly, "survived", "beak_pc1")
years <- temporal_landscape(prep, "survived", "beak_pc1", "year", landscape = FALSE)
plot_temporal_landscape(years, ncol = 4)
plot_temporal_landscape(years, type = "heatmap")

Plot Univariate Correlated Fitness Function

Description

Visualizes the univariate fitness function estimated via spline.

Usage

plot_univariate_fitness(
  uni,
  trait_col,
  title = NULL,
  show_points = TRUE,
  point_alpha = 0.25,
  ribbon_fill = NA,
  ribbon_linetype = "dashed",
  classic_plot = FALSE,
  ...
)

Arguments

Value

A ggplot object representing the univariate fitness function.

Examples

prep <- prepare_selection_data(bumpus, "survival", "total_length")
uni <- univariate_spline(prep, "survival", "total_length")
plot_univariate_fitness(uni, "total_length")
plot_univariate_fitness(uni, "total_length", classic_plot = TRUE)

Prepare data for selection analysis

Description

Cleans, standardizes, and calculates relative fitness for trait and fitness data. Standardizations and relative fitness calculations can be performed within groups.

Usage

prepare_selection_data(
  data,
  fitness_col,
  trait_cols,
  standardize = TRUE,
  group = NULL,
  add_relative = TRUE,
  na_action = c("warn", "drop", "none"),
  name_relative = "relative_fitness"
)

Arguments

Value

A modified data frame ready for selection analysis.

Examples

prep <- prepare_selection_data(bumpus, "survival", c("total_length", "weight"))
head(prep[, c("survival", "relative_fitness", "total_length", "weight")])

# standardised and relativised within each sex
by_sex <- prepare_selection_data(bumpus, "survival", "total_length", group = "sex")
tapply(by_sex$total_length, by_sex$sex, mean)

Print method for adaptive landscape

Description

Print method for adaptive landscape

Usage

## S3 method for class 'adaptive_landscape'
print(x, ...)

Arguments

Value

The input object x, invisibly.

Print the assumption checks

Description

Print the assumption checks

Usage

## S3 method for class 'selection_assumptions'
print(x, ...)

Arguments

Value

The input object x, invisibly.

Print a selection report

Description

Print a selection report

Usage

## S3 method for class 'selection_report'
print(x, ...)

Arguments

Value

The input object x, invisibly.

Open the Shiny app

Description

Starts the app that ships with the package: pick a bundled dataset or upload a CSV, choose the fitness column, the traits and an optional group, and read the gradients, fitness functions, surfaces and landscapes from the tabs. Needs the shiny package; plotly adds the rotatable 3D landscape and fields the static one.

Usage

run_app(...)

Arguments

Value

Whatever shiny::runApp() returns; called to open the app.

Examples

if (interactive()) run_app()

Calculate selection coefficients

Description

Main wrapper function for selection analysis following Lande & Arnold (1983). Extracts linear (beta), quadratic (gamma), and correlational selection gradients.

Usage

selection_coefficients(
  data,
  fitness_col,
  trait_cols,
  fitness_type = c("auto", "binary", "count", "continuous"),
  standardize = TRUE,
  group = NULL,
  use_relative_for_fit = TRUE,
  return_grouped = FALSE,
  se_type = c("ols", "hc3")
)

Arguments

Value

A data frame containing selection coefficients (Term, Type, Beta_Coefficient, Standard_Error, P_Value, Variance), with the standard errors used in the attribute "se_type".

References

MacKinnon, J. G. and White, H. (1985) Some heteroskedasticity-consistent covariance matrix estimators with improved finite sample properties. Journal of Econometrics 29, 305-325. Mitchell-Olds, T. and Shaw, R. G. (1987) Regression analysis of natural selection: statistical inference and biological interpretation. Evolution 41, 1149-1161.

Examples

selection_coefficients(bumpus, "survival", c("total_length", "weight"), fitness_type = "binary")

# one set of gradients per sex, each sex standardised on its own
selection_coefficients(bumpus, "survival", c("total_length", "weight"),
                       group = "sex", return_grouped = TRUE)

Calculate selection differential (S)

Description

Computes the selection differential S = Cov(z, w) for a single trait.

Usage

selection_differential(
  data,
  fitness_col,
  trait_col,
  standardized = TRUE,
  use_relative = TRUE,
  group = NULL,
  return_grouped = FALSE
)

Arguments

Value

A single numeric value representing overall S, or a data frame with S per group if return_grouped = TRUE.

Examples

prep <- prepare_selection_data(bumpus, "survival", "total_length")
selection_differential(prep, "survival", "total_length")
selection_differential(prep, "survival", "total_length", group = "sex", return_grouped = TRUE)

Standardised selection analysis report

Description

Runs a Lande-Arnold analysis and returns one table of selection differentials (S) and linear (beta), quadratic (gamma) and correlational (gamma_ij) gradients. All estimates use traits standardised to mean 0, SD 1 and relative fitness, so tables from different studies can be compared.

Usage

selection_report(
  data,
  fitness_col,
  trait_cols,
  fitness_type = c("auto", "binary", "count", "continuous"),
  standardize = TRUE,
  group = NULL,
  use_relative_for_fit = TRUE,
  include_differentials = TRUE,
  digits = 4,
  se_type = c("ols", "hc3")
)

Arguments

Details

Differentials and gradients are computed on the same individuals (those with complete fitness and trait values) and on the same trait and fitness scales. S is the population covariance (divides by n) while the OLS gradient on sd-standardised traits corresponds to the sample covariance (n - 1), so for a single trait the Differential and Linear rows differ by the factor (n - 1) / n.

Value

A data frame of class "selection_report" with columns Term, Type, Estimate, Std_Error, and P_Value.

Examples

selection_report(bumpus, "survival", c("total_length", "weight", "humerus"))

Fitness functions and adaptive landscapes over time

Description

Fits the fitness function (one trait) or fitness surface (two traits) separately for each level of a time column, usually year, and optionally the adaptive landscape for each, to show how selection changes over time. Beausoleil et al. (2019) did this for Darwin's finches, one year at a time.

Usage

temporal_landscape(
  data,
  fitness_col,
  trait_cols,
  time_col,
  fitness_type = c("auto", "binary", "count", "continuous"),
  min_n = 20,
  landscape = TRUE,
  k = NULL,
  bs = NULL,
  smoothing = NULL,
  bootstrap = FALSE,
  n_boot = 200,
  grid_n = 60,
  simulation_n = 300,
  mask = TRUE,
  too_far = NULL,
  count_family = c("poisson", "quasipoisson", "nb")
)

Arguments

Value

An object of class "temporal_landscape": fits and landscapes, one per period; summary, one row per period with n, mean fitness, trait means, the smooth's degrees of freedom, the position and height of the highest fitted fitness and whether it sits at the edge of the data (optimum_edge), the number of interior peaks (one trait) and the landscape optimum; grid, the fitted values of every period stacked with a time column; for one trait heat, each period's fitness function on one common trait grid, NA outside that period's data; landscape_grid, the landscapes stacked; data, the rows used; and skipped.

Examples

prep <- prepare_selection_data(finch_yearly, "survived", "beak_pc1")
years <- temporal_landscape(prep, "survived", "beak_pc1", "year", landscape = FALSE)
years$summary
plot_temporal_landscape(years, type = "heatmap")

Estimate univariate correlated fitness function

Description

This function calculates a univariate correlated fitness function using a GAM smooth term. The formula is w ~ f(z). Traits should be standardized BEFORE calling this function.

Usage

univariate_spline(
  data,
  fitness_col,
  trait_col,
  fitness_type = c("auto", "binary", "count", "continuous"),
  group = NULL,
  relative_col = NULL,
  k = 10,
  bs = c("cr", "tp", "ps"),
  smoothing = c("GCV.Cp", "REML", "ML"),
  bootstrap = FALSE,
  n_boot = 1000,
  by_group = FALSE,
  count_family = c("poisson", "quasipoisson", "nb")
)

Arguments

Details

By default the fitness function is a penalised cubic regression spline with the smoothing parameter chosen by generalised cross-validation (UBRE for binary and Poisson fitness), following Schluter (1988). bs and smoothing set the family of curves and the criterion, to match the smoother of another study; the amount of smoothing is still estimated from the data, though the criterion can change it a good deal. If mgcv's check suggests the basis dimension was too small a warning says so; the check works on the residuals, which for survival carry too little to test, so it is skipped there. With bootstrap = TRUE the result depends on the random seed; call set.seed() first for a reproducible ribbon.

Value

A list of class "univariate_fitness" containing the fitted GAM model, a prediction grid, and metadata.

Examples

prep <- prepare_selection_data(bumpus, "survival", "total_length")
uni <- univariate_spline(prep, "survival", "total_length")
plot_univariate_fitness(uni, "total_length")

# a thin-plate basis with REML, to match another study's smoother
univariate_spline(prep, "survival", "total_length", bs = "tp", smoothing = "REML")$spline_type