################################################################################ ### R BASICS WORKSHOP ### ### PRESENTATION 10: WRITING YOUR OWN FUNCTIONS ### ### ### ### Center for Conservation and Sustainable Development ### ### Missouri Botanical Garden ### ### Website: rbasicsworkshop.weebly.com ### ################################################################################ # In addition to using the functions in the packages that are available in R, # you can write your own functions! # These functions have the same properties that other functions in R. Let's # start checking the help for the function *function* # * * Function help ("function") # The basic structure of a function definition is: # function.name <- function(argument1, argument2) # { # code # result # } ## Easy EXAMPLE 1 ## my.fun # There is no function with this name yet my.fun <- function(x) { result <- x + 10 result } # Now use this feature: my.fun(4) my.fun(34) ## Easy EXAMPLE 2 ## # This function calculates the t statistic for comparing the average of two # Vectors: my.twosam # Still there is no function with this name my.twosam <- function(y1, y2) # two arguments without pre-determined values { # Sample Size n1 <- length(y1) n2 <- length(y2) # Averages b1 <- mean(y1) b2 <- mean(y2) # Variances s1 <- var(y1) s2 <- var(y2) # Joint Variance s <- ((n1-1) * s1 + (n2-1) * s2) / (n1 + n2-2) # t statistic tst <- (b1 - b2) / sqrt(s * (1 / n1 + 1 / n2)) # Result exported function tst } # Now use this function to calculate t statistics using two vectors: data(iris) PL.versi <- iris$Petal.Length[which(iris$Species == "versicolor")] PL.seto <- iris$Petal.Length[which(iris $ Species == "setosa")] boxplot(PL.versi, PL.seto) tstat <- my.twosam(PL.versi, PL.seto) tstat # Many functions in R are written in R, and the code behind these them can be # accessed by writing the name of the function in the console lm ## Example 3 ## # The following code creates a function to simulate the dynamics of a # population according to the Ricker model. This model includes the following # parameters used to define the function arguments: # "nzero": the initial population size # "r": the growth rate # k": carrying capacity # "time": is the total number of time units over which the population will be simulated. my.ricker.fun <- function(nzero, r, K, time) { N <- numeric(time + 1) N[1] <- nzero for(i in 1:time) { N[i + 1] <- N[i] * exp(r * (1 - (N[i] / K))) } Time <- 0:time plot(Time, N, type = "l", xlim = c(0, time), cex.axis = 1.5, cex.lab = 1.5, bty = "n", lwd = 2) abline(h = K, lty = 3) N } par(mfrow = c (1,2)) sim.abunds.1 <- my.ricker.fun(nzero = 1, r = 0.1, K = 30, time = 100) sim.abunds.2 <- my.ricker.fun(nzero = 1, r = 0.05, K = 500, time = 100)