require(devtools) install_github("MorphoFun/psa", build_vignettes = TRUE) library(psa)
To see all of the vignettes that are available in psa
, type browseVignettes("psa")
.
data(BumpusMales)
str(BumpusMales)
?BumpusMales
w <- BumpusMales[,1]
z <- BumpusMales[,3:11]
z$Weight_g <- z$Weight_g^(1/3)
The benefit of using psa::glam or psa::gradients over stats::lm or stats::glm is that the psa functions have built-in options to make the regression coefficients more directly comparable to selection gradients. There are four main benefits to using the psa functions:
1. glam and gradients will check whether the data are centered to a mean of zero and standardized to unit variance prior, and will proceed with these standardizations as a default if raw data are entered into the function.
1. Both linear and nonlinear selection are evaluated, and are output in separate lists for glam and a single data.frame for gradients. "GL" = gradients for linear selection, and "GNL" = gradients for nonlinear selection. The linear terms in GNL should be ignored because they represent biased estimates of linear selection and should be ignored. Linear selection gradients should only be evaluated in the GL, which only conducts a first order model.
1. Quadratic regression coefficients and their associated standard errors from standard statistical programs need to be doubled in order to generate appropriate quadratic selection gradients and standard errors. Such an oversight has resulted in substantial underestimation of the strength of quadratic selection (Stincombe et al. 2008). However, glam automatically accounts for this discrepancy, so that the outputted coefficients represent the selection gradients, and do not require further modification.
To demonstrate the functionality of glam, we will be using the male sparrow data from the classic Bumpus (1968) study to reproduce the results from the Lande-Arnold (1983) multiple linear regression method.
# Estimating the gradients using the matrix algebra and OLS methods bm_all_g <- t(gradients(w, z, "all")) # Estimate the standard errors for the matrix algebra approach, using bootstrapping set.seed(123) bm_booted_gMatrix <- gradients_bootstats(w,z,method = 1) # Now, estimating phenotypic selection gradients with the OLS method, using the glam function (to get the se's) bm_gReg <- glam(w, z, "gaussian") # standard errors from the OLS method bm_gReg_se <- c(summary(bm_gReg$GL)$coefficients [-1,2], summary( bm_gReg$GNL)$coefficients [-c(1:(length(z)+1)),2]) # Combining both sets of standard errors for the selection gradients bm_all_g_se <- data.frame(cbind(gMatrix = bm_booted_gMatrix$se, gReg = bm_gReg_se), stringsAsFactors = FALSE) # melt the results based on the method used to compute the estimates library(reshape2) melt_bm_g <- melt(t(bm_all_g), id = names(bm_all_g)) melt_bm_gse <- melt(t(bm_all_g_se), id = names(bm_all_g_se)) limits <- aes(ymax = melt_bm_g$value + melt_bm_gse$value, ymin = melt_bm_g$value - melt_bm_gse$value)
# plot the data # see docs.ggplot2.org/current/geom_point.html ggplot(melt_bm_g, aes(Var1, value, fill = Var1)) + geom_point(aes(colour = factor(Var1))) + geom_errorbar(limits, width = 0.2) + facet_wrap(~Var2) + geom_abline(intercept = 0, slope = 0, colour = "gray") + labs(title = "Selection gradients for male Bumpus data") + theme_bw() + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + theme(legend.position = "none", panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black") )
References
Brodie ED III, Moore AJ, Janzen FJ. 1995. Visualizing and quantifying natural selection. TREE 10(8): 313-318.
Lande SJ, Arnold R. 1983. The measurement o fselection on correlated characters. Evolution 37(6): 1210-1226.
Phillips PC, Arnold SJ. 1989. Visualizing multivariate selection. Evolution 43(6): 1209-1222.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.