Last day

library(tidyverse)

Registered S3 method overwritten by 'dplyr':
  method           from
  print.rowwise_df     
── Attaching packages ───────────────────────────────────── tidyverse 1.2.1 ──
✔ ggplot2 3.2.1     ✔ purrr   0.3.2
✔ tibble  2.1.3     ✔ dplyr   0.8.3
✔ tidyr   0.8.3     ✔ stringr 1.4.0
✔ readr   1.3.1     ✔ forcats 0.4.0
── Conflicts ──────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()

df <- read_csv("~/Dropbox/data/mbp-stats-bootcamp/example-trial-data.csv")

Parsed with column specification:
cols(
  ID = col_double(),
  Group = col_character(),
  Sex = col_character(),
  Age = col_double(),
  BV = col_double(),
  striatum = col_double(),
  paramedian = col_double(),
  visual = col_double()
)

df %>% head

df %>% group_by(Group, Age) %>% summarise(n=n()) %>% spread(Age, n)

library(ggplot2)
library(RColorBrewer)
colours <- c("purple", "light green", "green", "black")
p <- position_jitterdodge(jitter.width = 2, dodge.width = 4)
df %>% ggplot() + aes(x=Age, y=BV, colour=Group) + geom_jitter(position=p, alpha=0.3) + stat_summary(fun.data=mean_cl_normal, geom="pointrange", position=p) + stat_summary(fun.y=mean, geom="line", position=p) + scale_colour_manual(values=colours)

library(broom)
df %>% filter(Age == 50) %>% lm(BV ~ Group, .) %>% tidy

df %>% filter(Age == 80) %>% lm(BV ~ Group, .) %>% tidy

library(lmerTest)

Attaching package: ‘lmerTest’

The following object is masked from ‘package:lme4’:

    lmer

The following object is masked from ‘package:stats’:

    step

df <- df %>% mutate(normAge = Age - 50)
l1 <- lmer(BV ~ Group + normAge + (1|ID), data=df)
l2 <- lmer(BV ~ Group * normAge + (1|ID), data=df)
l3 <- lmer(BV ~ Group * poly(normAge, 2) + (1|ID), data=df)

df %>% mutate(l1 = fitted(l1), 
              l2 = fitted(l2),
              l3 = fitted(l3)) %>% 
  ggplot() + 
  aes(x=normAge, y=BV, colour=Group) + 
  geom_jitter(position=p, alpha=0.1) + 
  #stat_summary(fun.data=mean_cl_normal, geom="pointrange", position=p) +
  #stat_summary(fun.y=mean, geom="line", position=p) +
  #scale_colour_manual(values=colours) + 
  stat_summary(aes(y=l1), fun.y=mean, geom="line", linetype=2) +
  stat_summary(aes(y=l2), fun.y=mean, geom="line", linetype=3) + 
  stat_summary(aes(y=l3), fun.y=mean, geom="line", linetype=1)

  #geom_smooth(method="lm", formula=y ~ poly(x,2), se=F)

summary(l1)

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: BV ~ Group + normAge + (1 | ID)
   Data: df

REML criterion at convergence: 1345.4

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.6188 -0.4130 -0.0026  0.3902  3.3898 

Random effects:
 Groups   Name        Variance Std.Dev.
 ID       (Intercept) 268.4    16.38   
 Residual             219.2    14.81   
Number of obs: 155, groups:  ID, 71

Fixed effects:
              Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)  424.57952    3.49766  93.44125 121.389  < 2e-16 ***
GroupReact50 -37.61437    6.38648  66.22088  -5.890 1.42e-07 ***
GroupReact80 -60.58780    8.39921  64.06362  -7.214 7.80e-10 ***
GroupSilent  -62.74092    5.78900  68.09797 -10.838  < 2e-16 ***
normAge        0.08443    0.06608  97.30119   1.278    0.204    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) GrpR50 GrpR80 GrpSln
GroupRect50 -0.466                     
GroupRect80 -0.300  0.183              
GroupSilent -0.494  0.271  0.207       
normAge     -0.423  0.039 -0.100 -0.004

summary(l2)

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: BV ~ Group * normAge + (1 | ID)
   Data: df

REML criterion at convergence: 1341.6

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.7596 -0.3587 -0.0264  0.3438  3.7657 

Random effects:
 Groups   Name        Variance Std.Dev.
 ID       (Intercept) 254.5    15.95   
 Residual             213.1    14.60   
Number of obs: 155, groups:  ID, 71

Fixed effects:
                      Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)          422.56638    3.68717 112.70622 114.604  < 2e-16 ***
GroupReact50         -39.10711    7.15606 104.44823  -5.465 3.17e-07 ***
GroupReact80         -53.53417   10.80727 133.24074  -4.954 2.17e-06 ***
GroupSilent          -52.53611    6.94224 117.89998  -7.568 9.22e-12 ***
normAge                0.17416    0.08935  95.06327   1.949   0.0542 .  
GroupReact50:normAge   0.09853    0.17903  87.54688   0.550   0.5835    
GroupReact80:normAge  -0.23350    0.21340  87.32792  -1.094   0.2769    
GroupSilent:normAge   -0.44994    0.17759 103.14429  -2.534   0.0128 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) GrpR50 GrpR80 GrpSln normAg GR50:A GR80:A
GroupRect50 -0.515                                          
GroupRect80 -0.341  0.176                                   
GroupSilent -0.531  0.274  0.181                            
normAge     -0.543  0.280  0.185  0.288                     
GrpRct50:nA  0.271 -0.489 -0.092 -0.144 -0.499              
GrpRct80:nA  0.227 -0.117 -0.648 -0.121 -0.419  0.209       
GrpSlnt:nrA  0.273 -0.141 -0.093 -0.580 -0.503  0.251  0.211

df %>% group_by(Age, Group) %>% summarise(m=mean(BV)) %>% spread(Age, m)

summary(l3)

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: BV ~ Group * poly(normAge, 2) + (1 | ID)
   Data: df

REML criterion at convergence: 1256.2

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.5917 -0.3364  0.0111  0.4035  3.7016 

Random effects:
 Groups   Name        Variance Std.Dev.
 ID       (Intercept) 252.6    15.89   
 Residual             210.3    14.50   
Number of obs: 155, groups:  ID, 71

Fixed effects:
                               Estimate Std. Error       df t value Pr(>|t|)
(Intercept)                     426.680      3.084   66.596 138.363  < 2e-16
GroupReact50                    -37.097      6.276   67.821  -5.911 1.22e-07
GroupReact80                    -58.773      8.441   71.450  -6.963 1.34e-09
GroupSilent                     -63.256      5.640   67.945 -11.216  < 2e-16
poly(normAge, 2)1                48.721     22.383   92.390   2.177  0.03205
poly(normAge, 2)2                42.339     21.382   83.504   1.980  0.05098
GroupReact50:poly(normAge, 2)1   13.980     46.744   87.034   0.299  0.76560
GroupReact80:poly(normAge, 2)1  -90.858     58.826   85.463  -1.545  0.12616
GroupSilent:poly(normAge, 2)1  -122.696     44.925  102.038  -2.731  0.00744
GroupReact50:poly(normAge, 2)2  -58.775     45.979   83.861  -1.278  0.20468
GroupReact80:poly(normAge, 2)2    8.207     51.350   79.754   0.160  0.87343
GroupSilent:poly(normAge, 2)2   -63.939     39.402   82.646  -1.623  0.10845
                                  
(Intercept)                    ***
GroupReact50                   ***
GroupReact80                   ***
GroupSilent                    ***
poly(normAge, 2)1              *  
poly(normAge, 2)2              .  
GroupReact50:poly(normAge, 2)1    
GroupReact80:poly(normAge, 2)1    
GroupSilent:poly(normAge, 2)1  ** 
GroupReact50:poly(normAge, 2)2    
GroupReact80:poly(normAge, 2)2    
GroupSilent:poly(normAge, 2)2     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) GrpR50 GrpR80 GrpSln p(A,2)1 p(A,2)2 GR50:(A,2)1
GroupRect50 -0.491                                                 
GroupRect80 -0.365  0.180                                          
GroupSilent -0.547  0.269  0.200                                   
ply(nrA,2)1  0.026 -0.013 -0.010 -0.014                            
ply(nrA,2)2  0.009 -0.005 -0.003 -0.005  0.116                     
GR50:(A,2)1 -0.013  0.128  0.005  0.007 -0.479  -0.055             
GR80:(A,2)1 -0.010  0.005 -0.228  0.006 -0.380  -0.044   0.182     
GrpS:(A,2)1 -0.013  0.006  0.005  0.028 -0.498  -0.058   0.239     
GR50:(A,2)2 -0.004  0.082  0.002  0.002 -0.054  -0.465   0.290     
GR80:(A,2)2 -0.004  0.002  0.027  0.002 -0.048  -0.416   0.023     
GrpS:(A,2)2 -0.005  0.002  0.002  0.045 -0.063  -0.543   0.030     
            GR80:(A,2)1 GS:(A,2)1 GR50:(A,2)2 GR80:(A,2)2
GroupRect50                                              
GroupRect80                                              
GroupSilent                                              
ply(nrA,2)1                                              
ply(nrA,2)2                                              
GR50:(A,2)1                                              
GR80:(A,2)1                                              
GrpS:(A,2)1  0.190                                       
GR50:(A,2)2  0.020       0.027                           
GR80:(A,2)2 -0.371       0.024     0.194                 
GrpS:(A,2)2  0.024       0.174     0.252       0.226     

anova(l3)

Type III Analysis of Variance Table with Satterthwaite's method
                       Sum Sq Mean Sq NumDF  DenDF F value  Pr(>F)    
Group                   32536 10845.4     3 69.234 51.5747 < 2e-16 ***
poly(normAge, 2)          117    58.7     2 85.781  0.2793 0.75698    
Group:poly(normAge, 2)   2737   456.2     6 84.880  2.1693 0.05383 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

library(forcats)
df2 <- df %>% filter(Age < 85 & Group != "React80")

df2 %>% ggplot() + aes(x=normAge, y=BV, colour=Group) + geom_jitter(position=p, alpha=0.3) + stat_summary(fun.data=mean_cl_normal, geom="pointrange", position=p) + stat_summary(fun.y=mean, geom="line", position=p) 

l1 <- lmer(BV ~ Group + normAge + (1|ID), df2)
l2 <- lmer(BV ~ Group * normAge + (1|ID), df2)

df2 %>% mutate(l1 = fitted(l1), 
              l2 = fitted(l2)) %>% 
  ggplot() + 
  aes(x=normAge, y=BV, colour=Group) + 
  geom_jitter(position=p, alpha=0.1) + 
  #stat_summary(fun.data=mean_cl_normal, geom="pointrange", position=p) +
  #stat_summary(fun.y=mean, geom="line", position=p) +
  #scale_colour_manual(values=colours) + 
  stat_summary(aes(y=l1), fun.y=mean, geom="line", linetype=2) +
  stat_summary(aes(y=l2), fun.y=mean, geom="line", linetype=3)

  #geom_smooth(method="lm", formula=y ~ poly(x,2), se=F)

summary(l2)

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: BV ~ Group * normAge + (1 | ID)
   Data: df2

REML criterion at convergence: 1041.4

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.2750 -0.4019  0.0660  0.4140  3.4173 

Random effects:
 Groups   Name        Variance Std.Dev.
 ID       (Intercept) 221.0    14.87   
 Residual             281.7    16.78   
Number of obs: 118, groups:  ID, 65

Fixed effects:
                      Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)          423.43853    3.96414 104.41875 106.817  < 2e-16 ***
GroupReact50         -40.44715    7.59023  99.55533  -5.329 6.16e-07 ***
GroupSilent          -53.49536    7.45567 106.77084  -7.175 9.98e-11 ***
normAge                0.02775    0.14177  60.16236   0.196    0.845    
GroupReact50:normAge   0.29484    0.26883  56.02561   1.097    0.277    
GroupSilent:normAge   -0.24952    0.26195  60.48506  -0.953    0.345    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) GrpR50 GrpSln normAg GR50:A
GroupRect50 -0.522                            
GroupSilent -0.532  0.278                     
normAge     -0.581  0.303  0.309              
GrpRct50:nA  0.306 -0.544 -0.163 -0.527       
GrpSlnt:nrA  0.314 -0.164 -0.612 -0.541  0.285

library(broom.mixed)
GroupFunc_Age0 <- tidy(l2)$estimate[1]
GroupReact50_Age0 <- tidy(l2)$estimate[1] + tidy(l2)$estimate[2]
GroupFunc_Age30 <- tidy(l2)$estimate[1] + (tidy(l2)$estimate[4]*30)
GroupReact50_Age30 <-  tidy(l2)$estimate[1] + tidy(l2)$estimate[2] + ((tidy(l2)$estimate[4]+tidy(l2)$estimate[5])*30)
c(GroupFunc_Age0, GroupFunc_Age30, GroupReact50_Age0, GroupReact50_Age30)

[1] 423.4385 424.2711 382.9914 392.6692

simTrialData <- function(
  func_at_0 = 423,
  react50_at_0 = 382,
  silent_at_0 = 372,
  sd_at_0 = 20,
  age_func = 0,
  age_react=0.2,
  age_silent=-0.2,
  sd_age=0.1,
  n_per_group=30) {
  
  outcome = data.frame(
    group=rep(c("Func", "React50", "Silent"), each=n_per_group),
    baseline = 0,
    followup = 0
  )
  outcome$baseline[outcome$group == "Func"] <- rnorm(n_per_group,
                                                     func_at_0, 
                                                     sd_at_0)
  outcome$baseline[outcome$group == "React50"] <- rnorm(n_per_group,
                                                     react50_at_0, 
                                                     sd_at_0)
  outcome$baseline[outcome$group == "Silent"] <- rnorm(n_per_group,
                                                     silent_at_0, 
                                                     sd_at_0)
  outcome$followup[outcome$group == "Func"] <- rnorm(n_per_group,
                                                     age_func, 
                                                     sd_age)
  outcome$followup[outcome$group == "React50"] <- rnorm(n_per_group,
                                                        age_react, 
                                                        sd_age)
  outcome$followup[outcome$group == "Silent"] <- rnorm(n_per_group,
                                                       age_silent, 
                                                       sd_age)
  outcome$followup <- outcome$followup + outcome$baseline
  return(outcome)
  }

library(tidyverse)
trialData <- simTrialData(react50_at_0 = 423, silent_at_0 = 423, sd_at_0 = 10, sd_age=5, age_silent = -5, age_react = 5)
longTrialdata <- trialData %>% mutate(id=1:n()) %>% gather(time, volume, baseline, followup) %>% mutate(age=ifelse(time == "baseline", 0, 30))

library(ggplot2)
p <- position_jitterdodge(jitter.width = 2, dodge.width = 3)
ggplot(longTrialdata) + aes(x=age, y=volume, colour=group) + geom_jitter(position=p) + stat_summary(fun.y=mean, geom="line", position=p)

library(lme4)
l1 <- lmer(volume ~ age*group+(1|id), data=longTrialdata)
l3 <- lm(followup ~ baseline + group, trialData)
l4 <- trialData %>% mutate(diff=followup-baseline) %>% lm(diff ~ group, .)
l5 <- lm(followup ~ group, trialData)

summary(l1)

Linear mixed model fit by REML ['lmerMod']
Formula: volume ~ age * group + (1 | id)
   Data: longTrialdata

REML criterion at convergence: 1220.9

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.7170 -0.5434 -0.0397  0.5448  2.3142 

Random effects:
 Groups   Name        Variance Std.Dev.
 id       (Intercept) 92.02    9.592   
 Residual             13.50    3.674   
Number of obs: 180, groups:  id, 90

Fixed effects:
                  Estimate Std. Error t value
(Intercept)      427.27609    1.87541 227.831
age                0.02358    0.03162   0.746
groupReact50      -4.77310    2.65223  -1.800
groupSilent       -6.43592    2.65223  -2.427
age:groupReact50   0.18707    0.04472   4.183
age:groupSilent   -0.19625    0.04472  -4.388

Correlation of Fixed Effects:
            (Intr) age    grpR50 grpSln ag:R50
age         -0.253                            
groupRect50 -0.707  0.179                     
groupSilent -0.707  0.179  0.500              
ag:grpRct50  0.179 -0.707 -0.253 -0.126       
age:grpSlnt  0.179 -0.707 -0.126 -0.253  0.500

summary(l3)

Call:
lm(formula = followup ~ baseline + group, data = trialData)

Residuals:
    Min      1Q  Median      3Q     Max 
-9.8442 -3.7822 -0.5094  3.4988 15.6967 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  44.46034   23.18816   1.917 0.058510 .  
baseline      0.89760    0.05423  16.553  < 2e-16 ***
groupReact50  5.12328    1.34733   3.803 0.000267 ***
groupSilent  -6.54640    1.36751  -4.787 6.97e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 5.121 on 86 degrees of freedom
Multiple R-squared:  0.8224,    Adjusted R-squared:  0.8162 
F-statistic: 132.7 on 3 and 86 DF,  p-value: < 2.2e-16

summary(l4)

Call:
lm(formula = diff ~ group, data = .)

Residuals:
     Min       1Q   Median       3Q      Max 
-10.2575  -4.1079  -0.3534   3.1162  14.8111 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)    0.7075     0.9486   0.746    0.458    
groupReact50   5.6120     1.3416   4.183 6.84e-05 ***
groupSilent   -5.8874     1.3416  -4.388 3.19e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 5.196 on 87 degrees of freedom
Multiple R-squared:  0.4579,    Adjusted R-squared:  0.4454 
F-statistic: 36.74 on 2 and 87 DF,  p-value: 2.709e-12

summary(l5)

Call:
lm(formula = followup ~ group, data = trialData)

Residuals:
     Min       1Q   Median       3Q      Max 
-29.0862  -6.7429  -0.0378   7.8986  23.4601 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  427.9836     1.9019 225.030  < 2e-16 ***
groupReact50   0.8389     2.6897   0.312    0.756    
groupSilent  -12.3233     2.6897  -4.582 1.53e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 10.42 on 87 degrees of freedom
Multiple R-squared:  0.2566,    Adjusted R-squared:  0.2395 
F-statistic: 15.01 on 2 and 87 DF,  p-value: 2.507e-06

summary(l4)

Call:
lm(formula = diff ~ group, data = .)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.245974 -0.052020  0.006899  0.055652  0.293808 

Coefficients:
              Estimate Std. Error   t value Pr(>|t|)    
(Intercept)   -0.01544    0.01810    -0.853    0.396    
groupReact50  20.02013    0.02559   782.285   <2e-16 ***
groupSilent  -29.96352    0.02559 -1170.822   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.09912 on 87 degrees of freedom
Multiple R-squared:      1, Adjusted R-squared:      1 
F-statistic: 1.932e+06 on 2 and 87 DF,  p-value: < 2.2e-16

volume ~ group*age volume ~ group+age+age^{2+group:age+group:age}2