FDA tutorial

• fda is a primitive package in R for dealing with functional data

#installing the package
if(!require(fda.usc)){
  install.packages("fda.usc", repos = "http://cran.us.r-project.org")
}

## Loading required package: fda.usc

## Warning: package 'fda.usc' was built under R version 3.2.4

## Loading required package: fda

## Warning: package 'fda' was built under R version 3.2.4

## Loading required package: splines
## Loading required package: Matrix
## 
## Attaching package: 'fda'
## 
## The following object is masked from 'package:graphics':
## 
##     matplot
## 
## Loading required package: MASS
## Loading required package: mgcv
## Loading required package: nlme
## This is mgcv 1.8-9. For overview type 'help("mgcv-package")'.
## Loading required package: rpart

#loading the package
library(fda.usc)

Typical Functional Data

data("tecator")
names(tecator)

## [1] "absorp.fdata" "y"

absorp <- tecator$absorp.fdata

Spectrometric Functional Data Curves

• Curves with Fat < 20 percent (red lines)
• Curves with Fat >= 20 percent (blue lines)

Fat20 <- ifelse(tecator$y$Fat < 20, 0, 1) * 2 + 2
plot(tecator$absorp.fdata, col = Fat20)

absorp.d1 <- fdata.deriv(absorp, nderiv = 1)
plot(absorp.d1, col = Fat20)

class(absorp.fd <- fdata2fd(absorp, type.basis= "fourier", nbasis= 15))

## [1] "fd"

class(absorp.fdata <- fdata(absorp.fd))

## [1] "fdata"

Smoothing using GCV criterion

data("phoneme")
learn <- phoneme$learn
l <- c(0, 2^seq(-2, 9, length.out = 30))
nb <- seq(7, 31, by = 2)
out0 <- min.basis(learn, lambda = l, numbasis = nb)
out1 <- min.np(learn, type.S = S.NW, par.CV = list(criteria = "GCV"))
out2 <- min.np(learn, type.S = S.LLR, par.CV = list(criteria = "GCV"))
out3 <- min.np(learn, type.S = S.KNN, h = 3:25, Ker = Ker.unif)

Exploring functional data

data("poblenou")
nox <- poblenou$nox
dd <- as.integer(poblenou$df$day.week)
working <- poblenou$nox[poblenou$df$day.festive == 0 & dd < 6]
nonworking <- poblenou$nox[poblenou$df$day.festive == 1 | dd > 5]

Plots

plot(func.mean(working), ylim = c(10, 170),main = "Centrality measures in working days")
legend(x = 11,y = 185, cex = 1, box.col = "white", lty = 1:5, col = 1:5,legend = c("mean", "trim.mode", "trim.RP", "median.mode", "median.RP"))
lines(func.trim.mode(working, trim = 0.15), col = 2, lty =2)
lines(func.trim.RP(working, trim = 0.15),col = 3,lty = 3)
lines(func.med.mode(working, trim = 0.15),col = 4,lty = 4)
lines(func.med.RP(working, trim = 0.15),col = 5, lty = 5)