stephenrkeller / accel

accel: R Functions to analyze accelerometer data

accel: R Functions to analyze accelerometer data forked from andy gougherty's github repo (agougher/accel)

This package includes functions to read in and analyze accelerometer data generated by Oregon Research Electronic AL100 Acceleration Loggers to estimate tree phenology.

Installation

#Install from github
devtools::install_github("stephenrkeller/accel")

#Load packages
require(accel)
require(forecast)
require(pracma)
require(gtools)

Raw data

To get a sense of what raw accelerometer data looks like, there are several files available here. If you download these files to your working directory, you can analyze them with the code below. The season-length example data below was processed in the same way, but for about 250 files.

#List files in order and calculates dominant period
x <- mixedsort(list.files("./", pattern="CSV", full.names=TRUE))
ex1 <- lapply(x,readAndDom)
ex1 <- data.frame(do.call(rbind, ex1))

#Plotting
plot(1:5, ex1$peakY,pch=20, col="red")
points(1:5, ex1$peakZ, pch=20)

Season-length example

The example data was generated by analyzing 250 acceleration files collected between March 1, 2016 and November 5, 2016, using the code below. The raw data was not included here because of its size, so the code is for illustrative purposes only. readAndDom reads in accelerometer files and calculates the dominant period of the three axes, as well as a weight. Reading in and calculating the dominant periods will take a while, so it is best to run in parallel when possible.

#NOT RUN
#x is a list of ordered csv file locations
al3 <- lapply(x,readAndDom)
al3 <- do.call(rbind, al3)

#Assign day of year
al3["day"] <- 61:310

The example data includes the dominant period values for each of three dimensional axes (x, y, z), and weights reflective of the difference between the dominant and second-most dominant periods for each axis, for a single accelerometer.

#Load example data and make long form dataframe for functions below
data(al3)
head(al3)

##   day    peakX    peakY    peakZ     wpeakX    wpeakY     wpeakZ
## 1  61 55.44444 52.52632 52.52632 0.20141298 0.1585592 0.13467467
## 2  62 52.52632 15.59375 17.50877 0.08484053 0.1573944 0.21712296
## 3  63 52.52632 52.52632 49.90000 0.16278577 0.1131661 0.04290269
## 4  64 19.19231 49.90000 18.48148 0.09507804 0.1368361 0.07151376
## 5  65 19.56863 15.12121 16.63333 0.27714685 0.4667288 0.40435821
## 6  66 19.56863 19.19231 19.96000 0.33783686 0.5037416 0.61549694

#Extract y and z axes, and bind them for the functions below
domY <- al3[,c("day","peakY","wpeakY")]
domZ <- al3[,c("day","peakZ","wpeakZ")]
colnames(domY) <- c("day","value","weights")
colnames(domZ) <- c("day","value","weights")
dat <- rbind(domY,domZ)

#Plot data, point size is proportional to the weight
plot(dat$day, dat$value, cex=dat$weights*2, xlab="Day of year",ylab="Dominant period")

Removing outliers

This function, rmOut, identifies and removes points that have residuals beyond 1.5x the interquartile range of a loess curve. In the plot below, points in purple are kept after removing outliers.

#The 'rmOut' function requires a dataframe with three columns named "day","value", and "weights"
dat2 <- rmOut(dat)
plot(dat$day, dat$value, xlab="Day of year",ylab="Dominant period")
points(dat2$day, dat2$value, pch=20, col="purple")

Fitting phenology model described by Elmore et al. 2012

The model is a dual logistic curve, with an additional parameter that controls the slope of the line between the two logistic curves. In the plot below, black points are cleaned dominant period values, and red points are fitted values from the model. The two vertical lines are the spring and autumn inflection points of the model.

#Remove NA's from data frame
dat2 <- na.omit(dat2)

#Set starting parameters
m <- c(20, 19, 145, 6, 275, 3, 1)

#Fit non-linear least squares model
mod = nls(value ~ a + (b- g*day)*((1/(1+exp((c-day)/d))) - (1/(1+exp((e-day)/f)))), 
             start = list(a = m[1],b = m[2],c = m[3], d = m[4], e = m[5], f = m[6], g = m[7] ), 
             algorithm="port",control = list(maxiter = 500), data=dat2, weights=dat2$weights^2)

#Plot data, and fitted values
#Units are in tenths of seconds, so divide by 10 to convert to seconds
plot(dat2$day, dat2$value/10, pch=20, xlab="Day of year",ylab="Dominant period (s)")
points(dat2$day, fitted(mod)/10, pch=20, col="red", lwd=4)
abline(v=c(coef(mod)[3],coef(mod)[5]))