Cortisol

1 Introduction

In a study on self-esteem and neuroendocrine response to challenge, cortisol levels were measured at 9 different occasions.

2 Data

# set path for data file
fLoc <- "https://quantdev.ssri.psu.edu/sites/qdev/files/TheCortisolData.csv"
# read in the .csv file using the url() function
dta <- read.csv(file=url(fLoc), header=TRUE)

head(dta)

  id cort_0 cort_1 cort_2 cort_3 cort_4 cort_5 cort_6 cort_7 cort_8
1  1    4.2    4.1    9.7   14.0   19.0   18.0   20.0   23.0   24.0
2  2    5.5    5.6   14.0   16.0   19.0   17.0   18.0   20.0   19.0
3  3    4.0    3.8    7.5   12.0   14.0   13.0    9.1    8.2    7.9
4  4    6.1    5.6   14.0   20.0   26.0   23.0   26.0   25.0   26.0
5  5    4.6    4.4    7.2   12.3   15.8   16.1   17.0   17.8   19.1
6  6    6.8    9.5   14.2   19.6   19.0   13.9   13.4   12.5   11.7

str(dta)

'data.frame':   34 obs. of  10 variables:
 $ id    : int  1 2 3 4 5 6 7 8 9 10 ...
 $ cort_0: num  4.2 5.5 4 6.1 4.6 6.8 7.4 9.2 3.9 9.3 ...
 $ cort_1: num  4.1 5.6 3.8 5.6 4.4 9.5 9.2 10 3.3 8.5 ...
 $ cort_2: num  9.7 14 7.5 14 7.2 14.2 14 16 9.4 11.5 ...
 $ cort_3: num  14 16 12 20 12.3 19.6 18 21 16 17 ...
 $ cort_4: num  19 19 14 26 15.8 19 19 24 18.1 21.6 ...
 $ cort_5: num  18 17 13 23 16.1 13.9 16 21 14.3 23.1 ...
 $ cort_6: num  20 18 9.1 26 17 13.4 16 19 13.7 23.7 ...
 $ cort_7: num  23 20 8.2 25 17.8 12.5 18 21 13.8 22.6 ...
 $ cort_8: num  24 19 7.9 26 19.1 11.7 18 25 13.9 24.7 ...

# make a copy of the original data in the long format
dtaL <- dta %>%                                        
 pivot_longer(cols=starts_with("cort"), names_to="time", values_to = "cort") %>%
 dplyr::mutate(zeit = substr(time, 6, 6),
               time = parse_number(time)/8.0, 
               id=as.factor(id)) %>%
 dplyr::arrange(id, time)                     # sort observations

head(dtaL)

# A tibble: 6 × 4
  id     time  cort zeit 
  <fct> <dbl> <dbl> <chr>
1 1     0       4.2 0    
2 1     0.125   4.1 1    
3 1     0.25    9.7 2    
4 1     0.375  14   3    
5 1     0.5    19   4    
6 1     0.625  18   5    

3 Summary statistics

dta %>%                
  furniture::table1("Baseline" = cort_0,
                    "Time 1" = cort_1,
                    "Time 2" = cort_2,
                    "Time 3" = cort_3,
                    "Time 4" = cort_4,
                    "Time 5" = cort_5,
                    "Time 6" = cort_6,
                    "Time 7" = cort_7,
                    "Time 8" = cort_8,
                    total = TRUE,
                    test = TRUE,
                    na.rm = FALSE,   # default: COMPLETE CASES ONLY
                    digits = 2,
                    align = "c",
                    output = "markdown",
                    caption = "Cortisol over time")

Cortisol over time

	Mean/Count (SD/%)
	n = 34
Baseline
	5.25 (2.45)
Time 1
	5.06 (2.59)
Time 2
	10.85 (3.11)
Time 3
	17.36 (3.69)
Time 4
	19.48 (3.72)
Time 5
	16.46 (4.18)
Time 6
	15.03 (4.57)
Time 7
	15.62 (4.81)
Time 8
	15.86 (5.61)

-在不同時間點下，平均Cortisol濃度不同。整體來說，隨著時間增加，平均Cortisol濃度增加。

3.1 Covariances and correlations

dta %>% 
  dplyr::select(starts_with("cort_")) %>%  
  var(use = "pairwise.complete.obs")  %>%  
  round(2)

       cort_0 cort_1 cort_2 cort_3 cort_4 cort_5 cort_6 cort_7 cort_8
cort_0   6.02   6.07   2.56   2.63   1.73   0.99   2.37   1.40   1.89
cort_1   6.07   6.73   3.34   3.00   1.65   0.88   2.14   1.37   1.70
cort_2   2.56   3.34   9.67   8.06   7.95   6.05   6.33   6.07   6.39
cort_3   2.63   3.00   8.06  13.58  10.24   5.98   4.20   2.91   4.28
cort_4   1.73   1.65   7.95  10.24  13.86  11.83  10.31   9.14  10.59
cort_5   0.99   0.88   6.05   5.98  11.83  17.48  16.86  17.21  18.84
cort_6   2.37   2.14   6.33   4.20  10.31  16.86  20.88  20.95  23.32
cort_7   1.40   1.37   6.07   2.91   9.14  17.21  20.95  23.13  25.19
 [ 達到了 getOption("max.print") -- 省略最後 1 列 ]]

-在不同時間點的預測平均Cortisol濃度的變異數及共變數矩陣

dta %>% 
  dplyr::select(starts_with("cort_")) %>% # just the outcome(s)
  cor(use = "pairwise.complete.obs") %>%   # correlation matrix
  corrplot::corrplot.mixed(upper = "ellipse")

-上半部為資料分布，下半部為相關係數，相近時間點的所測量的平均Cortisol濃度，相似性則越高，時間點相隔越久，相似性則越低。

psych::pairs.panels(dta[,-1])

-同上圖，為不同時間點測量平均Cortisol濃度的相關矩陣。左下角為觀察值的散布圖及迴歸線；對角線為根據觀察值繪製的直方圖；右上角為不同觀測時間的相關係數。

4 Visualization

# Use theme_set() to completely override the current theme.
old <- theme_set(theme_minimal())

ggplot(dtaL, 
       aes(x = time, y = cort, group = id)) +
  geom_point(size=rel(.8)) + 
  geom_line(alpha=.5, size=rel(.5)) +
  labs("Time (scaled) ", y="Cortisol") +
  scale_x_continuous(limits=c(0,1), breaks=seq(0, 1, by=.1)) 

-在不同觀測時間點下，每位受試者的平均Cortisol濃度。平均Cortisol濃度隨著時間增加，從Baseline到第1個觀測時間點增減幅度不大，但是從第2個觀測時間開始，平均Cortisol濃度顯著增加。到第5個觀測時間開始，平均Cortisol濃度增減幅度趨近於平緩，甚至有下降的趨勢。

5 Models

5.1 Quadratic trend

dtaL$timesq <- dtaL$time^2

m0 <- lme4::lmer(cort ~ 1 + time + timesq + (1 + time + timesq | id),
                    data = dtaL)
summary(m0, corr=FALSE)

Linear mixed model fit by REML ['lmerMod']
Formula: cort ~ 1 + time + timesq + (1 + time + timesq | id)
   Data: dtaL

REML criterion at convergence: 1649.9

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.2486 -0.6695  0.0349  0.5213  2.5559 

Random effects:
 Groups   Name        Variance Std.Dev. Corr       
 id       (Intercept)  1.31    1.14                
          time        28.45    5.33      0.63      
          timesq      26.58    5.16     -0.97 -0.44
 Residual              9.35    3.06                
Number of obs: 306, groups:  id, 34

Fixed effects:
            Estimate Std. Error t value
(Intercept)    3.601      0.469    7.67
time          41.051      2.188   18.76
timesq       -30.170      2.107  -14.32

-測量時間點預測平均Cortisol濃度的迴歸方程式為Y= -30.17 timesq + 41.051 time +3.601

texreg::knitreg(m0, 
                single.row = TRUE,
                stars = numeric(0),
                caption = "Model with Quadratic Trend",
                caption.above = TRUE,
                custom.note = "Model fit w/ REML")

Model with Quadratic Trend

	Model 1
(Intercept)	3.60 (0.47)
time	41.05 (2.19)
timesq	-30.17 (2.11)
AIC	1669.94
BIC	1707.18
Log Likelihood	-824.97
Num. obs.	306
Num. groups: id	34
Var: id (Intercept)	1.31
Var: id time	28.45
Var: id timesq	26.58
Cov: id (Intercept) time	3.86
Cov: id (Intercept) timesq	-5.73
Cov: id time timesq	-12.02
Var: Residual	9.35
Model fit w/ REML

-使用Quadratic Trend跑出模型的各項數值比較，包含截距、迴歸係數、變異數、共變數等。

5.1.1 Fitted values

#obtaining predicted scores for individuals and marginal prediction
dtaL %>%
  mutate(yhat = predict(m0),
         yhatb = predict(m0, re.form=NA)) %>%
 ggplot() +
  aes(x = time, 
      y = yhat, 
      group = id) +
  geom_point(alpha=.5, 
             size=rel(.5))+
  geom_line(size=rel(.5), 
            linetype="dotted") +
  geom_line(aes(x = time, 
                y = yhatb), 
            size=rel(.8))+
  labs(x = "Time (scaled)", 
       y = "Cortisol level")+
  scale_x_continuous(limits=c(0,1), breaks=seq(0, 1, by=.1)) 

-以上述迴歸方程式描繪的曲線。在不同觀測時間點下，每位受試者的平均Cortisol濃度。從Baseline開始平均Cortisol濃度隨時間增加而增加，但到第5個觀測時間之後，平均Cortisol濃度增加幅度趨近於平緩，甚至有下降的趨勢。

lme4::VarCorr(m0) %>% 
  print(., comp=c("Variance", "Std.Dev."), digits=3)

 Groups   Name        Variance Std.Dev. Corr       
 id       (Intercept)  1.31    1.14                
          time        28.45    5.33      0.63      
          timesq      26.58    5.16     -0.97 -0.44
 Residual              9.35    3.06                

-截距、time、timesq與殘差的變異數與相關係數比較。

5.2 Using orthogonal polynomials

m1 <- lme4::lmer(cort ~ poly(time, 2) + ( poly(time, 2) | id),
                    data = dtaL)
summary(m1, corr=FALSE)

Linear mixed model fit by REML ['lmerMod']
Formula: cort ~ poly(time, 2) + (poly(time, 2) | id)
   Data: dtaL

REML criterion at convergence: 1645.5

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.2485 -0.6695  0.0349  0.5213  2.5558 

Random effects:
 Groups   Name           Variance Std.Dev. Corr       
 id       (Intercept)      7.30    2.70               
          poly(time, 2)1 987.53   31.43     0.76      
          poly(time, 2)2  68.03    8.25    -0.17  0.51
 Residual                  9.35    3.06               
Number of obs: 306, groups:  id, 34

Fixed effects:
               Estimate Std. Error t value
(Intercept)      13.441      0.495   27.14
poly(time, 2)1   61.433      6.197    9.91
poly(time, 2)2  -48.240      3.369  -14.32

-使用orthogonal polynomials的迴歸式為Y= 61.433 poly(time, 2)1 - 48.240 poly(time, 2)2 +13.441

lme4::VarCorr(m1) %>% 
  print(., comp=c("Variance", "Std.Dev."), digits=3)

 Groups   Name           Variance Std.Dev. Corr       
 id       (Intercept)      7.30    2.70               
          poly(time, 2)1 987.53   31.43     0.76      
          poly(time, 2)2  68.03    8.25    -0.17  0.51
 Residual                  9.35    3.06               

-截距、poly(time, 2)1、poly(time, 2)2與殘差的變異數與相關係數比較

texreg::knitreg(list(m1,
                     m0),
                custom.model.names = c("Orthogonal polynomials",
                                       "Quadratic Trend"),
                single.row = TRUE,
                stars = numeric(0),
                caption = "Linear versus quadratic trend",
                caption.above = TRUE, 
                custom.note = "")

Linear versus quadratic trend

	Orthogonal polynomials	Quadratic Trend
(Intercept)	13.44 (0.50)	3.60 (0.47)
poly(time, 2)1	61.43 (6.20)
poly(time, 2)2	-48.24 (3.37)
time		41.05 (2.19)
timesq		-30.17 (2.11)
AIC	1665.54	1669.94
BIC	1702.78	1707.18
Log Likelihood	-822.77	-824.97
Num. obs.	306	306
Num. groups: id	34	34
Var: id (Intercept)	7.30	1.31
Var: id poly(time, 2)1	987.53
Var: id poly(time, 2)2	68.03
Cov: id (Intercept) poly(time, 2)1	64.95
Cov: id (Intercept) poly(time, 2)2	-3.72
Cov: id poly(time, 2)1 poly(time, 2)2	131.37
Var: Residual	9.35	9.35
Var: id time		28.45
Var: id timesq		26.58
Cov: id (Intercept) time		3.86
Cov: id (Intercept) timesq		-5.73
Cov: id time timesq		-12.02

-使用orthogonal polynomials與Quadratic Trend跑出模型的各項數值比較，包含截距、迴歸係數、變異數、共變數

anova(m1, m0)

Data: dtaL
Models:
m1: cort ~ poly(time, 2) + (poly(time, 2) | id)
m0: cort ~ 1 + time + timesq + (1 + time + timesq | id)
   npar  AIC  BIC logLik deviance Chisq Df Pr(>Chisq)
m1   10 1675 1712 -827.6     1655                    
m0   10 1675 1712 -827.6     1655     0  0           

-orthogonal polynomials與Quadratic Trend的迴歸式做變異數分析，兩者未達顯著差異

6 Piecewise models

dtaL <- dtaL %>%
  mutate(tk2 = ifelse(time < 0.13, 0, time),
         tk3 = ifelse(time < 0.55, 0, time))

m2 <- lme4::lmer(cort ~ 1 + poly(time, 2) + tk2 + tk3 + 
                 (1 | id) +
                 ( 0 + tk3 | id) + (0 + poly(time, 2) | id),
                 data = dtaL,
                 control = lmerControl(optimizer="optimx",    # get it to converge
                 calc.derivs = FALSE,
                 optCtrl = list(method = "nlminb",
                                starttests = FALSE,
                                kkt = FALSE)))
summary(m2, corr=FALSE)

Linear mixed model fit by REML ['lmerMod']
Formula: cort ~ 1 + poly(time, 2) + tk2 + tk3 + (1 | id) + (0 + tk3 |  
    id) + (0 + poly(time, 2) | id)
   Data: dtaL
Control: 
lmerControl(optimizer = "optimx", calc.derivs = FALSE, optCtrl = list(method = "nlminb",  
    starttests = FALSE, kkt = FALSE))

REML criterion at convergence: 1472

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-3.377 -0.496  0.005  0.484  2.716 

Random effects:
 Groups   Name           Variance Std.Dev. Corr 
 id       (Intercept)      5.68    2.38         
 id.1     tk3             21.92    4.68         
 id.2     poly(time, 2)1 184.70   13.59         
          poly(time, 2)2 160.03   12.65    -1.00
 Residual                  3.91    1.98         
Number of obs: 306, groups:  id, 34

Fixed effects:
               Estimate Std. Error t value
(Intercept)        1.65       1.67    0.98
poly(time, 2)1   -60.18      20.58   -2.92
poly(time, 2)2   -25.44       3.37   -7.54
tk2               31.82       3.21    9.91
tk3              -10.17       1.17   -8.67
optimizer (optimx) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')

-tk2為排除測量時間0.13以下預測平均Cortisol濃度的變項，從此時間點起每增加1單位平均Cortisol濃度增加31.82。tk3為排除測量時間0.55以下預測平均Cortisol濃度的變項，從此時間點起每增加1單位平均Cortisol濃度減少10.17

6.0.1 Fitted values

dtaL %>%
  mutate(yhat = predict(m2),
         yhatb = predict(m2, re.form=NA)) %>%
 ggplot() +
  aes(x = time, 
      y = yhat, 
      group = id) +
  geom_point(alpha=.5, 
             size=rel(.5))+
  geom_line(size=rel(.5), 
            linetype="dotted") +
  geom_line(aes(x = time, 
                y = yhatb), 
            size=rel(.8))+
  labs(x = "Time (scaled)", 
       y = "Cortisol level")+
  scale_x_continuous(limits=c(0,1), breaks=seq(0, 1, by=.1)) 

-在測量時間0.13至0.55期間，每增加1單位平均Cortisol濃度增加31.82。從測量時間0.55開始，平均Cortisol濃度出現下降趨勢，每增加1單位平均Cortisol濃度減少10.17。

texreg::knitreg(list(m1,
                     m2),
                custom.model.names = c("Polynomials",
                                       "Piecewise"),
                single.row = TRUE,
                stars = numeric(0),
                caption = "Polynomial versus Piecewise",
                caption.above = TRUE, 
                custom.note = "")

Polynomial versus Piecewise

	Polynomials	Piecewise
(Intercept)	13.44 (0.50)	1.65 (1.67)
poly(time, 2)1	61.43 (6.20)	-60.18 (20.58)
poly(time, 2)2	-48.24 (3.37)	-25.44 (3.37)
tk2		31.82 (3.21)
tk3		-10.17 (1.17)
AIC	1665.54	1493.99
BIC	1702.78	1534.94
Log Likelihood	-822.77	-735.99
Num. obs.	306	306
Num. groups: id	34	34
Var: id (Intercept)	7.30	5.68
Var: id poly(time, 2)1	987.53
Var: id poly(time, 2)2	68.03
Cov: id (Intercept) poly(time, 2)1	64.95
Cov: id (Intercept) poly(time, 2)2	-3.72
Cov: id poly(time, 2)1 poly(time, 2)2	131.37
Var: Residual	9.35	3.91
Var: id.1 tk3		21.92
Var: id.2 poly(time, 2)1		184.70
Var: id.2 poly(time, 2)2		160.03
Cov: id.2 poly(time, 2)1 poly(time, 2)2		-171.93

-使用Polynomials與Piecewise跑出模型的各項數值比較，包含迴歸係數、變異數、共變數等

7 References

Ram, N., & Grimm, K. (2007). Using simple and complex growth models to articulate developmental change: Matching theory to method. International Journal of Behavioral Development, 31, 303-316.

Seeman, T.E., Berkman, L.F. et al. (1995). Self-esteem and neuroendocrine response to challenge: MacArthur studies of successful aging Journal of Psychosomatic Research, 39(1), 69-84