Simulating data

library(GLMMcosinor)

Simulating rhythmic data

simulate_cosinor allows users to simulate circadian data from Gaussian, Gamma, Binomial, or Poisson distributions. It also supports generation of multi-component data, as well as simulation of grouped data with two levels.

Understanding the inputs for a simple model

n specifies the number of datapoints.

mesor, amp, and acro represent the parameters that will be used to simulate the dataset. Note that acro should be expressed in units of radians.

period determines the period of the dataset

n_components corresponds to the number of components in the simulated dataset. Details about how to specify a multi-component model are included later in this vignette

The family argument determines the distribution that the data is simulated from. Currently, simulate_cosinor supports simulations from Gaussian, Gamma, Binomial, and Poisson distributions:

  • family = 'gaussian'

  • family = 'gamma'

  • family = 'binomial'

  • family = 'poisson'

Note that the … parameter controls extra arguments such as standard deviation, and the shape parameter for a Gamma distribution:

  • sd controls the standard deviation when sampling from a normal distribution. sd is set to 1 by default

  • alpha controls the shape parameter for the Gamma distribution. alpha is set to 1 by default

n_period is the number of periods that are simulated. By default, the maximum period supplied defines the upper limit of the time vector used in the simulation. Thus, increasing n_period increases the number of cycles that are simulated.

Consider the following example of a single-component Poisson data-set with no grouping variable:

testdata <- simulate_cosinor(
  n = 200,
  mesor = 1,
  amp = 2,
  acro = 1.2,
  period = 12,
  n_period = 3,
  family = "poisson"
)

testdata

Now, let’s fit a cglmm() model to this simulated dataset to see how it matches with our original parameters:

object <- cglmm(
  Y ~ amp_acro(times,
    n_components = 1,
    period = 12
  ),
  data = testdata,
  family = poisson()
)
summary(object)
autoplot(object, superimpose.data = TRUE)

Simulating grouped cosinor data

The simulate_cosinor() function can simulate grouped data from two levels with their own parameters when beta.group = TRUE. The reference group is specified by the parameters mesor, amp, acro. The treatment group is specified in the same manner, but with the beta prefix. For example:

  • beta.mesor

  • beta.amp

  • beta.acro

Similarly, the standard deviation (sd) for the Gaussian distribution or the alpha parameter for the Gamma distribution are specified as:

  • beta.sd (1 by default)

  • beta.alpha (1 by default)

simulate_cosinor(
  n = 100,
  mesor = 1,
  amp = 2,
  acro = 1.2,
  period = 12,
  beta.group = TRUE,
  beta.mesor = 0.4,
  beta.amp = 0.5,
  beta.acro = 0.2,
  n_period = 3,
  n_components = 1,
  family = "poisson"
)
#> 
#>  Conditional Model 
#> Raw model coefficients:
#>                     estimate standard.error    lower.CI upper.CI    p.value    
#> (Intercept)       0.99064873     0.09024656  0.81376872  1.16753 < 2.22e-16 ***
#> group1           -0.68848245     0.12762791 -0.93862856 -0.43834 6.8734e-08 ***
#> group0:main_rrr1  0.85702693     0.12924737  0.60370673  1.11035 3.3364e-11 ***
#> group1:main_rrr1  0.58411241     0.12924737  0.33079221  0.83743 6.2034e-06 ***
#> group0:main_sss1  2.12764925     0.12633674  1.88003378  2.37526 < 2.22e-16 ***
#> group1:main_sss1  0.13867025     0.12633674 -0.10894522  0.38629    0.27237    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Transformed coefficients:
#>                   estimate standard.error    lower.CI upper.CI    p.value    
#> (Intercept)     0.99064873     0.09024656  0.81376872  1.16753 < 2.22e-16 ***
#> [group=1]      -0.68848245     0.12762791 -0.93862856 -0.43834 6.8734e-08 ***
#> [group=0]:amp1  2.29377123     0.12308361  2.05253179  2.53501 < 2.22e-16 ***
#> [group=1]:amp1  0.60034719     0.12677419  0.35187434  0.84882 2.1845e-06 ***
#> [group=0]:acr1  1.18787452     0.05769934  1.07478589  1.30096 < 2.22e-16 ***
#> [group=1]:acr1  0.23308831     0.21457304 -0.18746712  0.65364    0.27735    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Simulating multi-component cosinor data

To simulate multi-component data with the number of components corresponding to n_components , specify a vector of values for the parameter inputs, amp, acro, and period. Since only one mesor is provided for a multi-component cosinor curve, the mesor argument in simualte_cosinor will only accept single-element inputs. For example:

testdata <- simulate_cosinor(
  n = 200,
  mesor = 1,
  amp = c(0.2, 1),
  acro = c(1.2, 2),
  period = c(12, 6),
  n_components = 2,
  n_period = 2,
  family = "poisson"
)

testdata
object <- cglmm(
  Y ~ amp_acro(times,
    n_components = 2,
    period = c(12, 6)
  ),
  data = testdata,
  family = poisson()
)
summary(object)
autoplot(object, superimpose.data = TRUE)

In this example:

  • mesor = 1 suggests that the intercept of the overall cosinor curve (accounting for all components) is 1

  • amp=c(0.2,1) indicates that the amplitude of the first component is 0.2, and the second is 1

  • acro = c(1.2, 2) indicates that the acrophase of the first component is 1.2 radians, and the second is 2 radians

  • period = c(12, 6) indicates that the period of the first component is 12 units, and the second is 6 units.

  • n_components = 2, because we are specifying a two-component dataset

To simulate a dataset with more than two components, specify more elements in the vector inputs for these parameters. Ensure that the number of inputs for each parameter corresponds to the number of components specified in n_components. For example, if n_components = 3, then amp, acro, period must all have 3 elements corresponding to each of the three components.

The following are examples of a multi-component cosinor dataset with multiple groups. The first one is from a Poisson distribution, and the second is from a Gamma distribution

testdata <- simulate_cosinor(100,
  mesor = 7,
  amp = c(0.1, 0.4, 0.5),
  acro = c(1, 1.5, 0.1),
  beta.mesor = 4.4,
  beta.amp = c(2, 1, 0.4),
  beta.acro = c(1, -1.5, -1),
  family = "poisson",
  period = c(12, 6, 8),
  n_period = 2,
  n_components = 3
)

object <- cglmm(Y ~ group + amp_acro(times,
  n_components = 3,
  period = c(12, 6, 8),
  group = "group"
), data = testdata, family = poisson())
summary(object)
#> 
#>  Conditional Model 
#> Raw model coefficients:
#>                      estimate standard.error     lower.CI upper.CI    p.value
#> (Intercept)       7.004555472    0.003172728  6.998337039  7.01077 < 2.22e-16
#> group1           -2.602075479    0.017823662 -2.637009215 -2.56714 < 2.22e-16
#> group0:main_rrr1  0.054993419    0.004467515  0.046237251  0.06375 < 2.22e-16
#> group1:main_rrr1  1.083751367    0.019775924  1.044991268  1.12251 < 2.22e-16
#> group0:main_sss1  0.085732534    0.004274174  0.077355308  0.09411 < 2.22e-16
#> group1:main_sss1  1.661344858    0.022593537  1.617062338  1.70563 < 2.22e-16
#> group0:main_rrr2  0.029084034    0.004111333  0.021025969  0.03714 1.5041e-12
#> group1:main_rrr2  0.057735602    0.014159539  0.029983415  0.08549 4.5521e-05
#> group0:main_sss2  0.400739017    0.004521904  0.391876248  0.40960 < 2.22e-16
#> group1:main_sss2 -0.994410152    0.017340552 -1.028397010 -0.96042 < 2.22e-16
#> group0:main_rrr3  0.498451936    0.004352240  0.489921703  0.50698 < 2.22e-16
#> group1:main_rrr3  0.236123734    0.009842052  0.216833668  0.25541 < 2.22e-16
#> group0:main_sss3  0.053267902    0.004331736  0.044777854  0.06176 < 2.22e-16
#> group1:main_sss3 -0.346160630    0.013541908 -0.372702281 -0.31962 < 2.22e-16
#>                     
#> (Intercept)      ***
#> group1           ***
#> group0:main_rrr1 ***
#> group1:main_rrr1 ***
#> group0:main_sss1 ***
#> group1:main_sss1 ***
#> group0:main_rrr2 ***
#> group1:main_rrr2 ***
#> group0:main_sss2 ***
#> group1:main_sss2 ***
#> group0:main_rrr3 ***
#> group1:main_rrr3 ***
#> group0:main_sss3 ***
#> group1:main_sss3 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Transformed coefficients:
#>                    estimate standard.error     lower.CI upper.CI    p.value    
#> (Intercept)     7.004555472    0.003172728  6.998337039  7.01077 < 2.22e-16 ***
#> [group=1]      -2.602075479    0.017823662 -2.637009215 -2.56714 < 2.22e-16 ***
#> [group=0]:amp1  0.101854522    0.003943609  0.094125190  0.10958 < 2.22e-16 ***
#> [group=1]:amp1  1.983578524    0.026402665  1.931830251  2.03533 < 2.22e-16 ***
#> [group=0]:amp2  0.401793032    0.004535448  0.392903717  0.41068 < 2.22e-16 ***
#> [group=1]:amp2  0.996084811    0.017377572  0.962025395  1.03014 < 2.22e-16 ***
#> [group=0]:amp3  0.501290137    0.004335264  0.492793176  0.50979 < 2.22e-16 ***
#> [group=1]:amp3  0.419024581    0.011653132  0.396184862  0.44186 < 2.22e-16 ***
#> [group=0]:acr1  1.000452796    0.046751465  0.908821608  1.09208 < 2.22e-16 ***
#> [group=1]:acr1  0.992782238    0.007208565  0.978653711  1.00691 < 2.22e-16 ***
#> [group=0]:acr2  1.498347353    0.010195266  1.478364999  1.51833 < 2.22e-16 ***
#> [group=1]:acr2 -1.512801286    0.014169557 -1.540573108 -1.48503 < 2.22e-16 ***
#> [group=0]:acr3  0.106462619    0.008675068  0.089459799  0.12347 < 2.22e-16 ***
#> [group=1]:acr3 -0.972170153    0.028683022 -1.028387844 -0.91595 < 2.22e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
autoplot(object,
  superimpose.data = TRUE,
  x_str = "group",
  predict.ribbon = FALSE
)

testdata <- simulate_cosinor(500,
  mesor = 1,
  amp = c(0.5, 0.5, 0.5),
  acro = c(pi, pi / 2, pi),
  alpha = 2,
  beta.mesor = 2,
  beta.amp = c(0.2, 0.2, 0.2),
  beta.acro = c(pi / 2, pi, pi / 2),
  beta.alpha = 3,
  family = "gamma",
  period = c(12, 6, 8),
  n_period = 2,
  n_components = 3
)

object <- cglmm(Y ~ group + amp_acro(times,
  n_components = 3,
  period = c(12, 6, 8),
  group = "group"
), data = testdata, family = Gamma(link = "log"))
summary(object)
#> 
#>  Conditional Model 
#> Raw model coefficients:
#>                     estimate standard.error    lower.CI upper.CI    p.value    
#> (Intercept)       0.95997751     0.02836680  0.90437960  1.01558 < 2.22e-16 ***
#> group1            1.03396082     0.04011306  0.95534067  1.11258 < 2.22e-16 ***
#> group0:main_rrr1 -0.52265812     0.03815017 -0.59743108 -0.44789 < 2.22e-16 ***
#> group1:main_rrr1  0.05088299     0.03816745 -0.02392384  0.12569    0.18248    
#> group0:main_sss1  0.07016562     0.04179883 -0.01175857  0.15209    0.09322 .  
#> group1:main_sss1  0.21675887     0.04182011  0.13479297  0.29872 2.1820e-07 ***
#> group0:main_rrr2 -0.03505504     0.04075876 -0.11494073  0.04483    0.38976    
#> group1:main_rrr2 -0.21488726     0.03970967 -0.29271678 -0.13706 6.2513e-08 ***
#> group0:main_sss2  0.46166507     0.03921045  0.38481399  0.53852 < 2.22e-16 ***
#> group1:main_sss2 -0.03327023     0.04042366 -0.11249915  0.04596    0.41049    
#> group0:main_rrr3 -0.53016555     0.04035097 -0.60925201 -0.45108 < 2.22e-16 ***
#> group1:main_rrr3 -0.01546651     0.04046687 -0.09478012  0.06385    0.70231    
#> group0:main_sss3  0.04200748     0.03937756 -0.03517112  0.11919    0.28607    
#> group1:main_sss3  0.20325492     0.03943577  0.12596224  0.28055 2.5489e-07 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Transformed coefficients:
#>                   estimate standard.error    lower.CI upper.CI    p.value    
#> (Intercept)     0.95997751     0.02836680  0.90437960  1.01558 < 2.22e-16 ***
#> [group=1]       1.03396082     0.04011306  0.95534067  1.11258 < 2.22e-16 ***
#> [group=0]:amp1  0.52734687     0.03836159  0.45215955  0.60253 < 2.22e-16 ***
#> [group=1]:amp1  0.22265104     0.04142205  0.14146532  0.30384 7.6505e-08 ***
#> [group=0]:amp2  0.46299405     0.03925152  0.38606249  0.53993 < 2.22e-16 ***
#> [group=1]:amp2  0.21744757     0.03992544  0.13919515  0.29570 5.1416e-08 ***
#> [group=0]:amp3  0.53182717     0.04024444  0.45294952  0.61070 < 2.22e-16 ***
#> [group=1]:amp3  0.20384253     0.03931350  0.12678949  0.28090 2.1595e-07 ***
#> [group=0]:acr1  3.00814289     0.07889472  2.85351209  3.16277 < 2.22e-16 ***
#> [group=1]:acr1  1.34022636     0.17336135  1.00044435  1.68001 1.0685e-14 ***
#> [group=0]:acr2  1.64658265     0.08794759  1.47420854  1.81896 < 2.22e-16 ***
#> [group=1]:acr2 -2.98798584     0.18492073 -3.35042381 -2.62555 < 2.22e-16 ***
#> [group=0]:acr3  3.06252320     0.07424674  2.91700227  3.20804 < 2.22e-16 ***
#> [group=1]:acr3  1.64674410     0.19910305  1.25650930  2.03698 < 2.22e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
autoplot(object,
  superimpose.data = TRUE,
  x_str = "group",
  predict.ribbon = FALSE
)