Human Circadian Rhythms

Medical studies often collect physiological and/or psychological measurements over time from multiple subjects in order to study dynamics such as circadian rhythms and interactions between variables. Experiments are typically conducted in such a way that variables of interest are measured several times during a time period, e.g. 24 hours, from a group of normal (or ill) human subjects. The problems of interest are: (1) do circadian rhythms exist? (2) do demographic, environmental and psychological variables affect circadian rhythms, and if so, how? and (3) how are variables associated with each other?

In an experiment to study immunological responses in humans, blood samples were collected every two hours for 24 hours from 9 healthy normal volunteers, 11 patients with major depression and 16 patients with Cushing's disease. These blood samples were analyzed for parameters that measure immune functions and hormones of the HPA axis (Kronfol et al., 1997). We will concentrate on one hormone: cortisol (cort).The data frame horm.cort contains the following variables: ID of subjects, time when measurements are taken, hormone measurements conc, and subject type (normal, depression, cushing). For simplicity, the time variable of 24-hour period is transformed into [0,1].

Observations are shown in Figures , and . The data are noisy and it is difficult to identify patterns among subjects. The investigator hypothesizes that there is a common curve for all individuals. However, the time axis may be shifted and the magnitude of the values may differ between subjects; i.e., there may be phase and amplitude differences between subjects. Since the 24-hour periodicity is entrained, the cycle length is fixed. The common practice is to fit a sinusoidal function to each subject. Problems with this approach are: (1) the pattern over time may not be symmetric. That is, the peak and nadir may not be separated by 12 hours and/or the amplitude and width of the peak may differ from those of the nadir; (2) sometimes there are local minimum and maximum points (Wang and Brown, 1996). Although adding harmonics can improve the fit, it is difficult to decide how many harmonics to include in the model and the results are difficult to interpret.

We first show how to fit a SIM to this kind of data. Consider the following model (Wang and Brown, 1996)

$$y_{ij} = \phi_{1i} + \exp (\phi_{2i}) f(t_{ij}-\exp (\phi_{3i})/(1+\exp (\phi_{3i}))) + \epsilon_{ij},$$

$$i=1, \cdots, m, ~~ j=1, \cdots, n_i,$$ (95)

where $m$ is the total number of subjects, $n_i$ is the number of observations for subject $i$, $y_{ij}$ is the response of $i$th individual at the $j$th time point $t_{ij}$, $\phi_{1i}$ is the 24-hour mean of the $i$th individual, $\exp (\phi_{2i})$ is the amplitude of the $i$th individual where exponential transformation is used to enforce positive constraint, $\exp (\phi_{3i})/(1+\exp (\phi_{3i}))$ is the phase of the $i$th individual where the inverse logistic transformation forces the phase inside the interval [0,1]. $\epsilon_{ij}$ are random errors and $\epsilon_{ij} \stackrel{iid}{\sim} \mbox{N}(0,\sigma^2)$. The function $f$ is the common curve. Since it is a periodic function, we use periodic spline to model $f$. In order to make $f$ identifiable with $\phi_{1i}$, we use the side condition $\int_0^1 f=0$ which is equivalent to removing the constant from the model space. Thus we assume that $f \in W^0_2(per)$. In order to make $\phi_{2i}$ and $\phi_{3i}$ identifiable with $f$, we add constraints: $\phi_{21}=\phi_{31}=0$.

We now fit model () to cortisol measurements from normal subjects.

> options(contrasts=c("contr.treatment", "contr.poly"))
> data(horm.cort)
> cort.nor <- horm.cort[horm.cort$type=="normal",] 
> M <- model.matrix(~as.factor(ID), data=cort.nor)[,-1]
> cort.nor.snr.fit1 <- snr(conc~b1+exp(b2)*f(time-alogit(b3)), 
    func=f(u)~list(periodic(u)),params=list(b1~as.factor(ID), b2+b3~M-1),
    start=list(params=c(mean(cort.nor$conc),rep(0,24))), data=cort.nor, 
    spar="m", control=list(prec.out=0.001,converg="PRSS"))
> summary(cort.nor.snr.fit1)

Semi-parametric Nonlinear Regression Model Fit
  Model: horm ~ b1 + exp(b2) * f(time - alogit(b3)) 
  Data: cort.nor.dat 

       AIC     BIC   logLik 
  113.9554 183.449 -30.9777

Coefficients:
...

Standardized residuals:
      Min         Q1        Med        Q3      Max 
 -2.11257 -0.5740086 0.04326369 0.5185883 2.521769

GML estimate(s) of smoothing spline parameter(s): 1.505551e-06 
Equivalent Degrees of Freedom (DF) for spline function:  10.13926 
Degrees of freedom: 107 total; 71.86074 residual
Residual standard error: 0.4213113 

Converged after 5 iterations

Notice that we used the option converg="PRSS" instead of the default converg="COEF" because this option usually requires fewer number of iterations (the same is true for the snm function). We also relaxed the tolerance for convergence criterion. The potential risk of these options is that the algorithm may stop too early. The fits shown in Figure seems adequate.

Figure: Plots of the data and fitted curves. Circles are observations. Solid lines represent fits from cort.nor.snr.fit.1. Subjects' ID are shown in the strip.

Random errors in model () may be correlated. In the following we fit with an AR(1) within-subject correlation structure.

> cort.nor.snr.fit2 <- update(cort.nor.snr.fit1, cor=corAR1(form=~1|ID))
> summary(cort.nor.snr.fit2)
Semi-parametric Nonlinear Regression Model Fit
Model: conc ~ b1 + exp(b2) * f(time - alogit(b3)) 
Data: cort.nor 

       AIC      BIC    logLik
  117.6520 189.8184 -31.82602

Correlation Structure: AR(1)
 Formula: ~1 | ID 
 Parameter estimate(s):
       Phi 
-0.1641137 

Coefficients:
...

GML estimate(s) of smoothing spline parameter(s): 0.0001425255 
Equivalent Degrees of Freedom (DF) for spline function:  10.33217 
Residual standard error: 0.4311652 

Converged after 5 iterations

The lag 1 autocorrelation coefficient is small.

Parameters $\mbox{\boldmath$\phi$}_{1i}$, $\mbox{\boldmath$\phi$}_{2i}$ and $\mbox{\boldmath$\phi$}_{3i}$ in model () are deterministic. Thus model () has following drawbacks: (1) they ignore correlations among observations; (2) they have the number of parameters proportional to the number of subjects; and (3) it is difficult to investigate covariate effects on parameters and/or the common curve. In the remaining of this section we show how to fit hormone data using SNM models.

We first show how to fit individual hormone for each group. Consider the following model

$$y_{ij} = \mu + b_{1i} + \exp (b_{2i}) f(t_{ij}-\exp (b_{3i})/(1+\exp (b_{3i}))) + \epsilon_{ij},$$

$$i=1, \cdots, m, ~~ j=1, \cdots, n_i,$$ (96)

where the fixed effect $\mu$ represents 24-hour mean of the population, the random effects $b_{1i}$, $b_{2i}$ and $b_{3i}$ represent the $i$th subject's deviation of 24-hour mean, amplitude and phase. We assume that $f \in W^0_2(per)$ and $\mbox{\boldmath$b$}_i=(b_{1i},b_{2i},b_{3i})^T \stackrel{iid}{\sim} N(0,\sigma^2 \mbox{\boldmath$D$})$, where $\mbox{\boldmath$D$}$ is an unstructured positive-definite matrix. The assumption of zero population mean for amplitude and phase takes care of potential confounding between amplitude, phase and the nonparametric common function $f$ in a natural way.

We now fit model () to cortisol measurements from normal subjects.

> cort.nor.fit1 <- snm(conc~b1+exp(b2)*f(time-alogit(b3)), 
    func=f(u)~list(periodic(u)), data=cort.nor, fixed=list(b1~1), 
    random=list(b1+b2+b3~1), start=c(mean(cort.nor$conc)), groups=~ID, 
    spar="m", control=list(prec.out=0.005,converg="PRSS"))
> summary(cort.nor.fit1)
Semi-parametric Nonlinear Mixed Effects Model fit
  Model: conc ~ b1 + exp(b2) * f(time - alogit(b3)) 
  Data: cort.nor 

       AIC      BIC   logLik
  176.4104 224.4590 -70.2004

Random effects:
 Formula: list(b1 ~ 1, b2 ~ 1, b3 ~ 1)
 Level: ID
 Structure: General positive-definite, Log-Cholesky parametrization
         StdDev    Corr         
b1       0.2475622 b1     b2    
b2       0.1889174 -0.629       
b3       0.2483884  0.035 -0.479
Residual 0.3948469              

Fixed effects: list(b1 ~ 1) 
      Value  Std.Error DF  t-value p-value
b1 1.670316 0.07686368 98 21.73089  <.0001

GML estimate(s) of smoothing parameter(s): 0.0001218930 
Equivalent Degrees of Freedom (DF):  10.00480 

Converged after 4 iterations

> cort.dep <- horm.cort[horm.cort$type=="depression",] 
> cort.dep.fit1 <- snm(conc~b1+exp(b2)*f(time-alogit(b3)), 
    func=f(u)~list(periodic(u)), data=cort.dep, fixed=list(b1~1), 
    random=list(b1+b2+b3~1), start=c(mean(cort.dep$conc)), groups=~ID, 
    spar="m", control=list(prec.out=0.005,converg="PRSS"))
> summary(cort.dep.fit1)

Semi-parametric Nonlinear Mixed Effects Model fit
 Model: conc ~ b1 + exp(b2) * f(time - alogit(b3)) 
  Data: cort.dep 

       AIC      BIC    logLik
  248.2489 293.5048 -108.4047

Random effects:
 Formula: list(b1 ~ 1, b2 ~ 1, b3 ~ 1)
 Level: ID
 Structure: General positive-definite, Log-Cholesky parametrization
         StdDev    Corr         
b1       0.4139559 b1     b2    
b2       0.3883549 -0.815       
b3       0.3806670  0.195 -0.099
Residual 0.4390506              

Fixed effects: list(b1 ~ 1) 
      Value  Std.Error  DF  t-value p-value
b1 1.956112 0.09046153 121 21.62369  <.0001

GML estimate(s) of smoothing parameter(s): 0.0005380155 
Equivalent Degrees of Freedom (DF):  7.719702 

Converged after 5 iterations

> cort.cush <- horm.cort[horm.cort$type=="cushing",] 
> cort.cush.fit1 <- snm(conc~b1+exp(b2)*f(time-alogit(b3)), 
    func=f(u)~list(periodic(u)), data=cort.cush, fixed=list(b1~1), 
    random=list(b1+b2+b3~1), start=c(mean(cort.cush$conc)), groups=~ID, 
    spar="m", control=list(prec.out=0.005,converg="PRSS"))
> summary(cort.cush.fit1)
  Model: conc ~ b1 + exp(b2) * f(time - alogit(b3)) 
  Data: cort.cush 

        AIC       BIC   logLik
  -38.94984 -13.18616 27.47518

Random effects:
 Formula: list(b1 ~ 1, b2 ~ 1, b3 ~ 1)
 Level: ID
 Structure: General positive-definite, Log-Cholesky parametrization
         StdDev       Corr         
b1       3.429996e-01 b1     b2    
b2       1.028451e+04 -0.690       
b3       6.520918e+03 -0.107  0.593
Residual 1.685278e-01              

Fixed effects: list(b1 ~ 1) 
      Value  Std.Error  DF t-value p-value
b1 3.034334 0.07610914 170 39.8682  <.0001

GML estimate(s) of smoothing parameter(s): 999.9996 
Equivalent Degrees of Freedom (DF):  0.0002583048 

Converged after 2 iterations

The fits are shown in Figures , and .

Figure: Plots of the data and fitted curves. Circles are observations. Solid lines represent fits from cort.nordep.fit.3. Dotted lines represent fits from cort.nor.fit.1. Subjects' ID are shown in the strip.

Figure: Plots of the data and fitted curves. Circles are observations. Solid lines represent fits from cort.nordep.fit.3. Dotted lines represent fits from cort.dep.fit.1. Subjects' ID are shown in the strip.

Figure: Plots of the data and fitted curves. Circles are observations. Solid lines represent fits from cort.cush.fit.1. Subjects' ID are shown in the strip.

We calculate the posterior means and variances of the common function $f$ for all three groups:

u <- seq(0,1,len=50)
cort.nor.ci <- intervals(cort.nor.fit.1, newdata=data.frame(u=u))
cort.dep.ci <- intervals(cort.dep.fit.1, newdata=data.frame(u=u))
cort.cush.ci <- intervals(cort.cush.fit.1, newdata=data.frame(u=u))

The estimated common functions are shown in Figure together with their 95% Bayesian confidence intervals.

Figure: Solid lines are estimates of the common functions and dotted lines are 95% Bayesian confidence intervals. The left three panels are estimates from cort.nor.fit.1, cort.dep.fit.1 and cort.cush.fit.1 respectively. The right panel is the estimate from cort.nordep.fit.3.

It is obvious that the common function for Cushing group is almost zero, which suggests that circadian rhythms are lost for Cushing patients except subjects 3044 and 3069. It seems that the shape functions for normal and depression groups are similar. We now test this hypothesis by fitting data from these two groups simultaneous. Consider the following model

$$y_{ijk} = \mu_k + b_{1ik} + \exp (b_{2ik}) f(k,t_{ijk}-\exp (b_{3ik})/(1+\exp (b_{3ik}))) + \epsilon_{ijk},$$

$$i=1, \cdots, m, ~~j=1, \cdots, n_i, ~~k=1,2,$$ (97)

where $k$ represents group with $k=1$ and $k=2$ correspond to depression and normal groups respectively, fixed effects $\mu_k$ is the population 24-hour mean of group $k$, random effects $b_{1ik}$, $b_{2ik}$ and $b_{3ik}$ represent the $i$th subject's deviation of 24-hour mean, amplitude and phase. Note that subjects are nested within group, even though this is not made explicit in out notations. We allow different correlation structures for the random effects in each group. That is, we assume that $\mbox{\boldmath$b$}_{ik}=(b_{1ik},b_{2ik},b_{3ik})^T \stackrel{iid}{\sim} N(0,\sigma^2 \mbox{\boldmath$D$}_k)$, where $\mbox{\boldmath$D$}_k$'s are unstructured positive-definite matrices. We assume different common functions for each group. Thus $f$ is a function of both group (denoted as $k$) and time (denoted as $t$). We model group effect using an one-way ANOVA effect model and time effect using a periodic spline model without the constant term. That is, we assume that $f \in R^2 \otimes W^0_2(per)$. Writing $R^2=\{ 1 \} \oplus \{ g:~\sum_{k=1}^2 g(k) = 0 \}$, we have the following SS ANOVA decomposition

$$R^2 \otimes W^0_2(per) = W^0_2(per) \oplus \{ g:~\sum_{k=1}^2 g(k) = 0 \} \otimes W^0_2(per),$$

or equivalently,

$$f(k,t)=s(t)+ss(k,t),$$ (98)

where $s(t)$ is the main effect of time, and $ss(k,t)$ is the interaction between group and time. The hypothesis $H_0:~f(1,t)=f(2,t)$ is equivalent to $H_0:~ss(k,t)=0$. We now fit model (). Note that pdStrat is not available currently for the R version of nlme. The following results were based on the old S-Plus version.

> cort.nordep.dat <- horm[horm$type=="cort"&horm$group!="cushing",] 
> cort.nordep.dat$group <- as.vector(cort.nordep.dat$group)
> cort.nordep.fit.1 <- snm(horm~b1+exp(b2)*f(group,time-alogit(b3)), 
                         func=f(g,u)~list(list(periodic(u),
                         rk.prod(shrink1(g),periodic(u)))),
                         data=cort.nordep.dat, fixed=list(b1~group),
                         random=pdStrat(b1+b2+b3~1,strata=~as.factor(group)),	
                         weights=varIdent(form=~1|group), 
                         control=list(prec.out=0.005,converg="PRSS"),	
                         spar="m", start=c(1.8,-.2), groups=~ID)
> summary(cort.nordep.fit.1)
Semi-parametric Nonlinear Mixed Effects Model fit
  Model: cort ~ b1 + exp(b2) * f(group, time - alogit(b3)) 
  Data: cort.nordep.dat 

       AIC      BIC    logLik 
  430.7103 516.9962 -190.4965

Random effects:
 Formula: list(b1 ~ 1, b2 ~ 1, b3 ~ 1)
 Level: ID
 Structure: General positive-definite stratified by as.factor(group)
Stratum: depression 
                  StdDev   Corr        
b1.(Intercept) 0.4025420 b1.(I) b2    
            b2 0.3441774 -0.779       
            b3 0.3352246  0.080 -0.013
Stratum: normal 
                  StdDev   Corr        
b1.(Intercept) 0.2526908 b1.(I) b2    
            b2 0.2352516 -0.535       
            b3 0.2431394 -0.050 -0.384
 Within-group standard deviation: 0.4516581

Variance function:
 Structure: Different standard deviations per stratum
 Formula:  ~ 1 | group 
 Parameter estimates:
 depression    normal 
          1 0.9344799
Fixed effects: list(b1 ~ group) 
                   Value Std.Error  DF   t-value p-value 
b1.(Intercept)  1.860474 0.0962724 218  19.32511  <.0001
      b1.group -0.160235 0.1274092 218  -1.25764  0.2099
 Correlation: 
         b1.(I) 
b1.group -0.756

GML estimate(s) of smoothing parameter(s): 5.037821e-04 1.487554e+02 
Equivalent Degrees of Freedom (DF):  8.858698 

Converged after 4 iterations

> u <- seq(0,1,len=50)
> cort.nordep.ci <- intervals(cort.nordep.fit.1, newdata=data.frame(
                  g=rep(c("normal","depression"),c(50,50)),u=rep(u,2)),
                  terms=c(0,1))

The smoothing parameter for the interaction term $ss(k,t)$ is large, which means that the interaction is small. In fact, $ss(k,t)$ is essentially zero: the posterior means are on the magnitude of $10^{-6}$ while the posterior standard deviations are on the magnitude of $10^{-4}$. Therefore normal and depression groups have the same shape function.

Under the assumption of one shape function for both groups, we now can investigate differences of 24-hour mean, amplitude, and phase between two groups. For this purpose, consider the following model

$$y_{ijk} = \mu_k + b_{1ik} + \exp (b_{2ik}+d_1 \times I(k=2)) \times$$

$$f(t_{ijk}-\exp (b_{3ik}+d_2 \times I(k=2))/(1+\exp (b_{3ik}+d_2 \times I(k=2)))) + \epsilon_{ijk},$$

$$i=1, \cdots, m, ~~j=1, \cdots, n_i, ~~k=1,2,$$ (99)

where $d_1$ and $d_2$ measures the differences of amplitude and phase between normal and depression groups.

 
> cort.nordep.fit.2 <- snm(horm~b1+exp(b2+d1*I(group=="normal"))*
                        f(time-alogit(b3+d2*I(group=="normal"))), 
                        func=f(u)~list(periodic(u)),
                        data=cort.nordep.dat,
                        fixed=list(b1~group,d1+d2~1), 
                        random=pdStrat(b1+b2+b3~1,strata=~as.factor(group)),
                        weights=varIdent(form=~1|group), 
                        control=list(prec.out=0.005,converg="PRSS"),
                        spar="m", start=c(1.9,-0.3,0,0), groups=~ID)
> summary(cort.nordep.fit.2)
Semi-parametric Nonlinear Mixed Effects Model fit
  Model: cort ~ b1 + exp(b2 + d1 * I(group == "normal")) * 
            f(time - alogit(b3 + d2 * I(group == "normal"))) 
  Data: cort.nordep.dat 

       AIC      BIC    logLik 
  438.1875 531.3631 -192.2045

Random effects:
 Formula: list(b1 ~ 1, b2 ~ 1, b3 ~ 1)
 Level: ID
 Structure: General positive-definite stratified by as.factor(group)
Stratum: depression 
                  StdDev   Corr        
b1.(Intercept) 0.4072866 b1.(I) b2    
            b2 0.3789050 -0.784       
            b3 0.3324878  0.086 -0.046
Stratum: normal 
                  StdDev   Corr        
b1.(Intercept) 0.2448359 b1.(I) b2    
            b2 0.1626245 -0.522       
            b3 0.2661856 -0.032 -0.640
 Within-group standard deviation: 0.4515871

Variance function:
 Structure: Different standard deviations per stratum
 Formula:  ~ 1 | group 
 Parameter estimates:
 depression    normal 
          1 0.9358295
Fixed effects: list(b1 ~ group, d1 + d2 ~ 1) 
                   Value Std.Error  DF   t-value p-value 
b1.(Intercept)  1.907384 0.0956961 216  19.93168  <.0001
      b1.group -0.272381 0.1310609 216  -2.07828  0.0389
            d1  0.234996 0.0766247 216   3.06684  0.0024
            d2  0.029918 0.0915166 216   0.32691  0.7441
 Correlation: 
         b1.(I) b1.typ     d1 
b1.group -0.730              
     d1   0.000 -0.225       
     d2   0.000 -0.017 -0.428

GML estimate(s) of smoothing parameter(s): 0.0005858115 
Equivalent Degrees of Freedom (DF):  8.889222 

Converged after 4 iterations

The differences of 24-hour mean and amplitude are significant, while the difference of phase is not. We refit without the $d_2$ term.

> acth.nordep.fit.3 <- snm(horm~b1+exp(b2+d1*I(group=="normal"))*
                        f(time-alogit(b3)), 
                        func=f(u)~list(periodic(u)),
                        data=acth.nordep.dat,
                        fixed=list(b1~group,d1~1), 
                        random=pdStrat(b1+b2+b3~1,strata=~as.factor(group)),
                        weights=varIdent(form=~1|group), 
                        control=list(prec.out=0.005,converg="PRSS"),
                        spar="m", start=c(2.8,-0.4,0), groups=~ID)
> summary(acth.nordep.fit.3)
Semi-parametric Nonlinear Mixed Effects Model fit
  Model: acth ~ b1 + exp(b2 + d1 * I(group == "normal")) * 
                f(time - alogit(b3)) 
  Data: acth.nordep.dat 

       AIC      BIC    logLik 
  417.8258 518.0793 -180.2195

Random effects:
 Formula: list(b1 ~ 1, b2 ~ 1, b3 ~ 1)
 Level: ID
 Structure: General positive-definite stratified by as.factor(group)
Stratum: depression 
                  StdDev   Corr        
b1.(Intercept) 0.4414536 b1.(I) b2    
            b2 0.3572042 -0.606       
            b3 0.3117688  0.212  0.238
Stratum: normal 
                  StdDev   Corr        
b1.(Intercept) 0.3941497 b1.(I) b2    
            b2 0.3102877 -0.532       
            b3 0.1961919  0.067  0.113
 Within-group standard deviation: 0.4178081

Variance function:
 Structure: Different standard deviations per stratum
 Formula:  ~ 1 | group 
 Parameter estimates:
 depression    normal 
          1 0.9706348
Fixed effects: list(b1 ~ group, d1 ~ 1) 
                   Value Std.Error  DF   t-value p-value 
b1.(Intercept)  2.913696 0.1101288 222  26.45717  <.0001
      b1.group -0.509698 0.1732330 222  -2.94227  0.0036
            d1  0.340090 0.1283891 222   2.64890  0.0087
 Correlation: 
        b1.(I)  b1.group 
b1.group -0.636       
     d1   0.000 -0.314

GML estimate(s) of smoothing parameter(s): 0.0002518025 
Equivalent Degrees of Freedom (DF):  11.69337 

Converged after 6 iterations

The fits are shown in Figures and . The right panel of Figure shows the estimated common function and its 95% Bayesian confidence intervals. Data from two groups are pooled to estimate the common function which has narrower confidence intervals. We conclude that the depressed subjects have their mean cortisol level elevated and less profound circadian rhythm than normal subjects.