Hw Exercise2

1 Hw exercise 2

##load and see data
# trial by trial reaction time data


# load packages
pacman::p_load(languageR, lattice, tidyverse, nlme, lmtest)

# load data
data(lexdec, package="languageR")

# query data info
?lexdec

## starting httpd help server ... done

# see first 6 lines of data
head(lexdec)

##   Subject       RT Trial Sex NativeLanguage Correct PrevType PrevCorrect
## 1      A1 6.340359    23   F        English correct     word     correct
## 2      A1 6.308098    27   F        English correct  nonword     correct
## 3      A1 6.349139    29   F        English correct  nonword     correct
## 4      A1 6.186209    30   F        English correct     word     correct
## 5      A1 6.025866    32   F        English correct  nonword     correct
## 6      A1 6.180017    33   F        English correct     word     correct
##         Word Frequency FamilySize SynsetCount Length  Class FreqSingular
## 1        owl  4.859812  1.3862944   0.6931472      3 animal           54
## 2       mole  4.605170  1.0986123   1.9459101      4 animal           69
## 3     cherry  4.997212  0.6931472   1.6094379      6  plant           83
## 4       pear  4.727388  0.0000000   1.0986123      4  plant           44
## 5        dog  7.667626  3.1354942   2.0794415      3 animal         1233
## 6 blackberry  4.060443  0.6931472   1.3862944     10  plant           26
##   FreqPlural DerivEntropy Complex      rInfl meanRT SubjFreq meanSize
## 1         74       0.7912 simplex -0.3101549 6.3582     3.12   3.4758
## 2         30       0.6968 simplex  0.8145080 6.4150     2.40   2.9999
## 3         49       0.4754 simplex  0.5187938 6.3426     3.88   1.6278
## 4         68       0.0000 simplex -0.4274440 6.3353     4.52   1.9908
## 5        828       1.2129 simplex  0.3977961 6.2956     6.04   4.6429
## 6         31       0.3492 complex -0.1698990 6.3959     3.28   1.5831
##   meanWeight      BNCw      BNCc       BNCd BNCcRatio BNCdRatio
## 1     3.1806 12.057065  0.000000   6.175602  0.000000  0.512198
## 2     2.6112  5.738806  4.062251   2.850278  0.707856  0.496667
## 3     1.2081  5.716520  3.249801  12.588727  0.568493  2.202166
## 4     1.6114  2.050370  1.462410   7.363218  0.713242  3.591166
## 5     4.5167 74.838494 50.859385 241.561040  0.679589  3.227765
## 6     1.1365  1.270338  0.162490   1.187616  0.127911  0.934882

# check data structure
str(lexdec)

## 'data.frame':    1659 obs. of  28 variables:
##  $ Subject       : Factor w/ 21 levels "A1","A2","A3",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ RT            : num  6.34 6.31 6.35 6.19 6.03 ...
##  $ Trial         : int  23 27 29 30 32 33 34 38 41 42 ...
##  $ Sex           : Factor w/ 2 levels "F","M": 1 1 1 1 1 1 1 1 1 1 ...
##  $ NativeLanguage: Factor w/ 2 levels "English","Other": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Correct       : Factor w/ 2 levels "correct","incorrect": 1 1 1 1 1 1 1 1 1 1 ...
##  $ PrevType      : Factor w/ 2 levels "nonword","word": 2 1 1 2 1 2 2 1 1 2 ...
##  $ PrevCorrect   : Factor w/ 2 levels "correct","incorrect": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Word          : Factor w/ 79 levels "almond","ant",..: 55 47 20 58 25 12 71 69 62 1 ...
##  $ Frequency     : num  4.86 4.61 5 4.73 7.67 ...
##  $ FamilySize    : num  1.386 1.099 0.693 0 3.135 ...
##  $ SynsetCount   : num  0.693 1.946 1.609 1.099 2.079 ...
##  $ Length        : int  3 4 6 4 3 10 10 8 6 6 ...
##  $ Class         : Factor w/ 2 levels "animal","plant": 1 1 2 2 1 2 2 1 2 2 ...
##  $ FreqSingular  : int  54 69 83 44 1233 26 50 63 11 24 ...
##  $ FreqPlural    : int  74 30 49 68 828 31 65 47 9 42 ...
##  $ DerivEntropy  : num  0.791 0.697 0.475 0 1.213 ...
##  $ Complex       : Factor w/ 2 levels "complex","simplex": 2 2 2 2 2 1 1 2 2 2 ...
##  $ rInfl         : num  -0.31 0.815 0.519 -0.427 0.398 ...
##  $ meanRT        : num  6.36 6.42 6.34 6.34 6.3 ...
##  $ SubjFreq      : num  3.12 2.4 3.88 4.52 6.04 3.28 5.04 2.8 3.12 3.72 ...
##  $ meanSize      : num  3.48 3 1.63 1.99 4.64 ...
##  $ meanWeight    : num  3.18 2.61 1.21 1.61 4.52 ...
##  $ BNCw          : num  12.06 5.74 5.72 2.05 74.84 ...
##  $ BNCc          : num  0 4.06 3.25 1.46 50.86 ...
##  $ BNCd          : num  6.18 2.85 12.59 7.36 241.56 ...
##  $ BNCcRatio     : num  0 0.708 0.568 0.713 0.68 ...
##  $ BNCdRatio     : num  0.512 0.497 2.202 3.591 3.228 ...

# select a subset of variables for current example
dta <- lexdec %>% select(Subject, Trial, RT, Frequency)

# set plot theme
ot <- theme_set(theme_bw())

# word frequency effect
ggplot(data = dta, aes(x = Frequency, y = RT)) +
 geom_point() +
 stat_smooth(method="lm") +
 labs(x = "Word Frequency", y = "Reaction Time (transformed)")

## `geom_smooth()` using formula 'y ~ x'

# trial effect by subject
ggplot(data = dta, aes(x = Trial, y = RT)) +
 geom_point(pch =1, size = rel(.5)) +
 geom_line() +
 facet_wrap(~ Subject) +
 labs(x = "Trial ID", y = "Reaction Time (transformed)")

# autocorrelation by individual
acf.fnc(dta, group = "Subject", time = "Trial", x = "RT", plot = TRUE)

# simple linear model ignoring trial dependence - least-squares
summary(m0 <- lm(RT ~ Frequency, data = dta))

## 
## Call:
## lm(formula = RT ~ Frequency, data = dta)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.55407 -0.16153 -0.03494  0.11699  1.08768 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  6.588778   0.022296 295.515   <2e-16 ***
## Frequency   -0.042872   0.004533  -9.459   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2353 on 1657 degrees of freedom
## Multiple R-squared:  0.05123,    Adjusted R-squared:  0.05066 
## F-statistic: 89.47 on 1 and 1657 DF,  p-value: < 2.2e-16

# MLE of the model above
summary(m0a <- gls(RT ~ Frequency, data = dta, method = "ML"))

## Generalized least squares fit by maximum likelihood
##   Model: RT ~ Frequency 
##   Data: dta 
##         AIC       BIC   logLik
##   -88.57767 -72.33576 47.28883
## 
## Coefficients:
##                 Value   Std.Error   t-value p-value
## (Intercept)  6.588778 0.022295932 295.51482       0
## Frequency   -0.042872 0.004532505  -9.45874       0
## 
##  Correlation: 
##           (Intr)
## Frequency -0.966
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -2.3560145 -0.6868771 -0.1485727  0.4974772  4.6250713 
## 
## Residual standard error: 0.2351709 
## Degrees of freedom: 1659 total; 1657 residual

# simple linear model with autocorrelation error using gls
summary(m1 <- update(m0a, correlation = corAR1(form = ~ Trial | Subject)))

## Generalized least squares fit by maximum likelihood
##   Model: RT ~ Frequency 
##   Data: dta 
##         AIC       BIC   logLik
##   -514.6024 -492.9465 261.3012
## 
## Correlation Structure: ARMA(1,0)
##  Formula: ~Trial | Subject 
##  Parameter estimate(s):
##      Phi1 
## 0.6376182 
## 
## Coefficients:
##                 Value   Std.Error  t-value p-value
## (Intercept)  6.589575 0.018923164 348.2280       0
## Frequency   -0.041220 0.003458059 -11.9201       0
## 
##  Correlation: 
##           (Intr)
## Frequency -0.87 
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -2.3238200 -0.7053022 -0.1811588  0.4533196  4.4989014 
## 
## Residual standard error: 0.2408258 
## Degrees of freedom: 1659 total; 1657 residual

# compare models
anova(m0a, m1)

##     Model df       AIC       BIC    logLik   Test  L.Ratio p-value
## m0a     1  3  -88.5777  -72.3358  47.28883                        
## m1      2  4 -514.6024 -492.9465 261.30120 1 vs 2 428.0247  <.0001

#
# collect residuals from the two models above
dta <- dta %>% mutate(res0 = rstandard(m0),
                      res1 = resid(m1, type = "normalized"))

# Durbin-Watson test of serial correlation
dwtest(res0 ~ Frequency, data = dta)

## 
##  Durbin-Watson test
## 
## data:  res0 ~ Frequency
## DW = 0.99757, p-value < 2.2e-16
## alternative hypothesis: true autocorrelation is greater than 0

#
dwtest(res1 ~ Frequency, data = dta)

## 
##  Durbin-Watson test
## 
## data:  res1 ~ Frequency
## DW = 2.2241, p-value = 1
## alternative hypothesis: true autocorrelation is greater than 0

# check acf for residuals
acf.fnc(dta, group = "Subject", time = "Trial", x = "res1", plot = TRUE)