circacompare_mixed — circacompare

circacompare_mixed is similar to circacompare but allows for some simple, user-specified random-effects on the rhythmic parameters of choice.

Usage

circacompare_mixed(
  x,
  col_time,
  col_group,
  col_outcome,
  col_id,
  randomeffects = c(),
  period = 24,
  alpha_threshold = 0.05,
  nlme_control = list(),
  nlme_method = "REML",
  weights = NULL,
  suppress_all = FALSE,
  timeout_n = 10000,
  control = list()
)

Arguments

x: data.frame. This is the data.frame which contains the rhythmic data for two groups in a tidy format.
col_time: The name of the column within the data.frame, x, which contains time in hours at which the data were collected.
col_group: The name of the column within the data.frame, x, which contains the grouping variable. This should only have two levels.
col_outcome: The name of the column within the data.frame, x, which contains outcome measure of interest.
col_id: The name of the column within the data.frame, x, which contains the identifying values for the random effect, such as subject_id.
randomeffects: which rhythmic parameters to allow random effects. The default is to include no rhythmic parameters.
period: The period of the rhythm. For circadian rhythms, leave this as the default value, 24.
alpha_threshold: The level of alpha for which the presence of rhythmicity is considered. Default is to 0.05.
nlme_control: A list of control values for the estimation algorithm to replace the default values returned by the function nlme::nlmeControl. Defaults to an empty list.
nlme_method: A character string. If "REML" the model is fit by maximizing the restricted log-likelihood. If "ML" the log-likelihood is maximized. Defaults to "REML".
weights: An optional numeric vector of (fixed) weights internally passed to nlme::nlme() via nlme::varPower(). When present, the objective function is weighted least squares.
suppress_all: Logical. Set to TRUE to avoid seeing errors or messages during model fitting procedure. Default is FALSE. If FALSE, also runs nlme() with verbose = TRUE.
timeout_n: The upper limit for the model fitting attempts. Default is 10000.
control: list. Used to control the parameterization of the model.

Value

list

Examples

# Generate some data with within-id correlation for phase-shift (phi1)
# \donttest{
set.seed(99)
phi1_in <- 3.15

mixed_data <- function(n) {
  counter <- 1
  for (i in 1:n) {
    x <- make_data(k1 = 0, alpha1 = 0, phi1 = rnorm(1, phi1_in, 0.5), hours = 72, noise_sd = 1)
    x$id <- counter
    counter <- counter + 1
    if (i == 1) {
      res <- x
    } else {
      res <- rbind(res, x)
    }
  }
  return(res)
}
df <- mixed_data(20)
out <- circacompare_mixed(
  x = df,
  col_time = "time",
  col_group = "group",
  col_outcome = "measure",
  col_id = "id",
  control = list(grouped_params = c("phi"), random_params = c("phi1"))
)
#> 
#> **Iteration 1
#> LME step: Loglik: -5362.31, nlminb iterations: 1
#> reStruct  parameters:
#>       id 
#> 1.330561 
#> varStruct  parameters:
#> power 
#>     0 
#>  Beginning PNLS step: ..  completed fit_nlme() step.
#> PNLS step: RSS =  2914.277 
#>  fixed effects: -0.003046635  10.00689  0.001165668  3.036078  
#>  iterations: 5 
#> Convergence crit. (must all become <= tolerance = 1e-05):
#>     fixed  reStruct varStruct 
#>  2.841335  0.718343  0.000000 
#> 
#> **Iteration 2
#> LME step: Loglik: -4181.778, nlminb iterations: 3
#> reStruct  parameters:
#>        id 
#> 0.7491786 
#> varStruct  parameters:
#> power 
#>     0 
#>  Beginning PNLS step: ..  completed fit_nlme() step.
#> PNLS step: RSS =  2914.277 
#>  fixed effects: -0.003046635  10.00689  0.001165668  3.036078  
#>  iterations: 1 
#> Convergence crit. (must all become <= tolerance = 1e-05):
#>     fixed  reStruct varStruct 
#>         0         0         0 

# with sample weights (arbitrary weights for demonstration)
sw <- runif(n = nrow(df))
out2 <- circacompare_mixed(
  x = df,
  col_time = "time",
  col_group = "group",
  col_outcome = "measure",
  col_id = "id",
  control = list(grouped_params = c("phi"), random_params = c("phi1")),
  weights = sw
)
#> 
#> **Iteration 1
#> LME step: Loglik: -9253.195, nlminb iterations: 5
#> reStruct  parameters:
#>       id 
#> 3.094586 
#> varStruct  parameters:
#>        power 
#> -0.002104771 
#>  Beginning PNLS step: ..  completed fit_nlme() step.
#> PNLS step: RSS =  2902.204 
#>  fixed effects: -0.003036673  10.00696  0.001187815  2.721892  
#>  iterations: 7 
#> Convergence crit. (must all become <= tolerance = 1e-05):
#>     fixed  reStruct varStruct 
#>   7.58760  17.93843   1.00000 
#> 
#> **Iteration 2
#> LME step: Loglik: -4199.912, nlminb iterations: 7
#> reStruct  parameters:
#>         id 
#> -0.2066936 
#> varStruct  parameters:
#>         power 
#> -0.0008370184 
#>  Beginning PNLS step: ..  completed fit_nlme() step.
#> PNLS step: RSS =  2909.464 
#>  fixed effects: -0.00304237  10.00692  0.001174652  2.721905  
#>  iterations: 2 
#> Convergence crit. (must all become <= tolerance = 1e-05):
#>       fixed    reStruct   varStruct 
#> 0.011205726 0.006916916 1.514605343 
#> 
#> **Iteration 3
#> LME step: Loglik: -4199.912, nlminb iterations: 8
#> reStruct  parameters:
#>         id 
#> -0.2066936 
#> varStruct  parameters:
#>        power 
#> -0.000836894 
#>  Beginning PNLS step: ..  completed fit_nlme() step.
#> PNLS step: RSS =  2909.464 
#>  fixed effects: -0.00304237  10.00692  0.001174652  2.721905  
#>  iterations: 1 
#> Convergence crit. (must all become <= tolerance = 1e-05):
#>        fixed     reStruct    varStruct 
#> 0.000000e+00 6.778246e-07 1.485866e-04 
#> 
#> **Iteration 4
#> LME step: Loglik: -4199.912, nlminb iterations: 8
#> reStruct  parameters:
#>         id 
#> -0.2066936 
#> varStruct  parameters:
#>        power 
#> -0.000836894 
#>  Beginning PNLS step: ..  completed fit_nlme() step.
#> PNLS step: RSS =  2909.464 
#>  fixed effects: -0.00304237  10.00692  0.001174652  2.721905  
#>  iterations: 1 
#> Convergence crit. (must all become <= tolerance = 1e-05):
#>        fixed     reStruct    varStruct 
#> 0.000000e+00 9.486459e-11 2.079497e-08 
# }