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:
sdcontrols the standard deviation when sampling from a normal distribution.sdis set to 1 by defaultalphacontrols the shape parameter for the Gamma distribution.alphais 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"
)
testdataNow, let’s fit a cglmm() model to this simulated dataset
to see how it matches with our original parameters:
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.mesorbeta.ampbeta.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) 1.08099168 0.10459679 0.87598574 1.28600 < 2.22e-16 ***
#> group1 -0.82359493 0.14792220 -1.11351711 -0.53367 2.5804e-08 ***
#> group0:main_rrr1 0.61114573 0.14125166 0.33429755 0.88799 1.5140e-05 ***
#> group1:main_rrr1 0.34310534 0.14125166 0.06625716 0.61995 0.015139 *
#> group0:main_sss1 2.08416386 0.15405462 1.78222234 2.38611 < 2.22e-16 ***
#> group1:main_sss1 0.11360193 0.15405462 -0.18833959 0.41554 0.460871
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Transformed coefficients:
#> estimate standard.error lower.CI upper.CI p.value
#> (Intercept) 1.08099168 0.10459679 0.87598574 1.28600 < 2.22e-16 ***
#> [group=1] -0.82359493 0.14792220 -1.11351711 -0.53367 2.5804e-08 ***
#> [group=0]:amp1 2.17192037 0.15174164 1.87451222 2.46933 < 2.22e-16 ***
#> [group=1]:amp1 0.36142312 0.14097782 0.08511166 0.63773 0.010357 *
#> [group=0]:acr1 1.28555925 0.06617809 1.15585256 1.41527 < 2.22e-16 ***
#> [group=1]:acr1 0.31973856 0.42693799 -0.51704452 1.15652 0.453911
#> ---
#> 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 = 1suggests that the intercept of the overall cosinor curve (accounting for all components) is 1amp=c(0.2,1)indicates that the amplitude of the first component is 0.2, and the second is 1acro = c(1.2, 2)indicates that the acrophase of the first component is 1.2 radians, and the second is 2 radiansperiod = 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) 6.999948875 0.003286322 6.993507802 7.00639 < 2.22e-16
#> group1 -2.605275628 0.017849548 -2.640260098 -2.57029 < 2.22e-16
#> group0:main_rrr1 0.060964424 0.004100189 0.052928201 0.06900 < 2.22e-16
#> group1:main_rrr1 1.071531618 0.018907831 1.034472950 1.10859 < 2.22e-16
#> group0:main_sss1 0.083939692 0.004267111 0.075576309 0.09230 < 2.22e-16
#> group1:main_sss1 1.720709358 0.023004689 1.675620996 1.76580 < 2.22e-16
#> group0:main_rrr2 0.027060904 0.004121213 0.018983474 0.03514 5.1599e-11
#> group1:main_rrr2 0.072661231 0.014254587 0.044722754 0.10060 3.4436e-07
#> group0:main_sss2 0.390298026 0.004144309 0.382175330 0.39842 < 2.22e-16
#> group1:main_sss2 -1.017352137 0.016788806 -1.050257592 -0.98445 < 2.22e-16
#> group0:main_rrr3 0.502647738 0.004482093 0.493862997 0.51143 < 2.22e-16
#> group1:main_rrr3 0.216612361 0.010117265 0.196782885 0.23644 < 2.22e-16
#> group0:main_sss3 0.042588094 0.004282384 0.034194777 0.05098 < 2.22e-16
#> group1:main_sss3 -0.321959466 0.012445071 -0.346351357 -0.29757 < 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) 6.999948875 0.003286322 6.993507802 7.00639 < 2.22e-16 ***
#> [group=1] -2.605275628 0.017849548 -2.640260098 -2.57029 < 2.22e-16 ***
#> [group=0]:amp1 0.103742628 0.003987590 0.095927096 0.11156 < 2.22e-16 ***
#> [group=1]:amp1 2.027071954 0.027144164 1.973870370 2.08027 < 2.22e-16 ***
#> [group=0]:amp2 0.391235021 0.004168415 0.383065077 0.39940 < 2.22e-16 ***
#> [group=1]:amp2 1.019943638 0.016953192 0.986715992 1.05317 < 2.22e-16 ***
#> [group=0]:amp3 0.504448703 0.004437890 0.495750599 0.51315 < 2.22e-16 ***
#> [group=1]:amp3 0.388044859 0.011022265 0.366441617 0.40965 < 2.22e-16 ***
#> [group=0]:acr1 0.942644147 0.042147747 0.860036081 1.02525 < 2.22e-16 ***
#> [group=1]:acr1 1.013833432 0.006040253 1.001994754 1.02567 < 2.22e-16 ***
#> [group=0]:acr2 1.501573158 0.010471530 1.481049336 1.52210 < 2.22e-16 ***
#> [group=1]:acr2 -1.499495492 0.013783782 -1.526511209 -1.47248 < 2.22e-16 ***
#> [group=0]:acr3 0.084525639 0.008580012 0.067709125 0.10134 < 2.22e-16 ***
#> [group=1]:acr3 -0.978563788 0.030025180 -1.037412060 -0.91972 < 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.97361085 0.02945959 0.91587111 1.03135 < 2.22e-16 ***
#> group1 1.03114478 0.04165514 0.94950221 1.11279 < 2.22e-16 ***
#> group0:main_rrr1 -0.47204160 0.04276800 -0.55586535 -0.38822 < 2.22e-16 ***
#> group1:main_rrr1 -0.04988532 0.04289070 -0.13394955 0.03418 0.244797
#> group0:main_sss1 0.05493568 0.04030572 -0.02406208 0.13393 0.172890
#> group1:main_sss1 0.24223254 0.04054894 0.16275808 0.32171 2.3174e-09 ***
#> group0:main_rrr2 0.04009558 0.04004698 -0.03839506 0.11859 0.316724
#> group1:main_rrr2 -0.23749686 0.03993414 -0.31576634 -0.15923 2.7275e-09 ***
#> group0:main_sss2 0.55014372 0.04314622 0.46557869 0.63471 < 2.22e-16 ***
#> group1:main_sss2 -0.03215631 0.04391460 -0.11822735 0.05391 0.464018
#> group0:main_rrr3 -0.57371693 0.04173341 -0.65551291 -0.49192 < 2.22e-16 ***
#> group1:main_rrr3 0.06936524 0.04093543 -0.01086674 0.14960 0.090170 .
#> group0:main_sss3 -0.08303979 0.04112501 -0.16364332 -0.00244 0.043466 *
#> group1:main_sss3 0.23529109 0.04258099 0.15183388 0.31875 3.2812e-08 ***
#> ---
#> 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.97361085 0.02945959 0.91587111 1.03135 < 2.22e-16 ***
#> [group=1] 1.03114478 0.04165514 0.94950221 1.11279 < 2.22e-16 ***
#> [group=0]:amp1 0.47522753 0.04305650 0.39083834 0.55962 < 2.22e-16 ***
#> [group=1]:amp1 0.24731589 0.04116544 0.16663311 0.32800 1.8800e-09 ***
#> [group=0]:amp2 0.55160291 0.04309634 0.46713563 0.63607 < 2.22e-16 ***
#> [group=1]:amp2 0.23966391 0.04002317 0.16121994 0.31811 2.1227e-09 ***
#> [group=0]:amp3 0.57969537 0.04178409 0.49780006 0.66159 < 2.22e-16 ***
#> [group=1]:amp3 0.24530274 0.04321163 0.16060950 0.33000 1.3726e-08 ***
#> [group=0]:acr1 3.02573495 0.08416471 2.86077515 3.19069 < 2.22e-16 ***
#> [group=1]:acr1 1.77389663 0.17103370 1.43867674 2.10912 < 2.22e-16 ***
#> [group=0]:acr2 1.49804295 0.07269839 1.35555672 1.64053 < 2.22e-16 ***
#> [group=1]:acr2 -3.00701426 0.18289562 -3.36548308 -2.64855 < 2.22e-16 ***
#> [group=0]:acr3 -2.99785089 0.07085361 -3.13672142 -2.85898 < 2.22e-16 ***
#> [group=1]:acr3 1.28411142 0.16416108 0.96236162 1.60586 5.1881e-15 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
autoplot(object,
superimpose.data = TRUE,
x_str = "group",
predict.ribbon = FALSE
)