---
title: "'Accuracy and precision of energy expenditure, heart rate, and steps measured by combined-sensing Fitbits: Systematic Review and Meta-Analyses"
author: "Guillaume Chevance"
date: "last update July 2021"
output:
pdf_document: default
html_document: default
description: this code runs several meta-analyses of Bland-Altman studies
protocol & pre-registration: https://osf.io/preprints/sportrxiv/tgk7a/
contact: guillaume.chevance@isglobal.org
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
#Import library
library(dplyr)
library(outliers)
```
Import data
```{r}
#rm(list=ls())
data <- read.csv("data_r_4.csv", sep=";")
#library(readxl)
#df <- read_excel("Individual subjects raw data.xlsx")
```
1. Calculation of the mean bias and relative sd from means, sds and correlations
```{r}
selectionCriterion1 <- c(1:13, 33:36, 76:85, 106:107, 118:119, 127, 142, 154, 167, 262:266, 287:288)
data2 <- data[selectionCriterion1,]
data2$BA_bias <- data2$Mean_Fitbit - data2$Mean_Criterion
data2$BA_SD <- data2$SD_Fitbit^2 + data2$SD_Criterion^2 - 2*data2$Correlation*data2$SD_Fitbit*data2$SD_Criterion
data2$BA_SD <- sqrt(data2$BA_SD)
# For studies that have Two mean and two SD
data[selectionCriterion1,] <- data2
```
2. Calculation of the bias's sd from LOA (Limits of Agreements)
```{r}
selectionCriterion2<- c(37:38, 58:59, 62, 64:68, 104:105, 155:160, 210:218, 297:302)
data3 <- data[selectionCriterion2,]
data3$BA_SD <- (data3$LOA_U - data3$BA_bias)/1.96
# For studies that have bias and LOA
# We have two methods to calculate the agreement, according to what we have in our raw data.
data[selectionCriterion2,] <- data3
```
3. Rescale from absolute to relative values for steps and EE (energy expenditure, depedent) + remove NA for "duration_minutes"
(i.e., divide mean bias and relative sd by the time of the task/protocol in minutes)
```{r}
summary(data$BA_bias)
selectionCriterion3 <- !(data$Outcome=="HR") & !is.na(data$Duration_minutes) & !(data$relative_absolute=="relative")
data4 <- data[selectionCriterion3,]
data4$BA_bias <- data4$BA_bias / data4$Duration_minutes
# standirize procesdure to put all the bias in the "minute scale"
data4$BA_SD <- data4$BA_SD / data4$Duration_minutes
data[selectionCriterion3,] <- data4
data <- data[!is.na(data$BA_bias) & !is.na(data$BA_SD),]
#Until here we are computing the effect size of the MA
## Now we have all variables of interest, and we do teh analysis
```
Descriptives
```{r}
boxplot(BA_bias~Outcome, dat =data)
# Boxplot with Outcome as categories
n_distinct(data$Study)
table(data$Outcome)
unique(data$N_study) # 20 studies 3 172 subjects #PT
unique(data$ï..Study) # 20 studies 3 172 subjects #PT
table(data$ï..Study) # 20 studies 3 172 subjects #PT
```
Meta-analyses
```{r}
#to avoid N (number of data points) > n (number of observation)
# This statement is helping to re-define the unit of analysis.
# This should be revised to understand the step well.
#
# N_study number of subjects
# N_analysis number of observation (multiple measurements for a subject, always in HR and repeated measure)
data$N_study <- ifelse((data$N_study > data$N_analysis) & (data$Repeated_measures=="no"), data$N_analysis, data$N_study)
#data$N_study1 <- ifelse(data$N_study > data$N_analysis & data$Repeated_measures=="no", data$N_analysis, data$N_study)
#data$N_study3 <- ifelse((data$N_study > data$N_analysis) , data$N_analysis, data$N_study)
#Only for the row 127 and 128
#data$N_study-data$N_study2
#data$N_study3-data$N_study2
#data$N_study-data$N_study3
#This is for missingness: (data$N_study > data$N_analysis) CFM
#
#
#
#unique(data$N_analysis)
dfdf<-cbind(N_study=data$N_study,N_analysis= data$N_analysis) # PT
# call the BA_function to run the MA
source("BA_functions.R")
#MA loop
runMeta <- function(data){
data$needed <- ifelse(data$Repeated_measures == "no", 1, 0)
#------------------------------------------------
# Create imputed datasets by rho value
#------------------------------------------------
bias <- data$BA_bias
V_bias <- data$BA_SD^2/data$N_analysis
logs2 <- log(data$BA_SD^2) # Log variance of the Bias. Accroding to Tripen the transformation was applied because the pool variance is not normal, when the SS is small.
V_logs2 <- 2/(data$N_analysis - 1)
N <- data$N_analysis
n <- data$N_study
needed <- data$needed
# scenarios for the correlation for the repeated measures. Which ussually is unknown.
# We should have values base on the experts
new.8 <- rho_imputer(bias,V_bias,logs2,V_logs2,N,n,needed,.8) # Internal funcion inside BA_functions.R
new.5 <- rho_imputer(bias,V_bias,logs2,V_logs2,N,n,needed,.5) # Not used in the study
new1.0 <- rho_imputer(bias,V_bias,logs2,V_logs2,N,n,needed,1)
# At The end they choose the correlation value 0.8, because give similar results for the pool bias estimation.
#---------------------------------------------------
# Do analyses in various ways #This is creating scenarios for the tho assumption in the subject level data.
#---------------------------------------------------
#For only included data
inc <- subset(new.8, needed == "0")
bias <- inc$bias
V_bias <- inc$V_bias
logs2 <- inc$logs2
V_logs2 <- inc$V_logs2
# loa_maker is doing the meata-analysis for the Bias and for the log of variance of bias
#
results_inc <- loa_maker(bias,V_bias,logs2,V_logs2)
#For rho = .8
bias <- new.8$bias
V_bias <- new.8$V_bias
logs2 <- new.8$logs2
V_logs2 <- new.8$V_logs2
results_8 <- loa_maker(bias,V_bias,logs2,V_logs2)
#For rho = .5
#bias <- new.5$bias
#V_bias <- new.5$V_bias
#logs2 <- new.5$logs2
#V_logs2 <- new.5$V_logs2
# results_5 <- loa_maker(bias,V_bias,logs2,V_logs2)
#for rho = 1.0
bias <- new1.0$bias
V_bias <- new1.0$V_bias
logs2 <- new1.0$logs2
V_logs2 <- new1.0$V_logs2
results_1 <- loa_maker(bias,V_bias,logs2,V_logs2)
#-----------------------------------------------------
#Combine into table
#----------------------------------------------------
results <- rbind(results_inc,results_1,results_8)
results <- data.frame(results)
results$rho <- c("na","1.0",".8")
names(results) <- c("studies","bias","sd","tau2","LOA_L","LOA_U","CI_Lm","CI_Um","CI_Lr","CI_Ur","rho")
return(results)
}
```
Main analyses: full + without low quality + without outliers
```{r}
#heart rate
metaHR <- runMeta(data[data$Outcome== "HR",])
metaHR
tapply(data$Total, data$Total, length)
data_quality <- subset(data, Total>"1")
metaHR_quality <- runMeta(data_quality[data$Outcome== "HR",])
metaHR_quality
#energy expenditure
metaEE <- runMeta(data[data$Outcome== "EE",])
metaEE
metaEE_quality <- runMeta(data_quality[data$Outcome== "EE",])
metaEE_quality
#steps ("NA produce" - it is normal here - repeated measures always "no" for steps)
metaSTEPS <- runMeta(data[data$Outcome== "Steps",])
metaSTEPS
metaSTEPS_quality <- runMeta(data_quality[data$Outcome== "Steps",])
metaSTEPS_quality
###outliers detection###
boxplot(BA_bias~Outcome, dat =data)
data_HR <- subset(data, Outcome=="HR")
data_EE <- subset(data, Outcome=="EE")
data_STEPS <- subset(data, Outcome=="Steps")
grubbs.test(data_HR$BA_bias)
grubbs.test(data_EE$BA_bias)
grubbs.test(data_STEPS$BA_bias)
data_HR_outlier <- subset(data_HR, BA_bias>-30.5)
data_EE_outlier <- subset(data_EE, BA_bias<61)
data_STEPS_outlier <- subset(data_STEPS, BA_bias>-53.07)
boxplot(data_EE_outlier$BA_bias)
boxplot(data_HR_outlier$BA_bias)
boxplot(data_STEPS_outlier$BA_bias)
metaHR_outlier <- runMeta(data_HR_outlier)
metaHR_outlier
metaEE_outlier <- runMeta(data_EE_outlier)
metaEE_outlier
metaSTEPS_outlier <- runMeta(data_STEPS_outlier)
metaSTEPS_outlier
```
Subgroup analyses
```{r}
###participants' characteristics
################################
### Maybe we do meta-regression, we need to check for the distirbution of the continuous variables.
### If not Ratko/Marco/Issac should provide a cut-off point
### Experts should provide clinical LOA criterias.
#CHECK WHAT IS DOING THIS FUNCTION
tapply(data$Participants, data$Participants, length)
# Drop Unused Levels from Factors
data_healthy <- droplevels(data[!data$Participants == 'clinical',])
#data_clin <- subset(data, Participants == 'clinical')
metaHR_healthy <- runMeta(data_healthy[data_healthy$Outcome== "HR",])
metaHR_healthy
metaEE_healthy <- runMeta(data_healthy[data_healthy$Outcome== "EE",])
metaEE_healthy
metaSTEPS_healthy <- runMeta(data_healthy[data_healthy$Outcome== "Steps",])
metaSTEPS_healthy
tapply(data$Age_group, data$Age_group, length)
data_young <- droplevels(data[!data$Age_group == 'older',])
data_old <- subset(data, Age_group=="older")
metaHR_young <- runMeta(data_young[data$Outcome== "HR",])
metaHR_young
metaHR_old <- runMeta(data_old[data$Outcome== "HR",])
metaHR_old
metaEE_young <- runMeta(data_young[data$Outcome== "EE",])
metaEE_young
metaEE_old <- runMeta(data_old[data$Outcome== "EE",])
metaEE_old
metaSTEPS_young <- runMeta(data_young[data$Outcome== "Steps",])
metaSTEPS_young
metaSTEPS_old <- runMeta(data_old[data$Outcome== "Steps",])
metaSTEPS_old
### Fitbit devices
##################
################## Subgroup only the other device
##################
tapply(data$Fitbit, data$Fitbit, length)
data_blaze <- subset(data, Fitbit=="Blaze")
data_charge2 <- subset(data, Fitbit=="Charge 2")
data_chargeHR <- subset(data, Fitbit=="Charge HR")
data_Surge <- subset(data, Fitbit=="Surge")
metaHR_blaze <- runMeta(data_blaze[data$Outcome== "HR",])
metaHR_blaze
metaEE_blaze <- runMeta(data_blaze[data$Outcome== "EE",])
metaEE_blaze
metaSTEPS_blaze <- runMeta(data_blaze[data$Outcome== "Steps",])
metaSTEPS_blaze
metaHR_charge2 <- runMeta(data_charge2[data$Outcome== "HR",])
metaHR_charge2
metaEE_charge2 <- runMeta(data_charge2[data$Outcome== "EE",])
metaEE_charge2
metaSTEPS_charge2 <- runMeta(data_charge2[data$Outcome== "Steps",])
metaSTEPS_charge2
metaHR_chargeHR <- runMeta(data_chargeHR[data$Outcome== "HR",])
metaHR_chargeHR
metaEE_chargeHR <- runMeta(data_chargeHR[data$Outcome== "EE",])
metaEE_chargeHR
metaSTEPS_chargeHR <- runMeta(data_chargeHR[data$Outcome== "Steps",])
metaSTEPS_chargeHR
metaHR_Surge <- runMeta(data_Surge[data$Outcome== "HR",])
metaHR_Surge
metaEE_Surge <- runMeta(data_Surge[data$Outcome== "EE",])
metaEE_Surge
metaSTEPS_Surge <- runMeta(data_Surge[data$Outcome== "Steps",])
metaSTEPS_Surge
data_versa <- subset(data, Fitbit=="Versa")
metaHR_versa <- runMeta(data_versa[data$Outcome== "HR",])
metaHR_versa
metaEE_versa <- runMeta(data_versa[data$Outcome== "EE",])
metaEE_versa
metaSTEPS_versa <- runMeta(data_versa[data$Outcome== "Steps",])
metaSTEPS_versa
```
```{r}
### Device placement
####################
tapply(data$Fitbit, data$Fitbit, length)
table(data$Fitbit)
#data_placementOK <- subset(data, Q4_placement==1)
data_placementOK <- subset(data, Q4_multiple.devices==1)
metaHR_place <- runMeta(data_placementOK[data$Outcome== "HR",])
metaHR_place
metaEE_place <- runMeta(data_placementOK[data$Outcome== "EE",])
metaEE_place
metaSTEPS_place <- runMeta(data_placementOK[data$Outcome== "Steps",])
metaSTEPS_place
### Type of activity
####################
tapply(data$Activities, data$Activities, length)
data_cycling <- subset(data, Activities=="Cycling")
metaHR_cy <- runMeta(data_cycling[data$Outcome== "HR",])
metaHR_cy
metaEE_cy <- runMeta(data_cycling[data$Outcome== "EE",])
metaEE_cy
data_daily <- subset(data, Activities=="Daily living")
metaHR_dl <- runMeta(data_daily[data$Outcome== "HR",])
metaHR_dl
metaEE_dl <- runMeta(data_daily[data$Outcome== "EE",])
metaEE_dl
metaSTEPS_dl <- runMeta(data_daily[data$Outcome== "Steps",])
metaSTEPS_dl
data_trdm<- subset(data, Activities=="Trdm")
metaHR_trdm <- runMeta(data_trdm[data$Outcome== "HR",])
metaHR_trdm
metaEE_trdm <- runMeta(data_trdm[data$Outcome== "EE",])
metaEE_trdm
metaSTEPS_trdm<- runMeta(data_trdm[data$Outcome== "Steps",])
metaSTEPS_trdm
data_owalk<- subset(data, Activities=="Walking")
metaHR_owalk <- runMeta(data_owalk[data$Outcome== "HR",])
metaHR_owalk
metaEE_owalk <- runMeta(data_owalk[data$Outcome== "EE",])
metaEE_owalk
metaSTEPS_owalk <- runMeta(data_owalk[data$Outcome== "Steps",])
metaSTEPS_owalk
### Intensity
#############
tapply(data$Intensities, data$Intensities, length)
data_modvig<- subset(data, Intensities=="MV")
metaHR_modvig <- runMeta(data_modvig[data$Outcome== "HR",])
metaHR_modvig
metaEE_modvig <- runMeta(data_modvig[data$Outcome== "EE",])
metaEE_modvig
metaSTEPS_modvig <- runMeta(data_modvig[data$Outcome== "Steps",])
metaSTEPS_modvig
data_light <- subset(data, Intensities=="Light")
metaHR_light <- runMeta(data_light[data$Outcome== "HR",])
metaHR_light
metaEE_light <- runMeta(data_light[data$Outcome== "EE",])
metaEE_light
metaSTEPS_light <- runMeta(data_light[data$Outcome== "Steps",])
metaSTEPS_light
```
```{r}
sessionInfo()
#R version 4.0.3 (2020-10-10)
#Platform: x86_64-apple-darwin17.0 (64-bit)
#Running under: macOS Mojave 10.14.6
#Matrix products: default
#BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
#LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
#Random number generation:
# RNG: Mersenne-Twister
# Normal: Inversion
# Sample: Rounding
#locale:
#[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
#attached base packages:
#[1] stats graphics grDevices utils datasets methods base
#other attached packages:
#[1] stringr_1.4.0 outliers_0.14 dplyr_1.0.7 pracma_2.2.9 plyr_1.8.6 casnet_0.2.0 readxl_1.3.1 rio_0.5.16
#loaded via a namespace (and not attached):
# [1] tidyselect_1.1.0 xfun_0.18 purrr_0.3.4 haven_2.3.1 lattice_0.20-41 colorspace_1.4-1 vctrs_0.3.8 generics_0.0.2
# [9] utf8_1.1.4 blob_1.2.1 rlang_0.4.11 pillar_1.6.1 foreign_0.8-80 glue_1.4.2 DBI_1.1.0 lifecycle_1.0.0
#[17] munsell_0.5.0 gtable_0.3.0 cellranger_1.1.0 zip_2.1.1 knitr_1.30 forcats_0.5.0 curl_4.3 fansi_0.4.1
#[25] Rcpp_1.0.6 scales_1.1.1 farver_2.0.3 ggplot2_3.3.2 hms_0.5.3 digest_0.6.25 stringi_1.5.3 openxlsx_4.2.3
#[33] grid_4.0.3 cli_2.1.0 tools_4.0.3 magrittr_2.0.1 tibble_3.0.4 crayon_1.3.4 tidyr_1.1.3 pkgconfig_2.0.3
#[41] ellipsis_0.3.2 Matrix_1.2-18 data.table_1.13.2 assertthat_0.2.1 rstudioapi_0.11 R6_2.4.1 igraph_1.2.6 compiler_4.0.3
#[49] invctr_0.1.0
```
Sign in with Wallet
Connect another wallet