Skip to content

Instantly share code, notes, and snippets.

@tjmahr
Created February 26, 2019 20:45
Show Gist options
  • Save tjmahr/0d2b41ea1525205a99b19770fc916a90 to your computer and use it in GitHub Desktop.
Save tjmahr/0d2b41ea1525205a99b19770fc916a90 to your computer and use it in GitHub Desktop.
difference of two smooths
library(PupillometryR)
#> Loading required package: dplyr
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union
#> Loading required package: ggplot2
data("pupil_data")

#Check that IDs are not numeric
pupil_data$ID <- as.character(pupil_data$ID)
#remove participant number 8, who had problematic data
pupil_data <- subset(pupil_data, ID != 8)
#blinks were registered as -1, so replace with NAs
pupil_data$LPupil[pupil_data$LPupil == -1] <- NA
pupil_data$RPupil[pupil_data$RPupil == -1] <- NA

Sdata <- make_pupillometryr_data(data = pupil_data,
                                 subject = ID,
                                 trial = Trial,
                                 time = Time,
                                 condition = Type)


new_data <- replace_missing_data(data = Sdata,
                                 rpupil = RPupil,
                                 lpupil = LPupil)
#> Warning in replace_missing_data(data = Sdata, rpupil = RPupil, lpupil
#> = LPupil): replace_missing_data will only help if you have missing
#> timepoints, and a reliable time column.
plot(new_data, pupil = LPupil, group = 'condition')
#> Warning: Removed 3639 rows containing non-finite values (stat_summary).

regressed_data <- regress_data(data = new_data,
                               pupil1 = RPupil,
                               pupil2 = LPupil)
mean_data <- calculate_mean_pupil_size(data = regressed_data, rpupil = RPupil, lpupil = LPupil)
filtered_data <- filter_data(data = mean_data,
                             pupil = mpupil,
                             filter = 'hanning',
                             degree = 11)

int_data <- interpolate_data(data = filtered_data,
                             pupil = mpupil,
                             type = 'linear')

#plot
plot(int_data, pupil = mpupil, group = 'subject')

base_data <- baseline_data(data = int_data,
                           pupil = mpupil,
                           start = 0,
                           stop = 100)

plot(base_data, pupil = mpupil, group = 'subject')

library(mgcv)
#> Loading required package: nlme
#> 
#> Attaching package: 'nlme'
#> The following object is masked from 'package:dplyr':
#> 
#>     collapse
#> This is mgcv 1.8-26. For overview type 'help("mgcv-package")'.


base_data$IDn <- as.numeric(base_data$ID)
base_data$Typen <- ifelse(base_data$Type == 'Easy', .5, -.5)
base_data$Trialn <- as.numeric(substr(base_data$Trial, 5, 5))
base_data$Trialn <- ifelse(base_data$Typen == .5, base_data$Trialn, base_data$Trialn + 3)
base_data$ID <- as.factor(base_data$ID)

model_data <- base_data

m1 <- bam(mpupil ~ s(Time) + s(Time,  by = Typen),
          data = base_data,
          family = gaussian)

library(itsadug)
#> Loading required package: plotfunctions
#> 
#> Attaching package: 'plotfunctions'
#> The following object is masked from 'package:ggplot2':
#> 
#>     alpha
#> Loaded package itsadug 2.3 (see 'help("itsadug")' ).

diffs <- get_difference(
  m1,
  comp = list(Typen = c(-0.5, 0.5)),
  cond = list(Time = unique(base_data$Time)),
  print.summary = TRUE) 
#> Summary:
#>  * Time : numeric predictor; with 600 values ranging from 16.666667 to 10000.000000.

ggplot(diffs) + 
  aes(x = Time, y = difference) +
  geom_hline(yintercept = 0, color = "white", size = 2) + 
  geom_ribbon(
    aes(ymin = difference - CI, ymax = difference + CI), 
    alpha = .3) +
  geom_line() +
  ggtitle("This confidence interval is probably too narrow.") + 
  labs(subtitle = "The model should probably have by-subject random smooths.")

Created on 2019-02-26 by the reprex package (v0.2.1)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment