feeding.png

While on paternity leave I had an opportunity to test out periodic smoothing splines (within the framework of generalized additive models) on an interesting time-series-- an infant's feeding schedule.

load / format data and fit GAMs

library(plotrix)
library(mgcv)
library(plyr)
# read and format data
x <- read.csv('feeding.csv'stringsAsFactors=FALSE)
x$datetime <- as.POSIXct(strptime(x$datetimeformat='%Y-%m-%d %H%M'))
x$type <- factor(x$type)
x$d <- c(NAas.numeric(diff(x$datetime))) / 60
x$hour.fraction <- as.numeric(format(x$datetime"%H")) + (as.numeric(format(x$datetime"%M")) / 60)
x$hour <- as.numeric(format(round.POSIXt(x$datetime'hours')"%H"))
x$date <- as.Date(x$datetime)
# get the range of our data as a POSIXct object, rounded to hours
r <- as.POSIXct(round(range(x$datetime[-1]na.rm=TRUE)"hours"))
# generate sequences that we will use later
r.seq <- seq(r[1]r[2]by="12 hours") # for the x-axis of fig. 1
p.seq <- seq(r[1]r[2]by="1 hours") # used for model predictions in fig. 1
# fit a GAM to the time-series
l.ts <- gam(d ~ s(as.numeric(datetime))data=x)
# fit a GAM to feeding interval as smooth, periodic function of hour
l <- gam(d ~ s(hourbs='cc')data=x)
# generate predictions from our time-series model
d.ts <- data.frame(datetime=p.seq)
p.ts <- predict(l.tsd.ts)
p.ts <- data.frame(d.tsfit=as.numeric(p.ts))
# generate predictions from our hourly model
d <- data.frame(hour=seq(023length.out=100))
p <- predict(ldse.fit=TRUE)
p <- data.frame(dfit=as.numeric(p$fit)se.fit=as.numeric(p$se.fit))
# estimate 95% CI from standard error
p$upper <- p$fit + 1.96*p$se.fit
p$lower <- p$fit - 1.96*p$se.fit
# generate hourly predictions for use in fig. 1
p.ts.hourly <- predict(ldata.frame(hour=as.numeric(format(p.seq"%H"))))
# combine time-series model with hourly model predictions for fig. 1
p.ts.hourly.adjusted <- p.ts.hourly - mean(x$dna.rm=TRUE) + p.ts$fit
X

plot

# setup plot layout
layout(matrix(c(1,1,1,2,3,4)nrow=2ncol=3byrow=TRUE)widths=c(1111))
par(mar=c(3,4.5,1,0)cex.axis=0.6cex.lab=0.6)
# fig 1
plot(d ~ datetimedata=xtype='b'axes=FALSExlab=''ylab=''pch=c(1,16)[as.numeric(x$type)]cex=0.75col='RoyalBlue'lwd=1.5)
axis(2cex.axis=0.75line=-0.5las=1at=seq(0.550.5))
mtext('feeding interval (hours)'side=2cex=0.8font=2line=2)
lines(p.seqp.ts.hourly.adjustedlty=2)
lines(p.tslty=3)
axis.POSIXct(side=1at=r.seqcex.axis=0.75format="%m/%d\n%H:%M")
grid()
legend('topright'lty=c(NANA32)pch=c(116NANA)legend=c('breast milk''formula''trend''trend+model')bty='n'cex=0.8col=c('RoyalBlue''RoyalBlue''black''black')horiz=TRUE)
# fig 2
plot(d ~ hour.fractiondata=xaxes=FALSExlab=''ylab=''xlim=c(-1,24)cex=0.75col=rgb(0,0,0alpha=0.75))
lines(p$hourp$fitlwd=2col='RoyalBlue')
lines(p$hourp$lowerlty=2col='RoyalBlue')
lines(p$hourp$upperlty=2col='RoyalBlue')
axis(2cex.axis=0.75line=-0.5las=1at=seq(0.550.5))
mtext('feeding interval (hours)'side=2cex=0.8font=2line=2)
axis(1cex.axis=0.75at=seq(024by=4))
grid()
# fig 3
polar.plot(lengths=x$dpolar.pos=(x$hour.fraction)*360/23rp.type='s'clockwise=TRUEstart=0labels=0:23label.pos=1:24*360/24radial.lim=c(0,5)point.col=rgb(0,0,0alpha=0.75)cex=0.5)
polar.plot(lengths=p$fitpolar.pos=(p$hour)*360/23rp.type='p'clockwise=TRUEstart=0labels=0:23label.pos=1:24*360/24radial.lim=c(0,5)lwd=2line.col='RoyalBlue'add=TRUE)
polar.plot(lengths=p$lowerpolar.pos=(p$hour)*360/23rp.type='p'clockwise=TRUEstart=0labels=0:23label.pos=1:24*360/24add=TRUElty=2radial.lim=c(0,5)line.col='RoyalBlue')
polar.plot(lengths=p$upperpolar.pos=(p$hour)*360/23rp.type='p'clockwise=TRUEstart=0labels=0:23label.pos=1:24*360/24add=TRUElty=2radial.lim=c(0,5)line.col='RoyalBlue')
# fig 4
boxplot(d ~ hourdata=xhorizontal=TRUEborder=rgb(0,0,0alpha=0.75)axes=FALSEboxwex=0.5)
lines(p$fitp$hour+1lwd=2col='RoyalBlue')
lines(p$lowerp$hour+1lty=2col='RoyalBlue')
lines(p$upperp$hour+1lty=2col='RoyalBlue')
axis(1cex.axis=0.75las=1line=-0.5at=seq(0.550.5))
axis(2cex.axis=0.75las=1at=1:24labels=0:23tick=FALSEline=-1)
stripchart(x$hour+1method='stack'vertical=TRUEaxes=FALSEpch='|'cex=0.5add=TRUEat=5.75offset=-0.25)
X

estimate quantiles / ECDF

# estimate hourly quantiles of feeding intervals
q <- ddply(x'hour'function(i) {
        round(quantile(i$dna.rm=TRUE)2)
        })
# remove missing data and print table
q <- na.omit(q)
print(q)
# estimate ECDF by hour
l.ecdf <- list()
for(i in 0:23) {
        l.ecdf[[i+1]] <- ecdf(x$d[x$hour == i])
        }
# plot the ECDF for feeding intervals between 3-4 am
plot(l.ecdf[[3]]verticals=TRUE)
X

Attachment: feeding.csv