In-Class exercise7： Reproduce the results of the loudness example by completing the R markdown file

1 Introduction

Ten subjects were played tones at each of 5 loudnesses, in random order. Subjects were asked to draw a line on paper whose length matched the loudness of the tone. Each subject repeated each loudness 3 times, for a total of 30 trials per subject. Here is the mean of the 3 log-lengths for each loudness, the sd of the three log-lengths, and the number of replications, which is always 3.

A data frame with 50 observations on the following 5 variables. - subject - a factor with unique values for each subject - loudness - either 50, 60, 70, 80 or 90 db. Decibels are a logrithmic scale - y - a numeric vector giving the mean of the log-lengths of three lines drawn. - sd - a numeric vector, giving the sd of the three log lengths - n - a numeric vector, equal to the constant value 3

2 Data management

#install.packages("alr4")
library(alr4)
data(Stevens, package="alr4")
library(tidyverse)
library(lme4)
glimpse(Stevens)
Rows: 50
Columns: 5
$ subject  <fct> A, A, A, A, A, B, B, B, B, B, C, C, C, C, C, D, D, D, D, D...
$ loudness <dbl> 50, 60, 70, 80, 90, 50, 60, 70, 80, 90, 50, 60, 70, 80, 90...
$ y        <dbl> 0.379, 1.727, 3.016, 3.924, 4.413, 0.260, 0.833, 2.009, 2....
$ sd       <dbl> 0.507, 0.904, 0.553, 0.363, 0.092, 0.077, 0.287, 0.396, 0....
$ n        <dbl> 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3...

3 OLS regression lines over 10 subjects

coef(m0 <- nlme::lmList(y ~ loudness | subject, data=Stevens))
  (Intercept) loudness
A     -4.4937  0.10265
B     -4.5853  0.09355
C     -2.2936  0.06950
D     -2.9194  0.07506
E     -5.3851  0.10391
F     -2.4212  0.07740
G     -2.8109  0.06497
H     -3.4207  0.08223
I     -1.9487  0.06561
J     -1.7591  0.05525
ggplot(Stevens, 
aes(loudness, 
y, 
group=subject))+
geom_point(alpha=.5) + 
stat_smooth(method="lm", 
formula=y ~ x, 
se=F,
col='black',
lwd=rel(.5)) +
labs(x="loudness", 
y="Mean Line length") +
theme_minimal()

4 Visualization

ggplot(data=Stevens, aes(x=loudness, y=y)) +
  geom_point(alpha=0.5) +
  stat_smooth(method="lm", formula=y ~ x, se=F) +
  facet_grid(. ~ subject) +
  labs(x="Loudness (in db)",
       y="Mean Line length (in log)") +
  theme_minimal()

5 A random intercepts model

library(lme4)
summary(m1 <- lmer(y ~ loudness + (1 | subject), data=Stevens))
Linear mixed model fit by REML ['lmerMod']
Formula: y ~ loudness + (1 | subject)
   Data: Stevens

REML criterion at convergence: 57.1

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.8081 -0.5060  0.0184  0.5530  1.7721 

Random effects:
 Groups   Name        Variance Std.Dev.
 subject  (Intercept) 0.1428   0.378   
 Residual             0.0986   0.314   
Number of obs: 50, groups:  subject, 10

Fixed effects:
            Estimate Std. Error t value
(Intercept) -3.20377    0.25405   -12.6
loudness     0.07901    0.00314    25.2

Correlation of Fixed Effects:
         (Intr)
loudness -0.865

6 Random effects

library(lattice)
qqmath(ranef(m1))
$subject

7 Residual plot

plot(m1, resid(., scaled=TRUE) ~ fitted(.) | subject, 
     xlab="Fitted values", ylab= "Standardized residuals",
     abline=0, lty=3)

8 Normality QQ-plot

qqmath(m1, grid=TRUE)

##The end