Fitbit band data cluster

序言
Data and EDA
Clustering of the data in Aug
聚类应用

关于统计建模的理解: 统计建模是复杂规律的简化描述，把不能理解的因素理解成随机误差。总可以找到一种方式把事物变化理解成带约束最优化某个量。很多 regularized regression 就是一个约束变量个数的直观优化，带约束的优化可以转换为其他形式的不带约束的不直观优化。这种数学转化的 trade off 是获得了容易计算的特性，损失了可解释性。
统计建模是复杂规律的猜测，模型的可解释性支持了猜测的合理性。很多模型实质是用关于已知信息可测的投影预测未知变量，在特定模型假定下用已知的数据信息来给出预测，比如常见的线性回归，分位数回归。

序言

本数据是典型的函数数据，函数数据的聚类分析有若干中算法，本文将尝试用函数数据聚类分析包去对fitbit band data 进行聚类。

预期结果

周末和平常的心跳模式应该会有较大区别，可能会聚类成两类
周末心跳模式变化大，他们可能无法成为一类，如果分析能找出所有的工作日也挺好。

文献资料

R package Funclustering, 总体来说这一种类似EM算法
Functional data clustering: a survey, Clustering techniques for functional data are reviewed. Four groups of clustering algorithms for functional data are proposed.
fda.usc R package 中 kmeans.fd 函数可以实现K-means 聚类。
其他算法包括 funFEM, Classification of Functional Data: Unsupervised Curve Clustering When Curves are Misaligned, Distance for Functional Data Clustering Based on Smoothing Parameter Commutation, Bayesian Clustering of Functional Data Using Local Features

Data and EDA

疑问？

数据读入时，第一到第十二天的数据文件名称很好，其他文件几天的数据并成一天了，需要清理一下？铁哥你有没有按天清理好的数据呢？好把，数据清理的工作我去弄弄吧
Aug 28 和 Aug 31 日的数据好像有缺失

Data cleaning

我将会把每天的心跳输出整理成一个单独的文件,方便读取和处理。我们得到整个Aug的数据，数据清理过程的代码可打开源代码文件查看。为了保持美观，就不显示清理的代码了。

\(\\\)

数据读入

数据是fitbit band 记录的八月每天心跳，先读入数据成一个list
然后将这个数据存储成 .Rdata文件

url1 <- "/home/ghy/Dropbox/User1/cleaned_data/detailed-data-fitbit-aug"
url3 <- "_16-HR.csv"
day <- list()
for (i in 1:31) {
  url2 <- paste0(i,"_",i)
  url <- paste0(url1, url2, url3)
  day[[i]] = read.csv(url)
}

save(day, file = "/home/ghy/Dropbox/User1/cleaned_data/day.Rdata")
save(day, file = "/home/ghy/Dropbox/User1/gonghy_R_data_analysis/heartrate_app/day.Rdata")

我们画图看看每天心率变化。

var.rename <- function(x){
  names(x) <- c("time", "rate")
  x
}

# 画图函数

plot.data <- function(test){
  names(test) <- c("time", "rate")
  test.time <-strptime(test$time, format = "%Y-%m-%d %H:%M:%S")
  qplot(test.time, test$rate,xlab="时间",ylab="心率",main="心率随时间变化图", geom = "line")
}

# Aug 17 的心率图

plot.data(day[[17]])

更进一步，我们研究的是八月份的数据，我们希望有个 App 可以查看每天的心率图。

# 如果电脑安装了Rstudio直接运行如下代码
library(shiny)
shinyAppDir("/home/ghy/Dropbox/User1/gonghy_R_data_analysis/heartrate_app/")

此app的运行需要Rstudio。如果没有此环境，可以访问网址 https://heyanggong.shinyapps.io/heartrate_app/

异常数据: 可以发现 Aug 28, Aug 31 有某个时间段，完全没有数据，我们需要对该天的数据进行相应处理。

Clustering of the data in Aug

分析的基本步骤：

转为为函数数据类型
选取合适的方法进行聚类

这里我们用函数数据分析的方式来处理数据, 下面把数据转换为函数数据类型。

我们得到了函数对象的list rate.fd, 然后用R package 将它进行聚类。

方法一，用 Funclustering 包，代码繁杂，可以打开源文件查看源代码。
方法二，用fda.usc包，K-means方法进行聚类

set.seed(100)
library(fda.usc)
res <- kmeans.fd(fd, ncl = 2)

res$cluster

##  [1] 1 2 1 1 1 1 1 1 2 1 2 1 2 1 1 1 1 1 2 1 1 1 1 1 1 1 1

res <- kmeans.fd(fd, ncl = 3)

res$cluster

##  [1] 1 3 1 2 2 2 1 1 3 1 3 1 3 1 1 1 1 1 3 2 3 2 1 2 1 1 2

第二种方法简单方便。聚类分析的基本方法至此完成。

聚类应用

我们可以根据心率图来预测是周几？

Research Question: 我们将所有的周三与周五的数据进行放在一起，然后聚类分析

B样条基

minutebasis <- create.bspline.basis(rangeval = c(0, 1440), norder = 4, nbasis = 20)

wed <- list()
wed[[1]] <- day[[3]]
wed[[2]] <- day[[10]]
wed[[3]] <- day[[17]]
wed[[4]] <- day[[24]]

sat <- list()
sat[[1]] <- day[[5]]
sat[[2]] <- day[[12]]
sat[[3]] <- day[[19]]
sat[[4]] <- day[[26]]

test <- list()
coef <- c()

for (i in 1:4){
  test[[i]] <- with(wed[[i]] %>% var.rename, smooth.basis(time.to.minute(time), rate, minutebasis)$fd)
  coef <- cbind(coef, coef(test[[i]]))
}

for (i in 5:8){
  test[[i]] <- with(sat[[i-4]] %>% var.rename, smooth.basis(time.to.minute(time), rate, minutebasis)$fd)
  coef <- cbind(coef, coef(test[[i]]))
}

fd <- fd(coef, basisobj = minutebasis)

set.seed(10)
res <- kmeans.fd(fd, ncl = 2)

res$cluster

## [1] 1 2 2 2 2 1 1 2

Fourier 基

minutebasis <- create.fourier.basis(rangeval = c(0, 1440), nbasis = 10)

wed <- list()
wed[[1]] <- day[[3]]
wed[[2]] <- day[[10]]
wed[[3]] <- day[[17]]
wed[[4]] <- day[[24]]

sat <- list()
sat[[1]] <- day[[5]]
sat[[2]] <- day[[12]]
sat[[3]] <- day[[19]]
sat[[4]] <- day[[26]]

test <- list()
coef <- c()

for (i in 1:4){
  test[[i]] <- with(wed[[i]] %>% var.rename, smooth.basis(time.to.minute(time), rate, minutebasis)$fd)
  coef <- cbind(coef, coef(test[[i]]))
}

for (i in 5:8){
  test[[i]] <- with(sat[[i-4]] %>% var.rename, smooth.basis(time.to.minute(time), rate, minutebasis)$fd)
  coef <- cbind(coef, coef(test[[i]]))
}

fd <- fd(coef, basisobj = minutebasis)

set.seed(100)
res <- kmeans.fd(fd, ncl = 2)
res$cluster


par(mfrow = c(1,1))
plot(fd)

此方法下预测分类的正确率为 5/8。正确比较低，可能原因有

数据不够，每个类才4个样本
需要选取合适的基函数
本身心跳数据并没有过多的关于周几的信息

可以尝试如下途径改进结果

可能此女生活并不是特别规律，换成其他规律生活的人可能准确率大大提高
此种聚类分析是把函数图像相似的聚类，可以换成其他更加合理的相似性度量
既然一个解释变量心跳不够，我们可以引入更多解释变量。

我再试试其他的。

minutebasis <- create.fourier.basis(rangeval = c(0, 1440), nbasis = 10)

wed <- list()
wed[[1]] <- day[[1]]
wed[[2]] <- day[[8]]
wed[[3]] <- day[[15]]
wed[[4]] <- day[[22]]

sat <- list()
sat[[1]] <- day[[5]]
sat[[2]] <- day[[12]]
sat[[3]] <- day[[19]]
sat[[4]] <- day[[26]]

test <- list()
coef <- c()

for (i in 1:4){
  test[[i]] <- with(wed[[i]] %>% var.rename, smooth.basis(time.to.minute(time), rate, minutebasis)$fd)
  coef <- cbind(coef, coef(test[[i]]))
}

for (i in 5:8){
  test[[i]] <- with(sat[[i-4]] %>% var.rename, smooth.basis(time.to.minute(time), rate, minutebasis)$fd)
  coef <- cbind(coef, coef(test[[i]]))
}

fd <- fd(coef, basisobj = minutebasis);plot(fd)

set.seed(10000)
res <- kmeans.fd(fd, ncl = 2)
res$cluster


par(mfrow = c(1,1))
plot(fd)

Fitbit band data cluster

Heyang_Gong

Sat Oct 29 11:02:09 2016

序言

Data and EDA

Clustering of the data in Aug

聚类应用