Whelp Study
Version 2, August 2014
This is a worksheet for exploring activity metrics that can be drawn from whelp activity as measured by the Fitbit device.
The Fitbit is a very small wearable digital accelerometer worn by the pups unobtrusively on their collars. The Fitbit is designed for use by humans to measure physical activity, and promoted as a fitness aid. Here we are adapting it for use in an animal study, so that what are measured as “steps” in terms of human movement may not translate directly to puppy steps, but have been observed to be a good measure of relative activity for comparing a sample population of puppies.
Per-minute step counts from each Fitbit device are loaded wirelessly to an encrypted data store on a personal computer, then uploaded to the Fitbit website. The Fitbit company has generously provisioned us with a research-only Web account and Web Application Programming Interface (Web API) to retrieve this time series of minute by minute data collected by the Fitbit devices. The data is fetched by our Java program using the Web API, and loaded into a relational database. We are now using the statistical language R with some add-on packages to investigate activity metrics that can be calculated using empirical data collected by the Fitbit devices.
Connect to the Fitbit DB:
library("RODBC")
## Warning: package 'RODBC' was built under R version 2.15.3
conn <- odbcConnect("fitbit")
Get a summary of “steps” for one day of activity
sql1 <- "select dd.DateStr, dd.FBDevice, dd.PupName, SUM(s.Steps)"
sql1 <- paste(sql1, "from DayDevice dd, Steps s ")
sql1 <- paste(sql1, "where dd.DateStr = s.DateStr and dd.FBDevice = s.FBDevice ")
sql1 <- paste(sql1, "and dd.PupName = s.PupName and dd.FullDayData = 1 ")
sql1 <- paste(sql1, "and dd.PupName NOT LIKE 'CONTROL%' ") # exclude FITBIT measurement CONTROL tests
sql1 <- paste(sql1, "and dd.DateStr IN") # include only 1st full day of data
sql1 <- paste(sql1, " (select MIN(DateStr) from DayDevice ")
sql1 <- paste(sql1, " where FullDayData=1 group by PupName)")
sql1 <- paste(sql1, "group by dd.DateStr, dd.FBDevice, dd.PupName")
sql1 <- paste(sql1, "order by 1,2,3")
sum <- sqlQuery(conn, paste(sql1))
As of this run date we have 10 whelp/days of activity being measured.
Activity by whelp/day:
barplot(sum[, 4], col = "lightblue", ylim = c(0, 15000), ylab = "Fitbit 'steps'",
names.arg = as.character(sum[, 3]), las = 3)
mtext("Whelp name", side = 1, line = 4)
mtext("Activity per whelp/day", side = 3, line = 1)
# A simpler plot, same info: plot(sum[,3], sum[,4], type='h', col='red',
# ylim=c(0,15000), ylab='One day activity', xlab='', las=3) mtext('Whelp
# name', side=1, line=4)
Samples are taken two at a time, simultaneously, using two fitbit devices: Fitbit1 and Fitbit2. The samples are taken this way because we test two pups at a time: siblings from a litter whose names (by convention of the breeding kennel) begin with the same letter. We noticed after the first eight samples (four pairs of pups) that Fitbit1 was registering a slightly higher number of steps than FITBIT2, and suspected a measurement bias. Here are the raw step data from those first eight samples, showing the higher measurements pairwise by Fitbit device 1 over device 2:
plot(sum[1:8, 4], type = "b", ylim = c(0, 20000), ylab = "Fitbit 'steps'", xlab = "",
xaxt = "n")
mtext("Fitbit Device Number", side = 1, line = 2.5)
mtext("A measurement bias by Device?", side = 3, line = 1)
axis(1, at = c(1:length(sum[, 2])), labels = substring(as.character(sum[, 2]),
7, 7))
To explore this we set up a comparative test to detect measurement bias between the two Fitbit devices. The test is simply to record a full day of activity during which both devices are worn by the same subject. In a world with zero measurement bias, the devices would have precisely the same minute-by-minute step count. However we did indeed find differences in measurement as illustrated here:
For CONTROL pairs of Fitbit1 and Fitbit2 data, get 24 hours detail “steps” for each day
sqlCtrl <- "select dd.DateStr, dd.FBDevice, dd.PupName, s.TimeStr, s.Steps"
sqlCtrl <- paste(sqlCtrl, "from DayDevice dd, Steps s ")
sqlCtrl <- paste(sqlCtrl, "where dd.DateStr = s.DateStr and dd.FBDevice = s.FBDevice ")
sqlCtrl <- paste(sqlCtrl, "and dd.PupName = s.PupName and dd.FullDayData = 1 ")
sqlCtrl <- paste(sqlCtrl, "and dd.PupName LIKE 'CONTROL%' ") # FITBIT measurement CONTROL tests ONLY
sqlCtrl <- paste(sqlCtrl, "and dd.DateStr IN") # include only 1st full day of data
sqlCtrl <- paste(sqlCtrl, " (select MIN(DateStr) from DayDevice ")
sqlCtrl <- paste(sqlCtrl, " where FullDayData=1 group by PupName)")
sqlCtrl <- paste(sqlCtrl, "order by dd.DateStr, dd.FBDevice, dd.PupName, s.TimeStr")
dtlCtrl <- sqlQuery(conn, paste(sqlCtrl))
# Transform to n cases x 1440 steps matrix
adsCtrl <- as.matrix(xtabs(Steps ~ PupName + TimeStr, dtlCtrl))
A Measurement Contingency Table
Browsing a contingency table of activity vs no activity, we see that both devices had some periods during which they measured activity when the other device didn't measure any. One reassuring point we see here is that periods measuring zero activity (both devices showing FALSE when pups are at rest) are generally in agreement between the two devices.
library(xtable)
## Warning: package 'xtable' was built under R version 2.15.3
contable <- as.data.frame(ftable(c(adsCtrl[1, ] > 0), adsCtrl[2, ] > 0))
contable$ratio = signif(contable[, 3]/length(adsCtrl[1, ]), 2)
colnames(contable) <- c("Fitbit 1", "Fitbit 2", " Minutes (total 1440)", " Ratio")
rownames(contable) <- paste(1:4, "Steps recorded > 0")
print(xtable(contable), type = "html", html.table.attributes = "border=1, bgcolor='lightblue'")
| Fitbit 1 | Fitbit 2 | Minutes (total 1440) | Ratio | |
|---|---|---|---|---|
| 1 Steps recorded > 0 | FALSE | FALSE | 1240 | 0.86 |
| 2 Steps recorded > 0 | TRUE | FALSE | 28 | 0.02 |
| 3 Steps recorded > 0 | FALSE | TRUE | 11 | 0.01 |
| 4 Steps recorded > 0 | TRUE | TRUE | 161 | 0.11 |
Boxplots of Fitbit1 vs Fitbit2 measurement
The boxplots show some variation in the distribution of measurement between the devices. The overlap of the notches around the means of the box plots illustrate that the medians do not differ significantly (p = 0.05) as calculated by the boxplot notch parameter (Chambers et al., 1983, p. 62).
par(mfrow = c(1, 2))
# Get distributions of minutes in which at least one device recorded steps
# > 0
ctrl1 <- adsCtrl[1, adsCtrl[1, ] > 0 | adsCtrl[2, ] > 0]
ctrl2 <- adsCtrl[2, adsCtrl[1, ] > 0 | adsCtrl[2, ] > 0]
boxplot(ctrl1, ctrl2, width = c(rep(5, 2)), notch = T, names = c("Fitbit1",
"Fitbit2"))
Another visual to compare the variation in measurement
In this visual we see how the variances in measurement between Fitbit1 and Fitbit2 are spread across the histogram quantiles of measured 'steps' per minute. Across the spectrum, sometimes Fitbit1 registered more steps, and about as often Fitbit2 registered more steps. This adds intuitive confidence that the measurement variance between the two devices is by chance rather than a systematic measurement bias, such as a difference in senstivity of the accelerometers inside the two Fitbits.
par(mfrow = c(1, 1))
colors = c(rgb(1, 0, 0, 0.5), rgb(0, 0, 1, 0.5), rgb(0.5, 0, 1, 0.8))
hist(ctrl1, col = colors[1], main = "Fitbit1 vs Fitbit2 Measurement Variation",
ylab = "Fitbit 'steps' per minute", xlab = "Minutes in which at least one device recorded steps > 0")
hist(ctrl2, col = colors[2], add = T)
legend("topright", c("Fitbit1", "Fitbit2", "Overlap"), col = colors, pch = 15,
cex = 1.2, inset = 0.02)
SHOULD we adjust raw data for measurement variance? Is there a significant bias?
Test for normality Although we can see this in the visuals, we run the Shapiro test which gives a very high confidence for both Fitbit measurement test data sets that we can reject the NULL hypothesis of a normal distribution.
shapiro.test(ctrl1)
##
## Shapiro-Wilk normality test
##
## data: ctrl1
## W = 0.8488, p-value = 3.75e-13
shapiro.test(ctrl2)
##
## Shapiro-Wilk normality test
##
## data: ctrl2
## W = 0.8513, p-value = 4.983e-13
Wilcoxon Rank test The samples are independent since they are from two separate devices. Each measurement in a sample is independent rather than part of a time series, since each minute of measurement by a Fitbit device has no dependency on measurements in prior minutes.
Given these characteristics of the samples we run the Wilcoxon Rank Sum Test for unpaired samples to test the hypothesis that the distributions differ significantly with a confidence threshold of 95%.
wil = wilcox.test(ctrl1, ctrl2, conf.int = T, conf.level = 0.95)
wil
##
## Wilcoxon rank sum test with continuity correction
##
## data: ctrl1 and ctrl2
## W = 21246, p-value = 0.281
## alternative hypothesis: true location shift is not equal to 0
## 95 percent confidence interval:
## -1 4
## sample estimates:
## difference in location
## 1
With a p-value of 0.281 we cannot reject the null hypothesis that the true Fitbit1 and Fitbit2 distributions are the same.
Get 24 hours detail “steps” for each whelp/day (1440 minutes/day)
sql2 <- "select dd.DateStr, dd.FBDevice, dd.PupName, s.TimeStr, s.Steps"
sql2 <- paste(sql2, "from DayDevice dd, Steps s ")
sql2 <- paste(sql2, "where dd.DateStr = s.DateStr and dd.FBDevice = s.FBDevice ")
sql2 <- paste(sql2, "and dd.PupName = s.PupName and dd.FullDayData = 1 ")
sql2 <- paste(sql2, "and dd.PupName NOT LIKE 'CONTROL%' ") # exclude FITBIT measurement CONTROL tests
sql2 <- paste(sql2, "and dd.DateStr IN") # include only 1st full day of data
sql2 <- paste(sql2, " (select MIN(DateStr) from DayDevice ")
sql2 <- paste(sql2, " where FullDayData=1 group by PupName)")
sql2 <- paste(sql2, "order by dd.DateStr, dd.FBDevice, dd.PupName, s.TimeStr")
dtl <- sqlQuery(conn, paste(sql2))
# Transform to n cases x 1440 steps matrix
ads <- as.matrix(xtabs(Steps ~ PupName + TimeStr, dtl))
plot(0, 0, xlim = c(0, 1440), ylim = c(0, dim(ads)[1]), type = "n", xlab = "",
ylab = "")
rasterImage(as.raster(abs(ads[, 1:1440] - 255), max = 255), 0, 0, 1440, dim(ads)[1],
interpolate = FALSE)
mtext("1440 minutes from 12AM to 12PM", side = 1, line = 2.5)
mtext("Whelp number", side = 2, line = 2)
mtext("Activity Map", side = 3, line = 2)
mtext("Darker = more activity", side = 3, line = 1)
Defined as: Number of Naps per day, where a Nap is: a period of time between 6am and 10pm having at least 10 minutes of no steps counted.
library(TTR)
## Warning: package 'TTR' was built under R version 2.15.3
## Loading required package: xts
## Warning: package 'xts' was built under R version 2.15.3
## Loading required package: zoo
## Warning: package 'zoo' was built under R version 2.15.3
##
## Attaching package: 'zoo'
##
## The following object(s) are masked from 'package:base':
##
## as.Date, as.Date.numeric
source("C:\\R\\fitbit\\fitbit_lib.R")
naps <- rep(0, length(ads[, 1]))
for (i in 1:length(naps)) {
naps[i] <- napCount(ads[i, ], 10, 360, 1320) # params: 10=10minutes, 360=6am, 1320=10pm
}
barplot(naps, col = "lightblue", ylab = "Naps per Day", names.arg = as.character(sum[,
3]), las = 3)
mtext("Nap Count", side = 3, line = 1)
A better approach: build a ragged list of NapData vectors, from which more stats can calculated with a single set of assumptions about how “Nap” is defined.
Daily Nap Sleep Time - defined as: Sum of daily Nap times in hours.
Average Nap Time - defined as: Average of daily Nap times in minutes.
Standard Deviation of Nap Times - defined as: Standard deviation of daily Nap times in minutes.
napdat <- list()
for (i in 1:length(ads[, 1])) {
# params: 10=10minutes, 360=6am, 1320=10pm
napdat[[i]] <- napData(ads[i, ], 10, 360, 1320)
}
# Put nap stats into dataframe
napDf <- data.frame(ct = numeric(0), sum = numeric(0), mean = numeric(0), sd = numeric(0))
for (i in 1:length(napdat)) {
napDf[i, "ct"] = length(napdat[[i]])
napDf[i, "sum"] = sum(napdat[[i]])
napDf[i, "mean"] = mean(napdat[[i]])
napDf[i, "sd"] = sd(napdat[[i]])
}
# Bar plots
par(mfrow = c(2, 2))
barplot(napDf$ct, col = "lightblue", ylab = "Naps per Day", names.arg = as.character(sum[,
3]), las = 3)
mtext("Nap Count", side = 3, line = 1)
barplot(napDf$sum/60, col = "lightgreen", ylab = "Sum of Nap times in hours",
names.arg = as.character(sum[, 3]), las = 3)
mtext("Daily Nap Sleep Time", side = 3, line = 1)
barplot(napDf$mean, col = "pink", ylab = "Avg Nap time in minutes", names.arg = as.character(sum[,
3]), las = 3)
mtext("Average Nap Time", side = 3, line = 1)
barplot(napDf$sd, col = "orange", ylab = "Std Deviation in minutes", names.arg = as.character(sum[,
3]), las = 3)
mtext("Standard Deviation of Nap Times", side = 3, line = 1)
Metrics based on Nap Cycle Time were suggested by Clare's mentor - a possible interpretation:
In the previous section the napData() function returns a vector of nap lengths in minutes, with one element for each “Nap” as defined by the function params. For an interpretation of Nap Cycle Time, we have a new function napTimeElapsedSincePrevious() which will also return a vector of one element per nap: the number of minutes elapsed since the start of the previous nap. For the first element the value will be zero since there's no previous nap to compare. This way we can measure consistency of nap cycles, and optionally combine it with length of naps from the napData() function to derive other metrics.
Average Nap Cycle Time - defined as: Average of daily Nap cycle times in minutes.
Standard Deviation of Nap Cycle Times - defined as: Standard deviation of daily Nap cycle times in minutes.
napCycDat <- list()
for (i in 1:length(ads[, 1])) {
# params: 10=10minutes, 360=6am, 1320=10pm
napCycDat[[i]] <- napTimeElapsedSincePrevious(ads[i, ], 10, 360, 1320)
}
# Validate that napdat and napCycDat are same size
if (sum(summary(napdat[1:length(napdat)])[, 1] == summary(napCycDat[1:length(napdat)])[,
1]) != length(napdat)) +stop("An error occurred. Nap counts in napdat and napCycDat don't match!")
# Put nap cycle time stats into dataframe
napCycDf <- data.frame(ct = numeric(0), mean = numeric(0), sd = numeric(0))
for (i in 1:length(napCycDat)) {
napCycDf[i, "ct"] = length(napCycDat[[i]])
napCycDf[i, "mean"] = mean(napCycDat[[i]][2:length(napCycDat[[i]])]) # exclude first element which is set to zero
napCycDf[i, "sd"] = sd(napCycDat[[i]][2:length(napCycDat[[i]])]) # exclude first element which is set to zero
napCycDf[i, "median"] = median(napCycDat[[i]][2:length(napCycDat[[i]])]) # exclude first element which is set to zero
}
# Bar plots
par(mfrow = c(2, 2))
barplot(napCycDf$ct, col = "lightblue", ylab = "Naps per Day", names.arg = as.character(sum[,
3]), las = 3)
mtext("Nap Count", side = 3, line = 1)
barplot(napCycDf$mean, col = "pink", ylab = "Avg Nap cycle time in minutes",
names.arg = as.character(sum[, 3]), las = 3)
mtext("Average Nap Cycle Time", side = 3, line = 1)
barplot(napCycDf$sd, col = "orange", ylab = "Std Deviation in minutes", names.arg = as.character(sum[,
3]), las = 3)
mtext("Standard Deviation of Nap Cycle Times", side = 3, line = 1)
barplot(napCycDf$mean, col = "lightgreen", ylab = "Median Nap cycle time in minutes",
names.arg = as.character(sum[, 3]), las = 3)
mtext("Median Nap Cycle Time", side = 3, line = 1)
Since the pups we are able to monitor usually come to us in pairs of siblings, we could also look at times when a pup is active but the sibling is resting, and vice versa - that might be interesting too. Even when they're both active, it could be that one is consistently more active/feisty than the other. That way a confounding factor (the influence of the sibling degrading the independence of observations) could be treated as an asset - another behavior metric to study!