1
$\begingroup$

I have generated a time series data set of measurements that are a bit noisy and I want to apply kernel smoothing to the data. My time series data is not regular however, meaning that the time difference between consecutive observations can be variable so I want to account for that within the weighting of the kernel smoothing.

I give an example dataset below and apply kernel smoothing with a couple different bandwidth sizes as well as a rolling mean with different window sizes for comparison.

My impression was that kernel smoothing worked by taking a weighted average of the adjacent points weighted by how close they are (here in time) to the "focal observation". The bandwidth controls the relative contribution of adjacent points (essentially the weights given to adjacent points given how close they are) where larger bandwidths results in more smoothing (i.e. larger weights to closer observations). In contrast, a rolling mean takes a regular uniform average of observations within a window around the focal observation and so that all observations within the window contribute equally to the mean regardless of how far away they are from the focal observation.

However, when I applied this in practice, I don't see this in my results. Is my understanding wrong? Is my code wrong?

What I see is that when I apply kernel smoothing with a bandwidth of 1/2 or 1, it doesn't change the raw data, when I apply kernel smoothing with a bandwidth of 2 or 3, the results are equivalent, and when I apply a rolling mean of window size +/- 1 or +/- 2 these are equivalent to a kernel smooth with bandwidths of 2 or 4 respectively.

Can anyone explain to me what's going on?

Thanks!

#################
# Simulate data #
#################
# Simulate a time series data set with different underlying processes that
# last a random length with variable sampling frequencies
library(dplyr)
library(ggplot2)
library(tidyverse)
# Set random seed
set.seed(1)
# Total # of process bouts
n_mode = 5
# Average sampling rate
sample_rate = 10 #hz
# Durations can be irregularly long
durs = rgamma(n_mode, shape = 2.5, rate = .5)
# Total length of time collected
times = cumsum(durs) # Ending times
# Random order of modes
order = c(3,1,3,2,1)
# Stick in a data frame
df = data.frame(modes = order, 
                duration = durs, 
                end_time = times)
# Determine starting timestamp of bout
df$start_time = c(0, times[1:(length(times)-1)])
# Reorder columns
df = df %>% select(modes, duration, start_time, end_time)
# For each mode
for(i in 1:nrow(df)){
  # By each row
  df_sample = df[i,]
  # Generate an idealized sequence of data with perfect really fast sampling
  ts_ideal = seq(df$start_time[i], df$end_time[i], by = 1/10000)
  # Randomly sample observations so that the total sampling rate ends up
  # at an average of observations a second
  time_ss = sample(ts_ideal, size = sample_rate * df_sample$duration, replace = FALSE)
  # Put in order
  time_ss = time_ss[order(time_ss)]
  # Get "time stamp" second
  time_s = floor(time_ss)
  # Extract sub-second info
  ss = time_ss - time_s
  # Put into time df
  times_df = data.frame(index = i, mode = df_sample$mode, time_s, ss, time_ss)
  # If first iteration, generate combined df
  # Otherwise add to combined df
  if(i == 1){comb = times_df}else(comb = rbind(comb, times_df))
}
# Generate variable by process type (mode)
for(i in 1:length(unique(comb$index))){
  # subset by index
  df_index = comb[comb$index == i, ]
  # given the mode
  if(unique(df_index$mode) == 1){
    # signals
    df_index$x_raw = rep(0, length(df_index$time_ss))}
  if(unique(df_index$mode) == 2){
# signals
df_index$x_raw = sin(2 * pi * df_index$time_ss)}
  if(unique(df_index$mode) == 3){
    # signals
    df_index$x_raw = rep(1, length(df_index$time_ss))}
  # Combine
  if(i == 1){df_comb = df_index}else(df_comb = rbind(df_comb, df_index))
}
# Add some error (noise)
df_comb$x_raw = df_comb$x_raw + rnorm(n = nrow(df_comb), mean = 0, sd = .05)
# Define as factor
df_comb$mode = as.factor(df_comb$mode)
# Pull out transition times of x for visualization
trans_times = df_comb$time_ss[df_comb$index!=lag(df_comb$index)]
# Remove NA
trans_times = trans_times[!is.na(trans_times)]

#############
# Visualize #
#############
df_comb %>% 
  ggplot(aes(x = time_ss, y = x_raw, color = mode, group = 1))+
  geom_point(aes(color = mode, group = 1))+
  geom_line(aes(color = mode, group = 1))+
  scale_color_manual(values = c("red", "blue","green"))+ #Color by mode
  geom_vline(xintercept = trans_times, linetype = "dashed") # Delineate transitions

####################
# Kernel smoothing #
####################
df_comb$x_ks1_2 = ksmooth(time(df_comb$time_ss), df_comb$x_raw, bandwidth = 0.5)$y
# Bandwidth is 1 (does not appear to work)- is same as raw
df_comb$x_ks1 = ksmooth(time(df_comb$time_ss), df_comb$x_raw, bandwidth = 1)$y
# Bandwidth is 2 (first that appears to work)
df_comb$x_ks2 = ksmooth(time(df_comb$time_ss), df_comb$x_raw, bandwidth = 2)$y
# Bandwidth is 3 (is same as 2)
df_comb$x_ks3 = ksmooth(time(df_comb$time_ss), df_comb$x_raw, bandwidth = 3)$y
# Bandwidth is 4 (higher smoothing)
df_comb$x_ks4 = ksmooth(time(df_comb$time_ss), df_comb$x_raw, bandwidth = 4)$y

########
# Plot #
########
df_comb %>% pivot_longer(cols = x_raw:x_ks4, names_to = "type", values_to = "x") %>% 
  ggplot(aes(x = time_ss, y = x, group = type, color = type))+
  geom_point()+
  geom_line()+
  scale_color_manual(values = c("x_raw" = "black", 
                                "x_ks1_2" = "red", 
                                "x_ks1" = "orange",
                                "x_ks2" = "green", 
                                "x_ks3" = "blue",
                                "x_ks4" = "purple"))+ #color by bandwidth
  geom_vline(xintercept = trans_times, linetype = "dashed") # Delineate transitions

##############################################################
# Why do different bandwidths produce the same numbers ??!! #
#############################################################
# Kernel smoothing bandwidth of 3 same as kernel smoothing bandwidth of 2?
sum(df_comb$x_ks2 == df_comb$x_ks3)/nrow(df_comb) # 100% same
# Kernel smoothing bandwidth of 1 and kernel smoothing bandwidth of 1/2 same as raw?
sum(df_comb$x_ks1 == df_comb$x_raw)/nrow(df_comb) # 100% same
sum(df_comb$x_ks1_2 == df_comb$x_raw)/nrow(df_comb) # 100% same

##################
# Rolling Window #
##################
# Quarter second window (+/- 1 observation)
setDT(df_comb)[, "x_rm1_4" := frollapply(x_raw,
                       n = (floor(sample_rate/4) + 1),
                       FUN = mean, 
                       align = "center",
                       adaptive = F,
                       na.rm = T)]
# Half second window (+/- 2 observations)
df_comb[, "x_rm1_2" := frollapply(x_raw,
                       n = sample_rate/2,
                       FUN = mean, 
                       align = "center", 
                       na.rm = T)]
########
# Plot #
########
df_comb %>% pivot_longer(cols = c(x_raw, 
                                  x_rm1_4,
                                  x_rm1_2), names_to = "type", values_to = "x") %>% 
  ggplot(aes(x = time_ss, y = x, group = type, color = type))+
  geom_point()+
  geom_line()+
  scale_color_manual(values = c("x_raw" = "black", 
                                "x_rm1_4" = "green", 
                                "x_rm1_2" = "purple"))+ #color by smoothing window width
  geom_vline(xintercept = trans_times, linetype = "dashed") # Delineate transitions

###################################################
# Plot both Kernel and Rolling Mean on same Graph #
###################################################
df_comb %>% pivot_longer(cols = c(x_raw, 
                                  x_rm1_4,
                                  x_ks2,
                                  x_rm1_2,
                                  x_ks4
                                  ), names_to = "type", values_to = "x") %>% 
  ggplot(aes(x = time_ss, y = x, group = type, color = type))+
  geom_point()+
  geom_line()+
  scale_color_manual(values = c("x_raw" = "black", 
                                "x_rm1_4" = "red",
                                "x_ks2" = "orange",
                                "x_rm1_2" = "green",
                                "x_ks4" = "blue"
                                ))+ #color by smoothing window width
  geom_vline(xintercept = trans_times, linetype = "dashed") # Delineate transitions

# Plot by each mode in case differences are very tiny
df_comb %>% filter(index == 1) %>% pivot_longer(cols = c(x_raw, 
                                  x_rm1_4,
                                  x_ks2,
                                  x_rm1_2,
                                  x_ks4
), names_to = "type", values_to = "x") %>% 
  ggplot(aes(x = time_ss, y = x, group = type, color = type))+
  geom_point()+
  geom_line()+
  scale_color_manual(values = c("x_raw" = "black", 
                                "x_rm1_4" = "red",
                                "x_ks2" = "orange",
                                "x_rm1_2" = "green",
                                "x_ks4" = "blue"
  ))
df_comb %>% filter(index == 2) %>% pivot_longer(cols = c(x_raw, 
                                                         x_rm1_4,
                                                         x_ks2,
                                                         x_rm1_2,
                                                         x_ks4
), names_to = "type", values_to = "x") %>% 
  ggplot(aes(x = time_ss, y = x, group = type, color = type))+
  geom_point()+
  geom_line()+
  scale_color_manual(values = c("x_raw" = "black", 
                                "x_rm1_4" = "red",
                                "x_ks2" = "orange",
                                "x_rm1_2" = "green",
                                "x_ks4" = "blue"
  ))
##########################################################
# Why is kernel smoothing the same as rolling means ??!! #
##########################################################
# Round to same number of decimal places so precision isn't the underlying difference
df_comb$x_rm1_4 = round(df_comb$x_rm1_4, 12)
df_comb$x_rm1_2 = round(df_comb$x_rm1_2, 12)
df_comb$x_ks2 = round(df_comb$x_ks2, 12)
df_comb$x_ks4 = round(df_comb$x_ks4, 12)
# Pull out observations where the rolling mean and kernel smooth is not the same
df_comb[df_comb$x_rm1_4!=df_comb$x_ks2,]
df_comb[df_comb$x_rm1_2!=df_comb$x_ks4,]
# Eyeball the data
cbind(df_comb[1:6, "x_rm1_4"], df_comb[1:6, "x_ks2"], df_comb[1:6, "x_rm1_4"] == df_comb[1:6, "x_ks2"])
```
Improve this question
$\endgroup$
1
  • $\begingroup$ It's very good that the OP solved their own question, but I find it hard to imagine that this is both on-topic and likely to be useful to anyone else. Too much depends on a reader wading through lengthy code -- or even on running it -- to find out exactly what was done. (R is, I guess, the most commonly used software here, but it's not universal, nor is any other software.) $\endgroup$
    – Nick Cox
    Commented Feb 16 at 10:43

1 Answer 1

Reset to default
0
$\begingroup$

So, it turns out that all I needed to do was to remove the time() function within the kernel smoothing function and use my variable by itself (it's already a numeric representation of time, not a posix_ct type object) and it performs as expected. So my code should been:

# THIS WORKS
df_comb$x_ks1_2 = ksmooth(df_comb$time_ss, df_comb$x_raw, bandwidth = 0.5)$y
# THIS DOES NOT
df_comb$x_ks4 = ksmooth(time(df_comb$time_ss), df_comb$x_raw, bandwidth = 4)$y

This is because the time() function converted by numeric time into a sequential time series treating each individual observation as an integer. I think this is why it wouldn't allow me to use an uneven number for my bandwidth or go below 1 for the bandwidth value.

Well, if anyone needs a practice data set (I created it to represent accelerometer data with different behavior types "modes") in which to test out methods and better understand how they work...here you go!

Improve this answer
$\endgroup$
1
  • $\begingroup$ Assuming this is the right answer -the implication to me is that the thread is about using the right code. $\endgroup$
    – Nick Cox
    Commented Feb 16 at 10:44

Not the answer you're looking for? Browse other questions tagged or ask your own question.