Recognising human activity

# Loading packages
library(tidyverse) 
library(dplyr)
library(readr)
library(caret)
library(nnet)

# Creating a path to the files
list.files(path = "../input") 

# Coping all files to the working directory
system(paste0("cp -r ", list.files("../input", pattern = "recognition", full.names=TRUE), "/* ./"))

# List files in the current working directory
list.files()

# Read the text file containing the activity labels
act_labels = read_delim("activity_labels.txt"," ",col_names=FALSE,trim_ws=TRUE) 
act_labels = act_labels %>% select(X1,X2)
act_labels

# Loading the signals from the training data and adding column names
labels = read_delim("./RawData/Train/labels_train.txt", " ", col_names = F)
colnames(labels) <- c('trial', 'userid', 'activity', 'start', 'end')

# Matching the labels with the activites
labels = labels %>% mutate(activity = act_labels$X2[activity])

# Look at our new dataframe
print(labels)

# Add the sequence start:end to each row in a list, creating a nester tibble
sample_labels_nested = 
    labels %>% 
    rowwise() %>% # do next operation(s) rowwise
    mutate(sampleid = list(start:end)) %>%
    ungroup()

# Check the resulting table
print(sample_labels_nested)

# Unnest the nested tabel.
sample_labels = 
    sample_labels_nested %>% 

    # Rows are segments, we need to keep track of different segements
    mutate(segment = row_number() ) %>% 

    # Expand the data frame to one sample per row
    unnest(cols = c(sampleid)) %>% 

    # Remove columns we don't need anymore
    select(-start, -end) 

# Check the result (first few rows are not interesting; rows 977-990 are)
print(sample_labels[0:1000, ])

# Identify the file name and extract the 'username' (participant ID) and 'expname' (experimental trial)
filename_acc = "RawData/Train/acc_exp01_user01.txt"
username = gsub(".+user(\\d+).+", "\\1", filename_acc) %>% as.integer()
expname  = gsub(".+exp(\\d+).+", "\\1", filename_acc) %>% as.integer()

# Import the data from the file
acc_data = read.delim(filename_acc, " ", header = FALSE) %>%
  setNames(c("X1", "X2", "X3"))

# Creating a plot
options(repr.plot.width=20)
plot.ts(acc_data, xlab="Sample number")

# Identify the file name and extract the 'username' (participant ID) and 'expname' (experimental trial)
filename_gyro = "RawData/Train/gyro_exp01_user01.txt"
username = gsub(".+user(\\d+).+", "\\1", filename_gyro) %>% as.integer()
expname  = gsub(".+exp(\\d+).+", "\\1", filename_gyro) %>% as.integer()

# Read gyro_data and give names to variables
gyro_data <- read.delim(filename_gyro, " ", header = FALSE) %>%
  setNames(c("X1", "X2", "X3"))

# Creating a plot
options(repr.plot.width=20)
plot.ts(gyro_data, xlab="Sample number")

# Combine accelometer and gyroscope data
combined_data <- cbind(acc_data, gyro_data)

# Store signals of the first participant for both acc and gyro data in a data frame, with 'userid' and 'trial'
user_df = 
    data.frame(userid = username, trial = expname, combined_data) %>%
    mutate(sampleid = 0:(nrow(combined_data)-1) ) %>%
    left_join(sample_labels) 

# Check the result (first few rows are not interesting; the following are)
print(user_df[1227:1239, ])

# Set plot width options for the visualization
options(repr.plot.width=20) 

# Create a ggplot visualization using user_df data
user_df %>%
  ggplot(aes(x = sampleid, y = X1, col = factor(activity), group = segment)) + 

# Add a line plot to the ggplot visualization
  geom_line()

# Set plot width options for the visualization
options(repr.plot.width=20) 

# Create a ggplot visualization using user_df data
user_df %>%
  ggplot(aes(x = sampleid, y = X1.1, col = factor(activity), group = segment)) + 

# Add a line plot to the ggplot visualization
  geom_line()

# Split data into epochs and add epoch ID
user_df %>%
  mutate(epoch = sampleid %/% 128) %>% 
  group_by(epoch) %>%
  head(5)

# Lagged cor feature
lagged_cor = function(x, y=x, lag=0) {
   lagged_cor = cor(x, dplyr::lag(y, lag), use='pairwise')
    return(lagged_cor)
}

# Skewness feature
skewness <- function(x) {
  n <- length(x)
  mean_x <- mean(x, na.rm = TRUE)
  s <- sqrt(sum((x - mean_x)^2, na.rm = TRUE) / (n - 1))
  skewness <- sum(((x - mean_x) / s)^3, na.rm = TRUE) / n
  return(skewness)
}

# Kurtosis feature
kurtosis <- function(x) {
  n <- length(x)
  mean_x <- mean(x, na.rm = TRUE)
  s <- sqrt(sum((x - mean_x)^2, na.rm = TRUE) / (n - 1))
  kurtosis <- sum(((x - mean_x) / s)^4, na.rm = TRUE) / n - 3
  return(kurtosis)
}

# Energy feature
energy = function(x) {
  mean(x^2)
}

# Spectral mean
spe_mean <- function(var) {
    spec <- spectrum(var, plot = FALSE)
    df <- spec$freq[2] - spec$freq[1]
    spe_mean <- sum(spec$freq * spec$spec * df)
    return(spe_mean)
}

# Spectral variance
spe_var <- function(var) {
    spec <- spectrum(var, plot = FALSE)
    df <- spec$freq[2] - spec$freq[1]
    spe_var <- sum((spec$freq - mean(var))^2 * spec$spec * df)
    return(spe_var)
}

# Spectral peak frequency function
spectral_peak = function(x) {
  spec = spectrum(x, plot = FALSE) 
  peak = spec$freq[which.max(spec$spec)]
  return(peak)                      
}

# Most common value feature
most_common_value = function(x) {
    counts = table(x, useNA='no')
    most_frequent = which.max(counts)
    return(names(most_frequent))
}

# Magnitude feature
rms <- function(x){
  sqrt((1/length(x)*(sum(x^2))))
}

# Listing all raw data files
list.files("./RawData/Train/", pattern = "^acc")
list.files("./RawData/Train/", pattern = "^gyro")

# Feature function
extractFeatures <- function(filename, sample_labels) {
    
    # Extract user and experimental run ID's from file name
    username = gsub(".+user(\\d+).+", "\\1", filename) %>% as.numeric()
    expname  = gsub( ".+exp(\\d+).+", "\\1", filename) %>% as.numeric()
    
    # Import the sensor signals from the file
    user01 <- read_delim(filename, " ", col_names = F, progress = TRUE, 
                 col_types = "ddd")
    
    # Merge signals with labels 
    user_df <- 
        data.frame(userid = username, trial = expname, user01) %>%
        mutate(sampleid = 0:(nrow(user01)-1) ) %>%
        left_join(sample_labels, by = c('userid','trial','sampleid')) 
    
    # Split in epochs of 128 samples and compute features per epoch
    usertimedom <-  user_df %>%
    
          # Add an epoch ID variable (on epoch = 2.56 sec)
          mutate(epoch = sampleid %/% 128) %>% 

          # Extract statistical features from each epoch
          group_by(epoch) %>%
          summarise(
            # Keep track of user and experiment information
            user_id = username, 
            exp_id = expname,   
              
            # Epoch's activity labels and start sample
            activity = most_common_value(c("-", activity)),
            sampleid = sampleid[1],
              
            # Adding Features  ---------------------------
            # Mean
            m1 = mean(X1), 
            m2 = mean(X2),
            m3 = mean(X3), 
              
            # Mean difference
            mdif1_2 = m1 - m2, 
            mdif1_3 = m1 - m3,
            mdif2_3 = m2 - m3,
            
            # Minimum
            min1 = min(X1),
            min2 = min(X2),
            min3 = min(X3),
              
            # Maximum
            max1 = max(X1),
            max2 = max(X2),
            max3 = max(X3),
              
            # Standard deviation
            sd1 = sd(X1), 
            sd2 = sd(X2), 
            sd3 = sd(X3), 
      
            # Median
            median1 = median(X1),
            median2 = median(X2),
            median3 = median(X3),
              
            # Range
            range1 = max(X1) - min (X1),
            range2 = max(X2) - min (X2),
            range3 = max(X3) - min (X3),
            
            # 25th percentile
            q1_25 = quantile(X1, .25),
            q2_25 = quantile(X2, .25),
            q3_25 = quantile(X3, .25),
            
            # 50th percentile
            q1_50 = quantile(X1, .50),
            q2_50 = quantile(X2, .50),
            q3_50 = quantile(X3, .50),
              
            # 75th percentile
            q1_75 = quantile(X1, .75),
            q2_75 = quantile(X2, .75),
            q3_75 = quantile(X3, .75),
            
            # Autocorrelation for X1(t1) and X1(t2) etc.
            AR1.1 = cor(X1, lag(X1), use = "pairwise"),
            AR1.2 = cor(X1, lag(X1, n = 2), use = "pairwise"),
            AR2.1 = cor(X1, lag(X2), use = "pairwise"),
            AR2.2 = cor(X1, lag(X2, n = 2), use = "pairwise"),
            AR3.1 = cor(X1, lag(X3), use = "pairwise"),
            AR3.2 = cor(X1, lag(X3, n = 2), use = "pairwise"),
        
            # Cross-lagged correlation for X1(t1) and X2(t0) etc.
            AR12.1 = cor(X1, lag(X2), use = "pairwise"),
            AR12.2 = cor(X1, lag(X2), use = "pairwise"),
            AR13.1 = cor(X1, lag(X2), use = "pairwise"),
            AR13.2 = cor(X1, lag(X2), use = "pairwise"),
            AR23.1 = cor(X1, lag(X2), use = "pairwise"),
            AR23.2 = cor(X1, lag(X2), use = "pairwise"),
              
            # Skewness
            skew1 = e1071::skewness(X1),
            skew2 = e1071::skewness(X2),
            skew3 = e1071::skewness(X3),
              
            # Kurtosis
            kurt1 = e1071::kurtosis(X1),
            kurt2 = e1071::kurtosis(X2),
            kurt3 = e1071::kurtosis(X3),
      
            # Energy
            energy_1 = mean(X1^2),
            energy_2 = mean(X2^2),
            energy_3 = mean(X3^2),
              
            # Spectral mean
            spec_mean_1 = spe_mean(X1),
            spec_mean_2 = spe_mean(X2),
            spec_mean_3 = spe_mean(X3),
              
            # Spectral variance
            spec_var_1 = spe_var(X1),
            spec_var_2 = spe_var(X2),
            spec_var_3 = spe_var(X3),
                            
            # Spectral peak
            spectral_peak1 = spectral_peak(X1),
            spectral_peak2 = spectral_peak(X2),
            spectral_peak3 = spectral_peak(X3),
                            
            # RMS
            rms_1 = rms(X1),
            rms_1 = rms(X2),
            rms_1 = rms(X3),
              
            n_samples = n()  
          ) 
    
    usertimedom 
}


# Loading both acc and gyro data ---------------------------
filename_acc <- dir("./RawData/Train/", "^acc", full.names = TRUE)
filename_gyro <- dir("./RawData/Train/", "^gyro", full.names = TRUE)

# Using map_dfr to run the function `extractFeatures` on all elements
myData_acc = map_dfr(filename_acc, extractFeatures, sample_labels)
myData_gyro = map_dfr(filename_gyro, extractFeatures, sample_labels)

# Check the results
print(myData_acc)
print(myData_gyro)

# Create a combined dataset
myData <-  myData_acc %>%
    left_join(myData_gyro, by = c("user_id", "exp_id", "activity", "sampleid", "epoch", "n_samples"))

# Print the combined dataset
head(myData)

# #Check dim
dim(myData)

# Remove rows with NA values
myData1 <- myData[myData$activity != "-", ]

# Remove rows with epoch not 128
myData2 <- myData1 %>% filter(n_samples == 128)

# Check dim
dim(myData2)

# Coloums we want to keep
drop.cols <- c("epoch", "user_id", "exp_id", "activity", "sampleid", "n_samples")
check_data <- myData2 %>% 
            select(-one_of(drop.cols))

# Check whether any predictors are at or near zero variance
nzv <- nearZeroVar(check_data)
nzv_names <- colnames(check_data)[nzv]
nzv_names

# Remove any predictors near zero variance
myData3 <- myData2 %>% 
    select(-one_of(nzv_names))

# Check dim
dim(myData3)

# Create a dataframe containing only numeric values
myData3_num <- myData3 %>% select_if(is.numeric)

# Remove constant variables
myData3_num <- myData3_num %>% 
  select(-where(~ sd(.) == 0))

# Check for missing values
if (anyNA(myData3_num)) {
  myData3_num <- myData3_num %>% na.omit()  # This removes rows with missing values
}

# Calculate the correlation matrix
R <- cor(myData3_num, use = "pairwise.complete.obs")

# Specify the variables to keep
variables_to_keep <- c("epoch", "user_id", "exp_id", "activity", "sampleid", "n_samples")

# Identify the indices of variables to remove based on high correlation
cor_indices <- caret::findCorrelation(R, cutoff = 0.8)  

# Find the variable names to remove, excluding the specified variables to keep
vars_to_remove <- names(myData3_num)[cor_indices]
vars_to_remove <- setdiff(vars_to_remove, variables_to_keep)

# Remove the highly correlated variables from myData2
myData4 <- myData3 %>% select(-one_of(vars_to_remove))

# Check the resulting dataset
head(myData4)

# Check dim
dim(myData4)

# Creating a new dataframe only including the classes to predict and the predictors
cleaned_data <-  myData4 %>%
    select(-epoch, -user_id, -exp_id, -sampleid, -n_samples)

# Visualizing the clean data
head(cleaned_data)

# Check dim
dim(cleaned_data)

# Cross validate our models using 10 folds
trcntr <- trainControl('cv', number = 10, p = 0.8)

# Fitting Multinomial logistic regression
fit_multinom = caret::train(activity ~ ., data = cleaned_data, 
                       method = "multinom", 
                       trControl = trcntr, 
                       MaxNWts = 1500, 
                       verbose = FALSE, 
                       trace = FALSE)
# Storing Results
fit_multinom$result

# Fitting LDA
fit_lda <- caret::train(activity ~., data = cleaned_data, 
                        method = 'lda', 
                        trControl = trcntr)
# Storing Results
fit_lda$result

# FittingkNN 
fit_knn = caret::train(activity ~ ., data = cleaned_data, 
                       method="knn", 
                       trControl = trcntr)
# Storing Results
fit_knn$result

# Fitting KNN
fit_KNN1 <- caret::train(activity ~ ., data = cleaned_data, 
                       method = "knn", 
                       trControl = trcntr,
                       preProcess = "scale",
                       metric = "Accuracy")
# Storing Results
fit_KNN1$results

# Create a data frame with model names and accuracy scores
model_names <- c("multinom", "lda", "knn", "knn1")  # Update with your actual model names
Acc <- sapply(list(fit_multinom, fit_lda, fit_knn, fit_KNN1), function(mdl) max(mdl$results$Accuracy))
model_data <- data.frame(Model = model_names, Accuracy = Acc)

# Find the best performing model
best_model <- model_data[which.max(model_data$Accuracy), ]

# Create a bar plot with ggplot2 with pastel blue and green colors, increased font size, and thinner bars
ggplot(model_data, aes(x = Model, y = Accuracy, fill = Model)) +
  geom_bar(stat = "identity", width = 0.5) +  # Adjust the width parameter for thinner bars
  scale_fill_manual(values = c("multinom" = "#66c2a5", "lda" = "#fc8d62", "knn" = "#8da0cb", "knn1" = "#e78ac3")) +  # Specify pastel colors
  labs(title = "Model Performance for Physical Activity Classifier",
       x = "Model",
       y = "Accuracy") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 30, hjust = 1, size = 30), 
        axis.text.y = element_text(size = 30),  
        axis.title = element_text(size = 30),   
        plot.title = element_text(size = 30))

# Iterate through accelerometer filenames
filenames_test_acc = list.files("./RawData/Test/", "^acc", full.names = TRUE)
test_acc <- map_dfr(filenames_test_acc, extractFeatures, sample_labels) 

# Iterate through gyroscope filenames
filenames_test_gyro = list.files("./RawData/Test/", "^gyro", full.names = TRUE)
test_gyro <- map_dfr(filenames_test_gyro, extractFeatures, sample_labels)

# Merge together
myData_test <-  test_acc %>%
    left_join(test_gyro, 
              by = c("user_id", "exp_id", "activity", "sampleid", "epoch", "n_samples"))

# Take a look
head(myData_test)

# Create a df to use for creating predictions
pred_df <- myData_test %>%
  select(-epoch, -user_id, -exp_id, -sampleid, n_samples)

# We create predictinos based on fit_lda_clean
predictions <- predict(object = fit_multinom, newdata = myData_test[,-1])

# Put the predictions into test_df
myData_test["activity"] <- predictions

# Take a look
head(myData_test)

# Formating our data for the competition
myData_test %>%

    # Prepend "user" and "exp" to user_id and exp_id
    mutate(
        user_id = paste(ifelse(user_id < 10, "user0", "user"), user_id, sep=""), 
        exp_id = paste(ifelse(exp_id < 10, "exp0", "exp"), exp_id, sep="")
    ) %>% 

     # Unite columns user_id, exp_id and sample_id into a string 
     # Separated by "_" and store it in the new variable `Id`
    unite(Id, user_id, exp_id, sampleid) %>%

    # Retain only the `Id` and  predictions
    select(Id, Predicted = activity) %>%

    # Write to file
    write_csv("submission_final.csv")

# Check the result: print first 20 lines in the submission file
cat(readLines("submission_final.csv",20), sep="\n")