In robschick/smokevis: Makes Simulated Smoking Data for Analysis and Visualisation

knitr::opts_chunk$set(warning=FALSE, message=FALSE, fig.height = 7, fig.width = 7)
library(magrittr)

Goal of this is to simply explore making some simulated data that we can visualise with a flexidashboard.

Data Info

The idea is to take some various geocoded locations, and add them in to a data frame. These will be master ones to be recycled. These will contain a few candidate locations like: home, work, grocery, pub1, pub2, etc. This will just be a starting point.

library(smokevis)
locs1

Those are the raw data, and we can visualise them as follows

library(leaflet)
m <- leaflet(data = locs1) %>% setView(lng = -2.9, lat = 56.403, zoom = 11)
m %>% addTiles() %>% addMarkers(lng = locs1$lng, lat = locs1$lat, popup = locs1$location)

Next we might imagine that we want to separate out craving events with smoking events. That might look like this:

pal <- colorFactor(c("navy", "red"), domain = c("crave", "smoke"))

locs1$type <- factor(
    ifelse(runif(nrow(locs1)) > 0.75, "crave", "smoke"),
    c("smoke", "crave"))

m %>% 
  addTiles() %>%
  addCircleMarkers(
  lng = ~lng, lat = ~lat, popup = ~location,
    radius = ~ifelse(locs1$type == "smoke", 6, 10),
    color = ~pal(locs1$type),
    stroke = FALSE, fillOpacity = 0.5
  )

Now just experimenting with some different functionality to toggle smoking events on and off:

m %>% 
  addTiles() %>%
  addCircleMarkers(
    data = subset(locs1, type == 'smoke'),
    lng = ~lng, lat = ~lat, popup = ~location,
    radius = 10,
    color = 'red',
    stroke = FALSE, fillOpacity = 0.5,
    group = 'Smoking'
  ) %>%
  addCircleMarkers(
    data = subset(locs1, type != 'smoke'),
    lng = ~lng, lat = ~lat, popup = ~location,
    radius = 6,
    color = 'blue',
    stroke = FALSE, fillOpacity = 0.5,
    group = 'Craving'
  ) %>%
  # Layers control
  addLayersControl(
    overlayGroups = c("Smoking", "Craving"),
    options = layersControlOptions(collapsed = FALSE)
  )

Mulitple Smoking Events

Ok, let's start building up a data frame of smoking events. I'm going to say this person is an average smoker smoking ~14 cigarettes per day (from this 2014 Scottish Health Survey). I'll simulate a standard week, comprised of 5 workdays, and 2 weekend days.

Here are some simulated data for smoking events. I'll start with a daily tally:

library(lubridate)
nweeks <- 1
nweekday <- nweeks * 5
nwkndday <- nweeks * 2

cigsweek    <- floor(rnorm(nweekday, 16, 3))
cigsweekend <- floor(rnorm(nwkndday, 12, 5))
dfNum <- data.frame(ncigs = c(cigsweekend[1], cigsweek, cigsweekend[2]),
                 dayofweek = wday(1:7, label = TRUE),
                 dayInt = wday(1:7),
                 date = seq.Date(as.Date('2016-06-19'), as.Date('2016-06-25'), by = 'day'))

Let's see the raw data then that summarise the smoking events per day of the week:

library(formattable)
as.htmlwidget(formattable(dfNum, list(
  # ncigs = color_tile("white", "orange"),
  ncigs = color_tile("#FEE5D9", "#EF3B2C"),
  dayofweek = formatter("span",
                    style = x ~ ifelse(x != "Sun" & x != "Sat", style(color = "green", font.weight = "bold"), NA)),
  date = formatter("span",
                    style = x ~ ifelse(dfNum$dayofweek != "Sun" & dfNum$dayofweek != "Sat", style(color = "green", font.weight = "bold"), NA))
)))

With that made, I can now simulate data at the daily scale. This builds off the number of cigarettes in dfNum:

smokeDat <- vector(mode = 'list', (nrow(dfNum)))

for(i in seq_along(lubridate::wday(1:(nweekday + nwkndday)))){

  ncigs <- dfNum$ncigs[i]
  dayVal <- dfNum$dayInt[i]
  dayLab <- dfNum$dayofweek[i]
  dayDate <- dfNum$date[i]

  # On weekend days subject has lower stress, can't smoke at work
  if (dayVal == 1 | dayVal == 7) { 
      stress <- round(tnorm(ncigs, lo = 0, hi = 10, mu = 5, sig = 2), 1)  
      idx <- sample((1:nrow(locs1))[locs1$location != 'work'], ncigs, replace = TRUE) 
    } else {
      stress <- round(tnorm(ncigs, lo = 0, hi = 10, mu = 8, sig = 3), 1)
      idx <- sample(nrow(locs1), ncigs, replace = TRUE)
    }

  dfDay <- data.frame(
    eventID = 1:ncigs,
    craving = round(tnorm(ncigs, lo = 0, hi = 10, mu = 7, sig = 3), 1),
    satisfaction = round(tnorm(ncigs, lo = 0, hi = 10, mu = 6, sig = 3), 1),
    taste = round(tnorm(ncigs, lo = 0, hi = 10, mu = 5, sig = 3), 1),
    reduction = round(tnorm(ncigs, lo = 0, hi = 10, mu = 3, sig = 1), 1),
    stress = stress,
    time = randTimes(ncigs, 
                     st = paste(dayDate, ' 06:30:00', sep = ''), 
                     et = paste(dayDate, ' 22:30:00', sep = '')),
    dayofweek = dayLab,
    lat = jitter(locs1$lat[idx]),
    lng = jitter(locs1$lng[idx]),
    qualLoc = locs1$location[idx],
    smoked = TRUE
  )  
  smokeDat[[i]] <- dfDay
}
sdf <- as.data.frame(data.table::rbindlist(smokeDat))

Here's what the data look like:

# head(sdf, 2)
df <- subset(sdf, select = -c(lat, lng, smoked))
as.htmlwidget(formattable(df, list(
  eventID = color_tile("#FEE5D9", "#EF3B2C", 0.2),
  craving = color_bar("#FCBBA1", 0.2),
  stress = color_bar("#FCBBA1", 0.2),
  satisfaction = color_bar("#A1D99B", 0.2),
  taste = color_bar("#A1D99B", 0.2),
  reduction = formatter("span",
                          style = x ~ style(color = ifelse(rank(-x) <= 5, "green", "gray"),
                                            font.weight = ifelse(rank(-x) <= 5, "bold", 'normal')),
                          x ~ sprintf("%.1f (rank: %2.0f)", x, rank(-x))),
  dayofweek = formatter("span",
                    style = x ~ ifelse(x != "Sun" & x != "Sat", style(color = "green", font.weight = "bold"), NA)),
  time = formatter("span",
                    style = x ~ ifelse(sdf$dayofweek != "Sun" & sdf$dayofweek != "Sat", style(color = "green", font.weight = "bold"), NA))
) ,
  align = c("r","r","r","l","c","r","r","c", 'c'))
)

Let's look at those mulitple smoking events in the map, because we want to be able to try clustering by location.

m %>% 
  addTiles() %>%
  addCircleMarkers(
    data = sdf,
    lng = ~lng, lat = ~lat, popup = ~qualLoc,
    radius = ~craving,
    color = 'red',
    stroke = TRUE, fillOpacity = 0.5,
    group = 'Smoking',
    clusterOptions = markerClusterOptions()
  ) %>%
  addCircleMarkers(
    data = subset(locs1, type != 'smoke'),
    lng = ~lng, lat = ~lat, popup = ~location,
    radius = 6,
    color = 'blue',
    stroke = FALSE, fillOpacity = 0.5,
    group = 'Craving',
    clusterOptions = markerClusterOptions()
  ) %>%
  # Layers control
  addLayersControl(
    overlayGroups = c("Smoking", "Craving"),
    options = layersControlOptions(collapsed = FALSE)
  )

Craving Events

Now let's make some craving data, so we can plot the symbols according to amount of craving. The simulation is essentially as above, but I only simulate stress and craving to mimic what the app currently asks

craveDat <- vector(mode = 'list', (nrow(dfNum)))

for(i in seq_along(lubridate::wday(1:(nweekday + nwkndday)))){

  ncigs <- floor(dfNum$ncigs[i] * 0.75)
  dayVal <- dfNum$dayInt[i]
  dayLab <- dfNum$dayofweek[i]
  dayDate <- dfNum$date[i]

  # On weekend days subject has lower stress, can't smoke at work
  if (dayVal == 1 | dayVal == 7) { 
      stress <- round(tnorm(ncigs, lo = 0, hi = 10, mu = 5, sig = 2), 1)
      idx <- sample((1:nrow(locs1))[locs1$location != 'work'], ncigs, replace = TRUE) 
    } else {
      stress <- round(tnorm(ncigs, lo = 0, hi = 10, mu = 8, sig = 3), 1)
      idx <- sample(nrow(locs1), ncigs, replace = TRUE)
    }

  dfDay <- data.frame(
    eventID = 1:ncigs,
    craving = tnorm(ncigs, lo = 0, hi = 10, mu = 8, sig = 3),
    stress = stress,
    time = randTimes(ncigs, 
                     st = paste(dayDate, ' 06:30:00', sep = ''), 
                     et = paste(dayDate, ' 22:30:00', sep = '')),
    dayofweek = dayLab,
    lat = jitter(locs1$lat[idx]),
    lng = jitter(locs1$lng[idx]),
    qualLoc = locs1$location[idx],
    smoked = FALSE
  )  
  craveDat[[i]] <- dfDay
}
cdf <- as.data.frame(data.table::rbindlist(craveDat))
head(cdf, 2)

Now we can see them on the map along with the smoking data. Here the smoking events are symbol-coded by satisfaction of the cigarette and the craving events are symbol-coded by craving. The events also have popups to return a bit more information:

mypal <- RColorBrewer::brewer.pal(2, 'Dark2')
m %>% 
  addTiles() %>%
  addCircleMarkers(
    data = sdf,
    lng = ~lng, lat = ~lat, popup = paste('Location: ', sdf$qualLoc, '<br>',
                                          'Satisfaction: ', round(sdf$satisfaction, 2)),
    radius = ~satisfaction,
    color = mypal[1],
    stroke = TRUE, fillOpacity = 0.5,
    group = 'Smoking'#, clusterOptions = markerClusterOptions()
  ) %>%
  addCircleMarkers(
    data = cdf,
    lng = ~lng, lat = ~lat, popup = paste('Location: ', cdf$qualLoc, '<br>',
                                          'Craving: ', round(cdf$craving, 2)),
    radius = ~craving,
    color = mypal[2],
    stroke = TRUE, fillOpacity = 0.5,
    group = 'Craving'#, clusterOptions = markerClusterOptions()
  ) %>%
  # Layers control
  addLayersControl(
    overlayGroups = c("Smoking", "Craving"),
    options = layersControlOptions(collapsed = FALSE)
  )

Time Series of Craving

Now I want to add in the day/night time. I need maptools to do this. This will show a static time series of craving.

library(sp)
library(maptools)
home <- matrix(c(locs1$lng[1], locs1$lat[1]), nrow = 1)
Home <- SpatialPoints(home, proj4string=CRS("+proj=longlat +datum=WGS84"))
tseq <- seq(from = as.POSIXct(dfNum$date[1], tz = 'BST'), 
            to = as.POSIXct(dfNum$date[length(dfNum$date)], tz = 'BST'), by = 'days')
sunup <- sunriset(Home, tseq, direction = 'sunrise', POSIXct.out=TRUE)
sundn <- sunriset(Home, tseq, direction = 'sunset', POSIXct.out=TRUE)
xrect <- data.frame(xmin = sundn$time[1:(nrow(sundn)-1)], xmax = sunup$time[2:nrow(sunup)], utc.date = 1)

library(ggplot2)
library(scales)
ggplot(data = sdf, aes(time, craving, group = dayofweek))+
  geom_line(aes(group = dayofweek))+
  annotate('rect', xmin = xrect$xmin, xmax = xrect$xmax, ymin = 0, ymax = Inf, fill = grey(.75), alpha = 0.2)+
  scale_x_datetime(minor_breaks = date_breaks('2 hour'), 
                   limits = c(as.POSIXct(dfNum$date[1], tz = 'BST'), 
                              as.POSIXct(dfNum$date[length(dfNum$date)] + 1, tz = 'BST'))) + 
  theme_bw()

And here's a timeline of stress, with points indicating smoking events (symbol-coded by craving - bigger point = more craving), and highlighting the critical cigarette, which here I've defined as the one with the biggest reduction in craving:

library(dplyr)
keyCigs <- sdf %>% 
  group_by(dayofweek) %>% 
  top_n(n = 1, wt = reduction)

ggplot(data = sdf, aes(time, stress, group = dayofweek))+
  geom_line(aes(group = dayofweek))+
  geom_point(aes(y = 10, size = craving), alpha = 1 / 5)+
  geom_point(data = keyCigs, aes(x = time, y = 10, size = reduction), 
             shape = 21, alpha = 0.9, colour = 'black', fill = mypal[2], stroke = 1)+
  scale_size(range = c(1, 5))+
  annotate('rect', xmin = xrect$xmin, xmax = xrect$xmax, ymin = 0, ymax = Inf, fill = grey(.75), alpha = 0.2)+
  scale_x_datetime(minor_breaks = date_breaks('2 hour'), 
                   limits = c(as.POSIXct(dfNum$date[1], tz = 'BST'), 
                              as.POSIXct(dfNum$date[length(dfNum$date)] + 1, tz = 'BST'))) + 
  theme_bw()

Using dygraphs

So that was a static time series, let's try to make it a bit interactive.

library(dygraphs)
library(xts)
# All "negative" values: craving and stress
stressDat <- as.matrix(sdf[, c(2, 6)])
stressMat <- xts(stressDat, order.by = sdf$time)

# All three "positive" values: satisfcation, taste, reduction
postDat <- as.matrix(sdf[, c(3, 4, 5)])
postMat <- xts(postDat, order.by = sdf$time)

# Make individual time series objects
craving <- xts(sdf$craving, order.by = sdf$time)
satisfaction <- xts(sdf$satisfaction, order.by = sdf$time)
taste <- xts(sdf$taste, order.by = sdf$time)
reduction <- xts(sdf$reduction, order.by = sdf$time)
stress <- xts(sdf$stress, order.by = sdf$time)

dygraph(postMat) %>% 
  dyRangeSelector() %>% 
  dyShading(from = xrect$xmin[1], to = xrect$xmax[1]) %>% 
  dyShading(from = xrect$xmin[2], to = xrect$xmax[2]) %>% 
  dyShading(from = xrect$xmin[3], to = xrect$xmax[3]) %>% 
  dyShading(from = xrect$xmin[4], to = xrect$xmax[4]) %>% 
  dyShading(from = xrect$xmin[5], to = xrect$xmax[5]) %>% 
  dyShading(from = xrect$xmin[6], to = xrect$xmax[6]) 

dygraph(stressMat) %>%
  dyHighlight(highlightSeriesOpts = list(strokeWidth = 3))

dygraph(stressMat) %>%
  dyHighlight(highlightCircleSize = 5, 
              highlightSeriesBackgroundAlpha = 0.2,
              hideOnMouseOut = FALSE)

Grouped Zoom:

dygraph(craving, main = "Craving", group = "negative")
dygraph(stress, main = "Stress", group = "negative")

dygraph(satisfaction, main = "Satisfaction of Cigarette", group = "positive")
dygraph(reduction, main = "Reduction in Craving", group = "positive")
dygraph(taste, main = "Taste of Cigarette", group = "positive")

We can add a tooltip to return a bit more information behind the 'key' cigarette:

ttip <- paste('Location: ', keyCigs$qualLoc[1], 
        '; Satisfaction: ', round(keyCigs$satisfaction[1], 2), 
        '; Craving: ', round(keyCigs$craving[1], 2), 
        '; Reduction: ', round(keyCigs$reduction[1], 2), sep = '')

dygraph(stress, main = "Stress", group = "negative") %>% 
  dyAxis("y", label = "Self-reported Stress Level") %>% 
  dyAnnotation(keyCigs$time[1], text = keyCigs$eventID[1], tooltip = ttip ) %>% 
  dyAnnotation(keyCigs$time[2], text = keyCigs$eventID[2], tooltip = paste('Craving: ', round(keyCigs$craving[2], 2))) %>% 
  dyAnnotation(keyCigs$time[3], text = keyCigs$eventID[3], tooltip = paste('Craving: ', round(keyCigs$craving[3], 2))) %>% 
  dyAnnotation(keyCigs$time[4], text = keyCigs$eventID[4], tooltip = paste('Craving: ', round(keyCigs$craving[4], 2))) %>% 
  dyAnnotation(keyCigs$time[5], text = keyCigs$eventID[5], tooltip = paste('Craving: ', round(keyCigs$craving[5], 2))) %>% 
  dyAnnotation(keyCigs$time[6], text = keyCigs$eventID[6], tooltip = paste('Craving: ', round(keyCigs$craving[6], 2))) %>% 
  dyAnnotation(keyCigs$time[7], text = keyCigs$eventID[7], tooltip = paste('Craving: ', round(keyCigs$craving[7], 2))) %>% 
  dyShading(from = xrect$xmin[1], to = xrect$xmax[1]) %>% 
  dyShading(from = xrect$xmin[2], to = xrect$xmax[2]) %>% 
  dyShading(from = xrect$xmin[3], to = xrect$xmax[3]) %>% 
  dyShading(from = xrect$xmin[4], to = xrect$xmax[4]) %>% 
  dyShading(from = xrect$xmin[5], to = xrect$xmax[5]) %>% 
  dyShading(from = xrect$xmin[6], to = xrect$xmax[6])