1030_random_coef

#peijun

Refer to the questions here.

Q1 PIAT scores

dta <- read.table("http://140.116.183.121/~sheu/lmm/Data/reading_piat.txt", header = T)

dta <- dta %>%
        mutate(cAge = Age - 6.5,
               ID = factor(ID))

summary(m0 <- lmer(PIAT ~ cAge + (cAge | ID), data = dta))

Linear mixed model fit by REML ['lmerMod']
Formula: PIAT ~ cAge + (cAge | ID)
   Data: dta

REML criterion at convergence: 1804.1

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.0042 -0.4893 -0.1383  0.4067  3.6892 

Random effects:
 Groups   Name        Variance Std.Dev. Corr
 ID       (Intercept)  5.507   2.347        
          cAge         3.377   1.838    0.53
 Residual             27.400   5.235        
Number of obs: 267, groups:  ID, 89

Fixed effects:
            Estimate Std. Error t value
(Intercept)  21.0621     0.5630   37.41
cAge          4.5399     0.2622   17.32

Correlation of Fixed Effects:
     (Intr)
cAge -0.288

Q2 Visual search

dta <- read.table("E:/201709/multilevel_analysis/1030/visual_search.csv", header=T, sep=",")

#plot
ggplot(dta, aes(Set.Size, RT, group = Participant))+
  geom_point()+ geom_line() +
  stat_smooth(aes(group = Dx),color = "gray", method = "lm") +
  facet_grid(. ~ Dx)+
  labs(x = "Set Size", y = "Response Time (ms)") +
  theme(legend.position = c(.9, .2)) +
  theme_bw()

summary(m0 <- lmer(RT ~ Set.Size + Dx + Set.Size : Dx + (Set.Size| Participant), data = dta))

Linear mixed model fit by REML ['lmerMod']
Formula: RT ~ Set.Size + Dx + Set.Size:Dx + (Set.Size | Participant)
   Data: dta

REML criterion at convergence: 2186.1

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.7578 -0.3130 -0.0750  0.3117  6.1372 

Random effects:
 Groups      Name        Variance Std.Dev. Corr
 Participant (Intercept) 657552.9 810.90       
             Set.Size       407.4  20.18   1.00
 Residual                772433.1 878.88       
Number of obs: 132, groups:  Participant, 33

Fixed effects:
                   Estimate Std. Error t value
(Intercept)         2078.75     270.98   7.671
Set.Size              73.49      11.40   6.446
DxControl          -1106.05     366.90  -3.015
Set.Size:DxControl   -21.74      15.44  -1.408

Correlation of Fixed Effects:
            (Intr) Set.Sz DxCntr
Set.Size    -0.071              
DxControl   -0.739  0.053       
St.Sz:DxCnt  0.053 -0.739 -0.071

#plot interval  ????? mean not on the line
library(ggplot2)
ggplot(dta, aes(Set.Size, RT, shape = Dx, linetype = Dx)) +
  stat_summary(fun.y = mean, geom = "line") +
  stat_summary(fun.y = mean, geom = "point") +
  stat_summary(fun.data = mean_se, geom = "errorbar",
               linetype = "solid", width = 0.5) +
  scale_shape_manual(values = c(1, 2)) +
  labs(x = "Set Size", y = "Response Time (ms)", linetype = "DX", shape = "Dx") +
  theme(legend.justification = c(0, 1), legend.position = c(0, 1),
        legend.background = element_rect(fill = "white", color = "black"))

#plot residual
plot(m0, resid(., type = "pearson") ~ fitted(.) | Dx,
     layout = c(3,1), aspect = 2, abline = 0, pch = 20, cex = .8,
     xlab = "Ftted values", ylab = "Pearson residuals")

Q3

Q4 Alcohol use

dta <- read.table("http://140.116.183.121/~sheu/lmm/Data/alcohol_use.txt", header=T, sep="")

dta <- dta %>%
        mutate(sid = factor(sid),
               coa = factor(coa),
               sex = factor(sex))

summary(m0 <- lmer(alcuse ~ coa + peer + age14 + peer:age14 + (1| sid), data = dta))

Linear mixed model fit by REML ['lmerMod']
Formula: alcuse ~ coa + peer + age14 + peer:age14 + (1 | sid)
   Data: dta

REML criterion at convergence: 621.6

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.30588 -0.66394 -0.05233  0.55906  2.60999 

Random effects:
 Groups   Name        Variance Std.Dev.
 sid      (Intercept) 0.3377   0.5811  
 Residual             0.4823   0.6945  
Number of obs: 246, groups:  sid, 82

Fixed effects:
            Estimate Std. Error t value
(Intercept) -0.31177    0.17225  -1.810
coa1         0.56513    0.15878   3.559
peer         0.69586    0.13219   5.264
age14        0.42469    0.09346   4.544
peer:age14  -0.15138    0.07481  -2.024

Correlation of Fixed Effects:
           (Intr) coa1   peer   age14 
coa1       -0.312                     
peer       -0.725 -0.134              
age14      -0.543  0.000  0.461       
peer:age14  0.442  0.000 -0.566 -0.814

summary(m1 <- lmer(alcuse ~ coa + peer + age14 + peer:age14 + (age14| sid), data = dta))

Linear mixed model fit by REML ['lmerMod']
Formula: alcuse ~ coa + peer + age14 + peer:age14 + (age14 | sid)
   Data: dta

REML criterion at convergence: 603.7

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.59995 -0.40432 -0.07739  0.44372  2.27436 

Random effects:
 Groups   Name        Variance Std.Dev. Corr 
 sid      (Intercept) 0.2595   0.5094        
          age14       0.1469   0.3832   -0.05
 Residual             0.3373   0.5808        
Number of obs: 246, groups:  sid, 82

Fixed effects:
            Estimate Std. Error t value
(Intercept) -0.31382    0.14871  -2.110
coa1         0.57120    0.14898   3.834
peer         0.69518    0.11322   6.140
age14        0.42469    0.10690   3.973
peer:age14  -0.15138    0.08556  -1.769

Correlation of Fixed Effects:
           (Intr) coa1   peer   age14 
coa1       -0.339                     
peer       -0.708 -0.146              
age14      -0.408  0.000  0.350       
peer:age14  0.332  0.000 -0.429 -0.814

Q5 Sleep deprivation

dta <- sleepstudy

dta <- dta %>%
        mutate(Days = factor(Days))

summary(m0 <- lmer(Reaction ~ Days + (1| Subject), data = dta))

Linear mixed model fit by REML ['lmerMod']
Formula: Reaction ~ Days + (1 | Subject)
   Data: dta

REML criterion at convergence: 1729.5

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.3473 -0.5293  0.0317  0.5328  4.2570 

Random effects:
 Groups   Name        Variance Std.Dev.
 Subject  (Intercept) 1375.5   37.09   
 Residual              987.6   31.43   
Number of obs: 180, groups:  Subject, 18

Fixed effects:
            Estimate Std. Error t value
(Intercept)  256.652     11.458  22.400
Days1          7.844     10.475   0.749
Days2          8.710     10.475   0.831
Days3         26.340     10.475   2.515
Days4         31.998     10.475   3.055
Days5         51.867     10.475   4.951
Days6         55.526     10.475   5.301
Days7         62.099     10.475   5.928
Days8         79.978     10.475   7.635
Days9         94.199     10.475   8.993

Correlation of Fixed Effects:
      (Intr) Days1  Days2  Days3  Days4  Days5  Days6  Days7  Days8 
Days1 -0.457                                                        
Days2 -0.457  0.500                                                 
Days3 -0.457  0.500  0.500                                          
Days4 -0.457  0.500  0.500  0.500                                   
Days5 -0.457  0.500  0.500  0.500  0.500                            
Days6 -0.457  0.500  0.500  0.500  0.500  0.500                     
Days7 -0.457  0.500  0.500  0.500  0.500  0.500  0.500              
Days8 -0.457  0.500  0.500  0.500  0.500  0.500  0.500  0.500       
Days9 -0.457  0.500  0.500  0.500  0.500  0.500  0.500  0.500  0.500

# plot
ggplot(dta, aes(Days, Reaction, group = Subject, color = Subject))+
  geom_point()+ geom_line() +
  labs(x = "Day", y = "Response Time (ms)") +
  theme_bw()

ggplot(dta, aes(Days, Reaction, color = Subject))+
  geom_point() +
  stat_smooth(aes(group = Subject), method = "lm", se = F) +
  labs(x = "Day", y = "Response Time (ms)") +
  theme_bw()

Q6 Alzheimer patients and cognitive recall

dta <- read.table("http://140.116.183.121/~sheu/lmm/Data/lecithin.txt", header = T, sep = "")

summary(m0 <- lmer(Recall ~ Visit + Group + (Group | Patient), data = dta))

Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
$checkConv, : unable to evaluate scaled gradient

Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
$checkConv, : Model failed to converge: degenerate Hessian with 1 negative
eigenvalues

Linear mixed model fit by REML ['lmerMod']
Formula: Recall ~ Visit + Group + (Group | Patient)
   Data: dta

REML criterion at convergence: 1269.6

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.73663 -0.56359 -0.08844  0.56694  3.11373 

Random effects:
 Groups   Name        Variance Std.Dev. Corr 
 Patient  (Intercept) 20.358   4.512         
          GroupP      14.602   3.821    -0.67
 Residual              8.275   2.877         
Number of obs: 235, groups:  Patient, 48

Fixed effects:
            Estimate Std. Error t value
(Intercept)   9.9455     1.0771   9.234
Visit         0.5060     0.1331   3.801
GroupP       -3.0195     1.2341  -2.447

Correlation of Fixed Effects:
       (Intr) Visit 
Visit  -0.371       
GroupP -0.751 -0.004
convergence code: 0
unable to evaluate scaled gradient
Model failed to converge: degenerate  Hessian with 1 negative eigenvalues

#plot
ggplot(dta, aes(Visit, Recall, group = Patient , color = Patient))+
  geom_point()+ geom_line() +
  facet_grid(. ~ Group) +
  labs(x = "Visit", y = "Recall score)") +
  theme_bw()

#plot
ggplot(dta, aes(Visit, Recall, shape = Group, linetype = Group)) +
  stat_summary(fun.y = mean, geom = "line") +
  stat_summary(fun.y = mean, geom = "point") +
  scale_shape_manual(values = c("1", "2")) +
  labs(x = "Visit", y = "Recall Score", linetype = "Group", shape = "Group")

Q7 Borckardt2014, mt and values

dta <- Borckardt2014

#plot
ggplot(dta, aes(mt, values)) +
  stat_smooth(method = "lm", se = F, linetype = "dashed")+
  geom_point(alpha = .5) +
  geom_vline(xintercept = 7, linetype = "dotted") +
  stat_smooth(aes(group = phase), method = "lm", se = F) +
  labs(x = "mt", y = "values")

#fit a linear piece-wise model vs simple linear model
summary(m1 <- lm(values ~ mt + phase, data = dta))

Call:
lm(formula = values ~ mt + phase, data = dta)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.93151 -0.91977  0.07436  1.06458  2.07241 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  7.722114   0.527293  14.645 1.92e-11 ***
mt          -0.001957   0.072773  -0.027    0.979    
phaseB      -4.765166   0.934791  -5.098 7.52e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.163 on 18 degrees of freedom
Multiple R-squared:  0.8144,    Adjusted R-squared:  0.7938 
F-statistic: 39.49 on 2 and 18 DF,  p-value: 2.612e-07

summary(m0 <- lm(values ~ mt, data = dta))

Call:
lm(formula = values ~ mt, data = dta)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.91342 -1.13420  0.00216  1.69697  2.52814 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  7.88095    0.80088   9.840 6.82e-09 ***
mt          -0.30519    0.06378  -4.785 0.000129 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.77 on 19 degrees of freedom
Multiple R-squared:  0.5465,    Adjusted R-squared:  0.5226 
F-statistic:  22.9 on 1 and 19 DF,  p-value: 0.0001288