Skip to contents

Fit Julia regression models to each time point of a time series data

Usage

compute_timewise_statistics(
  jlmer_spec,
  family = c("gaussian", "binomial"),
  statistic = c("t", "chisq"),
  ...
)

Arguments

jlmer_spec

Data prepped for jlmer from make_jlmer_spec()

family

A GLM family. Currently supports "gaussian" and "binomial".

statistic

Test statistic for calculating cluster mass. Can be one of "t" (default) from the regression model output or "chisq" from a likelihood ratio test (takes about twice as long to calculate).

...

Optional arguments passed to Julia for model fitting. Defaults to fast = TRUE (when family = "binomial") and progress = FALSE.

Value

A predictor-by-time matrix of cluster statistics, of class timewise_statistics.

Examples

# \donttest{
# \dontshow{
options("jlmerclusterperm.nthreads" = 2)
jlmerclusterperm_setup(cache_dir = tempdir(), verbose = FALSE)
julia_progress(show = FALSE)
# }

library(dplyr, warn.conflicts = FALSE)

# Specification object
spec <- make_jlmer_spec(
  weight ~ 1 + Diet, filter(ChickWeight, Time <= 20),
  subject = "Chick", time = "Time"
)
spec
#> ── jlmer specification ───────────────────────────────────────── <jlmer_spec> ──
#> Formula: weight ~ 1 + Diet2 + Diet3 + Diet4
#> Predictors:
#>   Diet: Diet2, Diet3, Diet4
#> Groupings:
#>   Subject: Chick
#>   Trial:
#>   Time: Time
#> Data:
#> # A tibble: 533 × 6
#>   weight Diet2 Diet3 Diet4 Chick  Time
#>    <dbl> <dbl> <dbl> <dbl> <ord> <dbl>
#> 1     42     0     0     0 1         0
#> 2     51     0     0     0 1         2
#> 3     59     0     0     0 1         4
#> # ℹ 530 more rows
#> ────────────────────────────────────────────────────────────────────────────────

# Predictor x Time matrix of t-statistics from regression output
empirical_statistics <- compute_timewise_statistics(spec)
round(empirical_statistics, 2)
#>          Time
#> Predictor     0    2    4    6    8   10   12   14   16   18   20
#>     Diet2 -1.60 1.67 2.60 3.52 2.38 1.91 1.94 1.34 1.16 1.40 1.51
#>     Diet3 -1.37 2.45 4.48 4.54 3.70 2.97 3.06 2.98 3.05 3.60 3.80
#>     Diet4 -0.92 3.53 6.27 7.00 5.12 4.07 3.66 2.78 2.16 2.13 2.64

# Collect as dataframe with `tidy()`
empirical_statistics_df <- tidy(empirical_statistics)
empirical_statistics_df
#> # A tibble: 33 × 3
#>    predictor  time statistic
#>    <chr>     <dbl>     <dbl>
#>  1 Diet2         0    -1.60 
#>  2 Diet3         0    -1.37 
#>  3 Diet4         0    -0.916
#>  4 Diet2         2     1.67 
#>  5 Diet3         2     2.45 
#>  6 Diet4         2     3.53 
#>  7 Diet2         4     2.60 
#>  8 Diet3         4     4.48 
#>  9 Diet4         4     6.27 
#> 10 Diet2         6     3.52 
#> # ℹ 23 more rows

# Timewise statistics are from regression models fitted to each time point
# - Note the identical statistics at `Time == 0`
empirical_statistics_df %>%
  filter(time == 0)
#> # A tibble: 3 × 3
#>   predictor  time statistic
#>   <chr>     <dbl>     <dbl>
#> 1 Diet2         0    -1.60 
#> 2 Diet3         0    -1.37 
#> 3 Diet4         0    -0.916
to_jlmer(weight ~ 1 + Diet, filter(ChickWeight, Time == 0)) %>%
  tidy() %>%
  select(term, statistic)
#> # A tibble: 4 × 2
#>   term        statistic
#>   <chr>           <dbl>
#> 1 (Intercept)   164.   
#> 2 Diet2          -1.60 
#> 3 Diet3          -1.37 
#> 4 Diet4          -0.916

# \dontshow{
JuliaConnectoR::stopJulia()
# }
# }