car effects emmeans magrittr multcomp
TRUE TRUE TRUE TRUE TRUE
summary(d <-read.delim( # summarize d, the result of loadingfile="_input/rts.csv", # this filestringsAsFactors=TRUE)) # change categorical variables into factors
CASE RT LENGTH LANGUAGE GROUP
Min. : 1 Min. : 271.0 Min. :3.000 english:4023 english :3961
1st Qu.:2150 1st Qu.: 505.0 1st Qu.:4.000 spanish:3977 heritage:4039
Median :4310 Median : 595.0 Median :5.000
Mean :4303 Mean : 661.5 Mean :5.198
3rd Qu.:6450 3rd Qu.: 732.0 3rd Qu.:6.000
Max. :8610 Max. :4130.0 Max. :9.000
NA's :248
CORRECT FREQPMW ZIPFFREQ SITE
correct :7749 Min. : 1.00 Min. :3.000 a:2403
incorrect: 251 1st Qu.: 17.00 1st Qu.:4.230 b:2815
Median : 42.00 Median :4.623 c:2782
Mean : 81.14 Mean :4.591
3rd Qu.: 101.00 3rd Qu.:5.004
Max. :1152.00 Max. :6.061
You can see what the variables/columns contain in the file _input/rts.r. Specifically, we are asking, does the (logged) reaction time to a word (in ms) vary as a function of
the Zipf frequency of that word (ZIPFFREQ);
the length of that word (LENGTH);
the language that word was presented in (LANGUAGE: spanish vs. english);
the speaker group that words was presented to (GROUP: english vs. heritage);
any pairwise interaction of these predictors? (Specifically, we are expecting that the effect of LENGTH will be stronger for when LANGUAGE is spanish.)
1.1.1 Exploration & preparation
As always, a decent amount of exploration is required to prepare the data for a hopefully unproblematic modeling process.
1.1.1.1 The response variable
We begin our exploration with the response variable RT. An initial histogram shows that we will probably want to transform the response given the long right tail. A log transformation seems to work well enough so that’s what we’ll go with:
par(mfrow=c(1, 2)) # make the plotting window have 1 rows & 2 columnshist(d$RT, main="")hist(log2(d$RT), main="")par(mfrow=c(1, 1)) # reset to defaultd$RT_log <-log2(d$RT)
Figure 1: The distribution of RT
1.1.1.2 Numeric predictors
Let’s now look at the numeric predictors. We begin with ZIPFFREQ, which looks ok:
par(mfrow=c(1, 2)) # make the plotting window have 1 rows & 2 columnshist(d$ZIPFFREQ) # hist(d$ZIPFFREQ, breaks="FD") # usually, also check the bin widthplot(main="RT (in ms, logged) as a function of Zipf frequency",pch=16, col="#00000020",xlab="Zipf frequency" , x=d$ZIPFFREQ,ylab="RT (in ms, logged)", y=d$RT_log); grid()abline(lm(d$RT_log ~ d$ZIPFFREQ), col="green", lwd=2)par(mfrow=c(1, 1)) # reset to default
Figure 2: The distribution of ZIPFFREQ & RT
Let’s now look at LENGTH, which looks ok:
par(mfrow=c(1, 2)) # make the plotting window have 1 rows & 2 columnshist(d$LENGTH) # hist(d$LENGTH, breaks="FD") # usually, also check the bin widthplot(main="RT (in ms, logged) as a function of Length",pch=16, col="#00000020",xlab="Length" , x=jitter(d$LENGTH),ylab="RT (in ms, logged)", y=d$RT_log); grid()abline(lm(d$RT_log ~ d$LENGTH), col="green", lwd=2)par(mfrow=c(1, 1)) # reset to default
Figure 3: The distribution of LENGTH & RT
1.1.1.3 Categorical predictors
Let’s now look at the categorical predictors. We begin with LANGUAGE, which looks ok:
qwe <-boxplot(main="RT (in ms, logged) as a function of stimulus language", d$RT_log ~ d$LANGUAGE, notch=TRUE,xlab="Stimulus language",ylab="RT (in ms, logged)"); grid()
Figure 4: The distribution of LANGUAGE & RT
Let’s now look at GROUP, which looks ok:
qwe <-boxplot(main="RT (in ms, logged) as a function of speaker group", d$RT_log ~ d$GROUP, notch=TRUE,xlab="Speaker group",ylab="RT (in ms, logged)"); grid()
Figure 5: The distribution of GROUP & RT
Do we have a good distribution for the interaction of the categorical variables? Yes:
table(d$LANGUAGE, d$GROUP)
english heritage
english 2005 2018
spanish 1956 2021
1.1.1.4 Deviance
First, we’ll quickly compute the deviance of the null model:
deviance(m_00 <-lm(RT_log ~1, data=d)) # 1581.951
[1] 1581.951
1.1.2 Model selection
The initial model involves all predictors and all their interactions:
summary(m_01 <-lm( # make/summarize the linear model m_01: RT_log ~1+# RT_log ~ an overall intercept (1)# & these predictors & their pairwise interactions (ZIPFFREQ+LANGUAGE+GROUP+LENGTH)^2,data=d, # those vars are in dna.action=na.exclude)) # skip cases with NA/missing data (good habit)
Call:
lm(formula = RT_log ~ 1 + (ZIPFFREQ + LANGUAGE + GROUP + LENGTH)^2,
data = d, na.action = na.exclude)
Residuals:
Min 1Q Median 3Q Max
-1.25135 -0.27954 -0.06143 0.20982 2.55535
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 9.007342 0.213938 42.103 < 2e-16 ***
ZIPFFREQ 0.006682 0.043985 0.152 0.87925
LANGUAGEspanish 0.308147 0.103832 2.968 0.00301 **
GROUPheritage 0.102615 0.098841 1.038 0.29922
LENGTH 0.106338 0.040699 2.613 0.00900 **
ZIPFFREQ:LANGUAGEspanish -0.010640 0.018610 -0.572 0.56751
ZIPFFREQ:GROUPheritage -0.001619 0.017674 -0.092 0.92702
ZIPFFREQ:LENGTH -0.019935 0.008436 -2.363 0.01814 *
LANGUAGEspanish:GROUPheritage -0.206348 0.019821 -10.410 < 2e-16 ***
LANGUAGEspanish:LENGTH 0.018155 0.010490 1.731 0.08355 .
GROUPheritage:LENGTH 0.001545 0.009436 0.164 0.86999
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.4234 on 7741 degrees of freedom
(248 observations deleted due to missingness)
Multiple R-squared: 0.1227, Adjusted R-squared: 0.1216
F-statistic: 108.3 on 10 and 7741 DF, p-value: < 2.2e-16
drop1(m_01, # drop each droppable predictor at a time from m_01 &test="F") %>%# test its significance w/ an F-test"["(order(.[,6]),) # order by p-values
anova(m_01, m_02, # compare m_01 to m_02test="F") # using an F-test
Analysis of Variance Table
Model 1: RT_log ~ 1 + (ZIPFFREQ + LANGUAGE + GROUP + LENGTH)^2
Model 2: RT_log ~ ZIPFFREQ + LANGUAGE + GROUP + LENGTH + ZIPFFREQ:LANGUAGE +
ZIPFFREQ:LENGTH + LANGUAGE:GROUP + LANGUAGE:LENGTH + GROUP:LENGTH
Res.Df RSS Df Sum of Sq F Pr(>F)
1 7741 1387.8
2 7742 1387.8 -1 -0.0015043 0.0084 0.927
drop1(m_02, # drop each droppable predictor at a time from m_02 &test="F") %>%# test its significance w/ an F-test"["(order(.[,6]),) # order by p-values
Let’s fit the 3rd model with the ns interaction GROUP:LENGTH deleted:
summary(m_03 <-update( # make m_03 an update of m_02, .~. # m_02, namely all of it (.~.), but then- GROUP:LENGTH)) # remove this interaction
Call:
lm(formula = RT_log ~ ZIPFFREQ + LANGUAGE + GROUP + LENGTH +
ZIPFFREQ:LANGUAGE + ZIPFFREQ:LENGTH + LANGUAGE:GROUP + LANGUAGE:LENGTH,
data = d, na.action = na.exclude)
Residuals:
Min 1Q Median 3Q Max
-1.25124 -0.27896 -0.06132 0.20991 2.55518
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 9.006924 0.207917 43.320 < 2e-16 ***
ZIPFFREQ 0.005935 0.043061 0.138 0.89038
LANGUAGEspanish 0.307637 0.103692 2.967 0.00302 **
GROUPheritage 0.102765 0.013427 7.654 2.19e-14 ***
LENGTH 0.107184 0.040362 2.656 0.00793 **
ZIPFFREQ:LANGUAGEspanish -0.010645 0.018600 -0.572 0.56711
ZIPFFREQ:LENGTH -0.019949 0.008434 -2.365 0.01804 *
LANGUAGEspanish:GROUPheritage -0.205486 0.019254 -10.673 < 2e-16 ***
LANGUAGEspanish:LENGTH 0.018175 0.010487 1.733 0.08313 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.4234 on 7743 degrees of freedom
(248 observations deleted due to missingness)
Multiple R-squared: 0.1227, Adjusted R-squared: 0.1218
F-statistic: 135.4 on 8 and 7743 DF, p-value: < 2.2e-16
anova(m_02, m_03, # compare m_02 to m_03test="F") # using an F-test
Analysis of Variance Table
Model 1: RT_log ~ ZIPFFREQ + LANGUAGE + GROUP + LENGTH + ZIPFFREQ:LANGUAGE +
ZIPFFREQ:LENGTH + LANGUAGE:GROUP + LANGUAGE:LENGTH + GROUP:LENGTH
Model 2: RT_log ~ ZIPFFREQ + LANGUAGE + GROUP + LENGTH + ZIPFFREQ:LANGUAGE +
ZIPFFREQ:LENGTH + LANGUAGE:GROUP + LANGUAGE:LENGTH
Res.Df RSS Df Sum of Sq F Pr(>F)
1 7742 1387.8
2 7743 1387.8 -1 -0.0053296 0.0297 0.8631
drop1(m_03, # drop each droppable predictor at a time from m_03 &test="F") %>%# test its significance w/ an F-test"["(order(.[,6]),) # order by p-values
Single term deletions
Model:
RT_log ~ ZIPFFREQ + LANGUAGE + GROUP + LENGTH + ZIPFFREQ:LANGUAGE +
ZIPFFREQ:LENGTH + LANGUAGE:GROUP + LANGUAGE:LENGTH
Df Sum of Sq RSS AIC F value Pr(>F)
LANGUAGE:GROUP 1 20.4159 1408.2 -13206 113.9055 < 2e-16 ***
ZIPFFREQ:LENGTH 1 1.0028 1388.8 -13314 5.5948 0.01804 *
LANGUAGE:LENGTH 1 0.5383 1388.4 -13316 3.0035 0.08313 .
ZIPFFREQ:LANGUAGE 1 0.0587 1387.9 -13319 0.3276 0.56711
<none> 1387.8 -13317
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Let’s fit the 4th model with the ns interaction ZIPFFREQ:LANGUAGE deleted:
summary(m_04 <-update( # make m_04 an update of m_03, .~. # m_03, namely all of it (.~.), but then- ZIPFFREQ:LANGUAGE)) # remove this interaction
Call:
lm(formula = RT_log ~ ZIPFFREQ + LANGUAGE + GROUP + LENGTH +
ZIPFFREQ:LENGTH + LANGUAGE:GROUP + LANGUAGE:LENGTH, data = d,
na.action = na.exclude)
Residuals:
Min 1Q Median 3Q Max
-1.24964 -0.27871 -0.06224 0.21043 2.55420
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 9.002136 0.207740 43.334 < 2e-16 ***
ZIPFFREQ 0.007370 0.042986 0.171 0.86387
LANGUAGEspanish 0.257516 0.055524 4.638 3.58e-06 ***
GROUPheritage 0.102783 0.013426 7.655 2.16e-14 ***
LENGTH 0.113334 0.038903 2.913 0.00359 **
ZIPFFREQ:LENGTH -0.021360 0.008065 -2.648 0.00811 **
LANGUAGEspanish:GROUPheritage -0.205373 0.019252 -10.668 < 2e-16 ***
LANGUAGEspanish:LENGTH 0.018385 0.010480 1.754 0.07943 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.4233 on 7744 degrees of freedom
(248 observations deleted due to missingness)
Multiple R-squared: 0.1227, Adjusted R-squared: 0.1219
F-statistic: 154.7 on 7 and 7744 DF, p-value: < 2.2e-16
anova(m_03, m_04, # compare m_03 to m_04test="F") # using an F-test
Analysis of Variance Table
Model 1: RT_log ~ ZIPFFREQ + LANGUAGE + GROUP + LENGTH + ZIPFFREQ:LANGUAGE +
ZIPFFREQ:LENGTH + LANGUAGE:GROUP + LANGUAGE:LENGTH
Model 2: RT_log ~ ZIPFFREQ + LANGUAGE + GROUP + LENGTH + ZIPFFREQ:LENGTH +
LANGUAGE:GROUP + LANGUAGE:LENGTH
Res.Df RSS Df Sum of Sq F Pr(>F)
1 7743 1387.8
2 7744 1387.9 -1 -0.058712 0.3276 0.5671
drop1(m_04, # drop each droppable predictor at a time from m_04 &test="F") %>%# test its significance w/ an F-test"["(order(.[,6]),) # order by p-values
Single term deletions
Model:
RT_log ~ ZIPFFREQ + LANGUAGE + GROUP + LENGTH + ZIPFFREQ:LENGTH +
LANGUAGE:GROUP + LANGUAGE:LENGTH
Df Sum of Sq RSS AIC F value Pr(>F)
LANGUAGE:GROUP 1 20.3956 1408.3 -13208 113.8018 < 2.2e-16 ***
ZIPFFREQ:LENGTH 1 1.2570 1389.1 -13314 7.0135 0.008106 **
LANGUAGE:LENGTH 1 0.5515 1388.4 -13318 3.0773 0.079429 .
<none> 1387.9 -13319
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
We’re done now – why? Because remember that we said “we are expecting that the effect of LENGTH will be stronger for when LANGUAGE is spanish”, meaning if the coefficient for LANGUAGEspanish:LENGTH is positive – which it is – we can take half its p-value and then it is in fact significant! We call the final model m_final (for easy copying and pasting across analyses) and do some housekeeping:
m_final <- m_04; rm(m_04); invisible(gc())
Normally, we would do some regression diagnostics now – I have to skip those now to discuss how to interpret the model, but we will return to this later! For now, we compute the confidence intervals …
ZIPFFREQ LENGTH GROUP LANGUAGE PRED PRED_LOW PRED_UPP
1 0.000000 0 english english 9.002136 8.594910 9.409361
2 1.000000 0 english english 9.009505 8.684218 9.334792
3 4.590539 0 english english 9.035967 8.950325 9.121610
4 0.000000 1 english english 9.115470 8.782580 9.448360
5 1.000000 1 english english 9.101480 8.835596 9.367365
6 4.590539 1 english english 9.051249 8.981904 9.120594
A maybe (!) nicer way to represent this might be this, but whether this is nicer or not is probably subjective:
with(nd, # use the data in nd andtapply( # apply PRED, # to the predicted reaction times# a grouping like this: the numerics in the rows & columns, the factors in the layers/sliceslist(round(ZIPFFREQ, 2), round(LENGTH, 2), "GROUP"=GROUP, "LANGUAGE"=LANGUAGE),# just list the values, but rounded to 2 decimals c)) %>%round(2)
, , GROUP = english, LANGUAGE = english
0 1 5.2
0 9.00 9.12 9.59
1 9.01 9.10 9.49
4.59 9.04 9.05 9.12
, , GROUP = heritage, LANGUAGE = english
0 1 5.2
0 9.10 9.22 9.69
1 9.11 9.20 9.59
4.59 9.14 9.15 9.22
, , GROUP = english, LANGUAGE = spanish
0 1 5.2
0 9.26 9.39 9.94
1 9.27 9.38 9.84
4.59 9.29 9.33 9.47
, , GROUP = heritage, LANGUAGE = spanish
0 1 5.2
0 9.16 9.29 9.84
1 9.16 9.27 9.74
4.59 9.19 9.22 9.37
1.1.3 Model evaluation & interpretation
Let’s look at the results.
1.1.3.1LANGUAGE:GROUP
Let’s visualize the nature of the effect of LANGUAGE:GROUP based on the predictions from effects:
lagr_d <-data.frame( # make lagr_d a data frame of lagr <-effect( # the effect"LANGUAGE:GROUP", # of LANGUAGE:GROUP m_final)) # in m_finalplot(lagr, # plot the effect of the interactionylim=range(d$RT_log, na.rm=TRUE), # w/ these y-axis limitsgrid=TRUE) # & w/ a gridplot(lagr, # plot the effect of the interactionylim=range(d$RT_log, na.rm=TRUE), # w/ these y-axis limitsx.var="GROUP", # w/ this predictor on the x-axisgrid=TRUE) # & w/ a grid
Figure 6: The interaction LANGUAGE:GROUP
Figure 7: The interaction LANGUAGE:GROUP
And here’s a version (of the first of the two effects plots) with everything: observed values, predicted values, and the confidence intervals of the predictions:
# split data & predictions up by levels of GROUP:d_split <-split(d, d$GROUP)lagr_d_split <-split(lagr_d, lagr_d$GROUP, drop=TRUE)par(mfrow=c(1, 2))stripchart(main="RT ~ LANGUAGE for GROUP='english'", # stripchart w/ main headingpch=16, col="#00000020", ylim=c(8, 12), # filled semi-transparent grey circles, y limits d_split$english$RT_log ~ d_split$english$LANGUAGE, # data: RT_log ~ LANG (for eng)xlab="Language", ylab="RT (in ms, logged)", # axis labelsvertical=TRUE, method="jitter"); grid() # vertical plot & jittered observations, add grid# add predicted means:points(seq(lagr_d_split$english$fit), lagr_d_split$english$fit, pch="X", col=c("red", "blue"))# add confidence intervals:for (i inseq(lagr_d_split$english$fit)) { # as many arrows as neededarrows(i, lagr_d_split$english$lower[i], # each from here i, lagr_d_split$english$upper[i], # to herecode=3, angle=90, col=c("red", "blue")[i]) # both tips w/ 90 degrees}stripchart(main="RT ~ LANGUAGE for GROUP='heritage'", # stripchart w/ main headingpch=16, col="#00000020", ylim=c(8, 12), # filled semi-transparent grey circles, y limits d_split$heritage$RT_log ~ d_split$heritage$LANGUAGE, # data: RT_log ~ LANG (for her)xlab="Language", ylab="RT (in ms, logged)", # axis labelsvertical=TRUE, method="jitter"); grid() # vertical plot & jittered observations, add grid# add predicted means:points(seq(lagr_d_split$heritage$fit), lagr_d_split$heritage$fit, pch="X", col=c("red", "blue"))# add confidence intervals:for (i inseq(lagr_d_split$heritage$fit)) { # as many arrows as neededarrows(i, lagr_d_split$heritage$lower[i], # each from here i, lagr_d_split$heritage$upper[i], # to herecode=3, angle=90, col=c("red", "blue")[i]) # both tips w/ 90 degrees}par(mfrow=c(1, 1))
Figure 8: The interaction LANGUAGE:GROUP
1.1.3.2LANGUAGE:LENGTH
Let’s visualize the nature of the effect of LANGUAGE:LENGTH based on the predictions from effects:
lale_d <-data.frame( # make lale_d a data frame of lale <-effect( # the effect"LANGUAGE:LENGTH", # of LANGUAGE:LENGTH m_final, # in m_finalxlevels=list( # for whenLENGTH=3:9))) # LENGTH is these valueslale_d <- lale_d[order(lale_d$LANGUAGE, lale_d$LENGTH),]; lale_d
plot(lale, # plot the effect of the interactionylim=range(d$RT_log, na.rm=TRUE), # w/ these y-axis limitsgrid=TRUE) # & w/ a grid
Figure 9: The interaction LANGUAGE:LENGTH
And here a version with everything: observed values, predicted values (red for LANGUAGE: english, blue for LANGUAGE: spanish), and the confidence intervals of the predictions (I’m building this up very stepwise, more so than I might do in actual research):
d_split <-split(d, d$LANGUAGE)lale_d_split <-split(lale_d, lale_d$LANGUAGE, drop=TRUE)plot(main="RT ~ LANGUAGE:LENGTH", type="n", # plot nothingxlab="Length" , x=d$LENGTH, # x-axisylab="RT (in ms, logged)", y=d$RT_log); grid() # y-axis & grid# plot observed data points for LANGUAGE being eng:points(jitter(d_split$eng$LENGTH), d_split$eng$RT_log, pch=16, col="#FF000010")lines(lale_d_split$eng$LENGTH, lale_d_split$eng$fit, col="red")# plot observed data points for LANGUAGE being spa:points(jitter(d_split$spa$LENGTH), d_split$spa$RT_log, pch=16, col="#0000FF10")lines(lale_d_split$spa$LENGTH, lale_d_split$spa$fit, col="blue")# plot confidence band for LANGUAGE being eng:polygon(col="#FF000030", border="red", # semi-transparent grey, for once w/ borderx=c( lale_d_split$eng$LENGTH, # x-coordinates left to right, thenrev(lale_d_split$eng$LENGTH)), # the same values back (in reverse order)y=c( lale_d_split$eng$lower, # y-coordinates lower boundaryrev(lale_d_split$eng$upper))) # y-coordinates upper boundary# plot confidence band for LANGUAGE being spa:polygon(col="#0000FF30", border="blue", # semi-transparent grey, for once w/ borderx=c( lale_d_split$spa$LENGTH, # x-coordinates left to right, thenrev(lale_d_split$spa$LENGTH)), # the same values back (in reverse order)y=c( lale_d_split$spa$lower, # y-coordinates lower boundaryrev(lale_d_split$spa$upper))) # y-coordinates upper boundary# plot legendlegend(title="Language", bty="n",x=6, xjust=0.5, y=11.5, yjust=0.5, ncol=2,legend=levels(d$LANGUAGE), fill=c("red", "blue"))
Figure 10: The interaction LANGUAGE:LENGTH
1.1.3.3ZIPFFREQ:LENGTH
Let’s visualize the nature of the effect of ZIPFFREQ:LENGTH based on the predictions from effects:
(zile_d <-data.frame( # make zile_d a data frame of zile <-effect( # the effect"ZIPFFREQ:LENGTH", # of ZIPFFREQ:LENGTH m_final, # in m_finalxlevels=list( # for whenZIPFFREQ=seq(3, 6, 0.25), # ZIPFFREQ is these valuesLENGTH=seq(3, 9, 0.5))))) # LENGTH is these values
plot(zile, # plot the effect of the interactionylim=range(d$RT_log, na.rm=TRUE), # w/ these y-axis limitsgrid=TRUE) # & w/ a grid
Figure 11: The interaction ZIPFFREQ:LENGTH
Not that great … There are two main ways we could try and plot this better. One would be with an actual 3-D plot. The simplest version of this might use rgl::plot3D like this, which would open a new window with a rotatable 3-D plot (but this wouldn’t be rendered into this HTML properly, which is why I am not executing it here):
While this is nice, it’s also hard to publish in print. Let’s do this differently with a numeric heatmap kind of plot that I often use. For that, the numeric predictions are first cut into 9 ordinal bins with lower/higher bin numbers reflecting lower/higher numeric predictions …
zile_d$fitcat <-as.numeric( # create a numeric column fitcat in zile_dcut(zile_d$fit, # by cutting the predicted values upbreaks=seq( # at the values ranging frommin(zile_d$fit), # the smallest obs/pred RT tomax(zile_d$fit), # the largest obs/pred RTlength.out=10), # into 9 groupsinclude.lowest=TRUE, # include the lowest, andlabels=1:9)) # use the numbers 1:9 as labels
… and then we plot them into a coordinate system:
plot(main="RT ~ ZIPFFREQ:LENGTH", type="n", # plot nothing, justxlab="Zipf frequency", xlim=range(zile_d$ZIPFFREQ), x=0, # get the coordinateylab="Length" , ylim=range(zile_d$LENGTH) , y=0) # system set uptext( # plot textx=zile_d$ZIPFFREQ, y=zile_d$LENGTH, # at these coordinates, namelylabels=zile_d$fitcat, # the 'categorical/ordinal' predictions# make the physical size of the numbers dependent on the valuescex=0.5+zile_d$fitcat/7,# make the darkness of the grey reflect the predicted RT_logcol=grey(level=zile_d$fitcat/10)) # see ?grey
Figure 12: The interaction ZIPFFREQ:LENGTH
However, there’s one thing that’s possibly problematic about this – what is that? The problem is that we are binning with cut on the basis of zile_d$fit (check the above code), but the range of those predictions is much smaller than the range of the actually observed reaction times in d$RT_log:
range(zile_d$fit)
[1] 9.120337 9.674938
range(d$RT_log)
[1] NA NA
That in turn means the binning steps based on zile_d$fit will be based on much smaller differences between predicted values (of only 40, check levels(cut(zile_d$fit, breaks=seq(min(zile_d$fit), max(zile_d$fit), length.out=10)))), meaning a difference of predictions of only 81 ms could lead to a difference in 3 bins – if we were to bin based on the actual value range, that would not happen, because there the binning steps are 429 ms wide (check levels(cut(d$RT_log, breaks=seq(min(d$RT_log), max(d$RT_log), length.out=10)))). Thus, the above plot exaggerates the effects somewhat.
Accordingly, we determine the range of the values to be plotted in a potentially more representative way, by using
as a minimum, the minimum of the predictions and the actual values combined;
as a maximum, the maximum of the predictions and the actual values combined:
# define the minimum of the range to be used for binningmin_of_obs_n_pred <-min(c( d$RT_log, zile_d$fit),na.rm=TRUE)# define the maximum of the range to be used for binningmax_of_obs_n_pred <-max(c( d$RT_log, zile_d$fit),na.rm=TRUE)
We can see that the initial plot indeed exaggerated the range: the upper plot showed prediction bins from 1 to 9, but we now see that the predicctions are actually only in a really small range of values.
zile_d$fitcat2 <-as.numeric( # create a numeric column fitcat2 in zile_dcut(zile_d$fit, # by cutting the predicted values upbreaks=seq( # at the values ranging from min_of_obs_n_pred, # the smallest obs/pred RT_log to <--------- max_of_obs_n_pred, # the largest obs/pred RT_log <---------length.out=10), # into 9 groupsinclude.lowest=TRUE, # include the lowest, andlabels=1:9)) # use the numbers 1:9 as labelstable(zile_d$fitcat, zile_d$fitcat2)
Thus, we plot this again (with some additional info in the form of the decile-based grid), which here doesn’t help much:
plot(main="RT ~ ZIPFFREQ:LENGTH", type="n", # plot nothing, justxlab="Zipf frequency", xlim=range(zile_d$ZIPFFREQ), x=0, # get the coordinateylab="Length" , ylim=range(zile_d$LENGTH) , y=0) # system set uptext( # plot textx=zile_d$ZIPFFREQ, y=zile_d$LENGTH, # at these coordinates, namelylabels=zile_d$fitcat2, # the 'categorical/ordinal' predictions# make the physical size of the numbers dependent on the valuescex=0.5+zile_d$fitcat2/7,# make the darkness of the grey reflect the predicted RT_logcol=grey(level=zile_d$fitcat2/15)) # see ?greyabline(h=quantile(d$LENGTH , probs=seq(0, 1, 0.1)), lty=3, col="#00000020")abline(v=quantile(d$ZIPFFREQ, probs=seq(0, 1, 0.1)), lty=3, col="#00000020")
Figure 13: The interaction ZIPFFREQ:LENGTH
1.1.4 Write-up
To determine whether the reaction time to a word (in ms) varies as a function of
the Zipf frequency of that word (ZIPFFREQ);
the length of that word (LENGTH);
the language that word was presented in (LANGUAGE: english vs. spanish);
the speaker group that words was presented to (GROUP: english vs. heritage);
any pairwise interaction of these predictors with the specific expectation that the effect of LENGTH will be stronger for when LANGUAGE is spanish
a linear model selection process was undertaken. In order to deal with the very long right tail of the response variable the response variable was logged; then a backwards stepwise model selection of all predictors and controls (based on significance testing) was applied. The final model resulting from the elimination of ns predictors contained the effects of ZIPFFREQ:LENGTH, LANGUAGE:GROUP, and LANGUAGE:LENGTH and was highly significant (F=154.7, df1=7, df2=7744, p<0.0001), but it explains only little of the response variable (mult. R2=0.123, adj. R2=0.122). The summary table of the model is shown here.
as hypothesized, the slowing-down effect of LENGTH is stronger for Spanish than for English words [show plot];
the interaction LANGUAGE:GROUP has the following effects [show at least one plot]:
reaction times are always higher when LANGUAGE is spanish than when LANGUAGE is english, but that difference is higher when GROUP is english than when GROUP is spanish;
when LANGUAGE is english, heritage speakers need more time than learners, but when LANGUAGE is spanish, heritage speakers need less time than learners;
all 4 means of LANGUAGE:GROUP are significantly different from each other (non-overlapping CIs);
the interaction LANGUAGE:LENGTH has the hypothesized effect that the effect of LENGTH is stronger – in fact, nearly 2.25 times as strong – when LANGUAGE is spanish than when LANGUAGE is english [show at least one plot];
the interaction of ZIPFFREQ:LENGTH is not strong: both variables behave as expected – LENGTH slows down, ZIPFFREQ speeds up – but
the effect of LENGTH is much stronger with rarer words than with frequent words;
the effect of ZIPFFREQ is much stronger with longer words than with shorter words.
1.2 Binary logistic regression modeling
To do a quick recap of linear modeling, we will look at a data set in _input/genitives.csv involving genitive choices as the response variable:
car effects emmeans magrittr multcomp
TRUE TRUE TRUE TRUE TRUE
source("_helpers/C.score.r"); source("_helpers/R2s.r"); # helper functionssummary(d <-read.delim( # summarize d, the result of loadingfile="_input/genitives.csv", # this filestringsAsFactors=TRUE)) # change categorical variables into factors
CASE GENITIVE SPEAKER MODALITY POR_LENGTH
Min. : 2 of:2720 nns:2666 spoken :1685 Min. : 1.00
1st Qu.:1006 s : 880 ns : 934 written:1915 1st Qu.: 8.00
Median :2018 Median : 11.00
Mean :2012 Mean : 14.58
3rd Qu.:3017 3rd Qu.: 17.00
Max. :4040 Max. :204.00
PUM_LENGTH POR_ANIMACY POR_FINAL_SIB POR_DEF
Min. : 2.00 animate : 920 absent :2721 definite :2349
1st Qu.: 6.00 collective: 607 present: 879 indefinite:1251
Median : 9.00 inanimate :1671
Mean : 10.35 locative : 243
3rd Qu.: 13.00 temporal : 159
Max. :109.00
You can see what the variables/columns contain in the file _input/genitives.r. Specifically, we want to see whether the choice of a genitive varies as a function of
whether the speaker is a native speaker or not (SPEAKER);
the length difference between possessor and possessum and whether native and non-native speakers react to it in the same way;
the animacy of the possessor and whether native and non-native speakers react to it in the same way (POR_ANIMACY);
whether the possessor ends in a sibilant or not and whether native and non-native speakers react to it in the same way (POR_FINAL_SIB);
whether the modality is speaking or writing and whether native and non-native speakers react to it in the same way (MODALITY);
the interaction of POR_FINAL_SIB with MODALITY.
1.2.1 Exploration & preparation
1.2.1.1 All variables
As always, a decent amount of exploration is required to prepare the data for a hopefully unproblematic modeling process. We will only focus on the four variables of interest (and even there will take some didactically motivated shortcuts). We begin our exploration with the creation of the length difference predictor:
Let’s look at the numeric predictor first by exploring its ‘components’:
par(mfrow=c(1, 2)) # make the plotting window have 1 rows & 2 columnshist(d$POR_LENGTH, main="")hist(d$PUM_LENGTH, main="")par(mfrow=c(1, 1)) # reset to default
Figure 14: The lengths of possessors/-ums
This doesn’t look all that great because it’s not likely that the long tails will magically cancel each other out, and they don’t:
par(mfrow=c(1, 2)) # make the plotting window have 1 rows & 2 columnshist(d$POR_LENGTH - d$PUM_LENGTH, main="")hist(d$POR_LENGTH - d$PUM_LENGTH, breaks="FD", main="")par(mfrow=c(1, 1)) # reset to default
Figure 15: The difference of the logged lengths
What if we compute not the difference of the lengths, but the difference of the logged lengths?
of s
.
absent nns 1462 556
ns 500 203
present nns 562 86
ns 196 35
d$MODALITY %>%ftable(d$SPEAKER, d$GENITIVE)
of s
.
spoken nns 914 347
ns 318 106
written nns 1110 295
ns 378 132
Not too bad. There are some low-frequency cells and for the numeric predictor even some empty cells, but that problem will not surface in the model because there we won’t use a factorized version of that predictor.
1.2.1.2 Deviance
First, we’ll quickly compute the deviance of the null model:
summary(m_01 <-glm( # make/summarize the gen. linear model m_01: GENITIVE ~1+# GENITIVE ~ an overall intercept (1)# these predictors & their interaction, as specified above: SPEAKER * (LEN_PORmPUM_LOG + POR_ANIMACY + POR_FINAL_SIB + MODALITY) + POR_FINAL_SIB:MODALITY,family=binomial, # resp = binarydata=d, # those vars are in dna.action=na.exclude)) # skip cases with NA/missing data
pchisq( # compute the area under the chi-squared curveq=m_01$null.deviance-m_01$deviance, # for this chi-squared value: null - res. dev.df=m_01$df.null-m_01$df.residual, # at this df: null - residual dflower.tail=FALSE) # only the right/upper tail/side
[1] 0
drop1(m_01, # drop each droppable predictor at a time from m_01 &test="Chisq") %>%# test its significance w/ a LRT"["(order(.[,4]),) # order by p-values
Let’s fit the 2nd model with the 2-way interaction SPEAKER:LEN_PORmPUM_LOG deleted:
summary(m_02 <-update( # make m_02 an update of m_01, .~. # m_01, namely all of it (.~.), but then- SPEAKER:LEN_PORmPUM_LOG)) # remove this interaction
drop1(m_02, # drop each droppable predictor at a time from m_02 &test="Chisq") %>%# test its significance w/ a LRT"["(order(.[,4]),) # order by p-values
drop1(m_03, # drop each droppable predictor at a time from m_03 &test="Chisq") %>%# test its significance w/ a LRT"["(order(.[,4]),) # order by p-values
drop1(m_04, # drop each droppable predictor at a time from m_04 &test="Chisq") %>%# test its significance w/ a LRT"["(order(.[,4]),) # order by p-values
1.2.3 Model evaluation & interpretation: numerical & visual
We first assess the predictive power of the model in various ways. For that, we first compute all numeric and categorical predictions (as well as PRED_PP_OBS, which is the predicted probability of the level of the response variable that should have been predicted) …
d$PRED_PP_2 <-predict( # make d$PRED_PP_2 the result of predicting m_final, # from m_finaltype="response") # predicted probabilities of "s"d$PRED_LO_2 <-predict( # make d$PRED_LO_2 the result of predicting m_final) # from m_final -- this time the log oddsd$PRED_CAT <-factor( # d$PRED_CAT is determined by ifelseifelse(d$PRED_PP_2>=0.5, # if the predicted prob. of "s" is >=0.5levels(d$GENITIVE)[2], # predict "s"levels(d$GENITIVE)[1]), # otherwise, predict "of"levels=levels(d$GENITIVE)) # make sure both levels are registeredd$PRED_PP_OBS <-# d$PRED_PP_OBS is determined by ifelseifelse(d$GENITIVE==levels(d$GENITIVE)[2], # if the obs gen = 2nd level of GENITIVE d$PRED_PP_2, # take its predicted prob.1-d$PRED_PP_2) # otherwise take 1 minus its predicted probabilityhead(d)
CASE GENITIVE SPEAKER MODALITY POR_LENGTH PUM_LENGTH POR_ANIMACY
1 2 of nns spoken 13 7 collective
2 3 of nns spoken 22 7 animate
3 4 of nns spoken 11 8 animate
4 5 of nns spoken 26 4 collective
5 6 s nns spoken 8 4 animate
6 7 s nns spoken 7 3 animate
POR_FINAL_SIB POR_DEF LEN_PORmPUM_LOG PRED_PP_2 PRED_LO_2 PRED_CAT
1 absent definite 0.8930848 0.20413276 -1.3606618015 of
2 absent definite 1.6520767 0.50009947 0.0003978715 s
3 present definite 0.4594316 0.47394893 -0.1042987321 of
4 absent definite 2.7004397 0.06303183 -2.6990094172 of
5 absent definite 1.0000000 0.61851765 0.4832611329 s
6 absent definite 1.2223924 0.57897800 0.3185793854 s
PRED_PP_OBS
1 0.7958672
2 0.4999005
3 0.5260511
4 0.9369682
5 0.6185176
6 0.5789780
… and then we evaluate the confusion matrix:
(c_m <-table( # confusion matrix: cross-tabulateOBS =d$GENITIVE, # observed choices in the rowsPREDS=d$PRED_CAT)) # predicted choices in the columns
PREDS
OBS of s
of 2485 235
s 323 557
c( # evaluate the confusion matrix"Class. acc."=mean(d$GENITIVE==d$PRED_CAT, na.rm=TRUE), # (tp+tn) / (tp+tn+fp+fn)"Prec. for s"=c_m["s","s"] /sum(c_m[ ,"s" ]),"Acc./rec. for s"=c_m["s","s"] /sum(c_m["s" , ]),"Prec. for of"=c_m["of","of"] /sum(c_m[ ,"of"]),"Acc./rec. for of"=c_m["of","of"] /sum(c_m["of", ]))
Class. acc. Prec. for s Acc./rec. for s Prec. for of
0.8450000 0.7032828 0.6329545 0.8849715
Acc./rec. for of
0.9136029
Let’s now look at the effects/interactions one by one.
1.2.3.1 Main effect 1: LEN_PORmPUM_LOG
Here’s the effects object for the first interaction, LEN_PORmPUM_LOG:
ld_d <-data.frame(ld <-effect("LEN_PORmPUM_LOG", m_final,xlevels=list(LEN_PORmPUM_LOG=seq(-5, 5, 1))))# we add the predicted log odds as well:ld_d$fit.lo <-qlogis(ld_d$fit)# for later, we add the proportions of values in the data:ld_d$PROP <- d$LEN_PORmPUM_LOG %>%cut(breaks=seq(-5.5, 5.5, 1), include.lowest=TRUE) %>% table %>% prop.table %>% as.numericld_d %>%round(4)
op <-par(mar=c(4, 4, 4, 4)) # make the plotting margins big enoughplot(pch=16, col="#00000060", type="b", yaxp=c(0.1, 0.9, 4),cex=0.5+ (10* ld_d$PROP),xlab="Difference of logged lengths (or-um)" , x=ld_d$LEN_PORmPUM_LOG,ylab="Pred. prob. of s-genitives" , ylim=c(0, 1), y=ld_d$fit)abline(h=seq(0.1, 0.9, 0.2), lty=3, col="#00000030")abline(v=quantile(d$LEN_PORmPUM_LOG, prob=seq(0, 1, 0.1)), lty=3, col="#00000030")polygon(border=NA, col="#00000020",c(ld_d$LEN_PORmPUM_LOG, rev(ld_d$LEN_PORmPUM_LOG)),c(ld_d$lower , rev(ld_d$upper)))axis(4, at=seq(0.1, 0.9, 0.2), labels=round(qlogis(seq(0.1, 0.9, 0.2)), 2))mtext("Pred. log odds", 4, 2.5)par(op) # reset to defaults
Figure 19: The effect of length differences
Interpretation: a very clear short-before-long effect: when the possessor is longer than the possessum (on the right), s-genitives, which would have the long possessor in the front, are very unlikely and the other way round.
1.2.3.2 Main effect 2: POR_FINAL_SIB
Here’s the effects object for the second main effect, POR_FINAL_SIB:
fs_d <-data.frame(fs <-effect("POR_FINAL_SIB", m_final))# we add the predicted log odds as well:fs_d$fit.lo <-qlogis(fs_d$fit)# for later, we add the proportions of values in the data:fs_d$PROP <- d$POR_FINAL_SIB %>% table %>% prop.table %>% as.numericfs_d
Figure 20: The effect of final sibilants of the possessor
I would probably plot that effect like this:
op <-par(mar=c(4, 4, 4, 4)) # make the plotting margins big enoughplot(pch=16, col="#00000060", axes=FALSE,cex=4* fs_d$PROP,xlab="Sibilant at end of possessor" , x=1:2,ylab="Pred. prob. of s-genitives", ylim=c(0, 1), y=fs_d$fit)abline(h=seq(0.1, 0.9, 0.2), lty=3, col="#00000030")axis(1, at=1:2, labels=levels(d$POR_FINAL_SIB))axis(2, at=seq(0.1, 0.9, 0.2)); abline(h=0.5, lty=3)axis(4, at=seq(0.1, 0.9, 0.2), labels=round(qlogis(seq(0.1, 0.9, 0.2)), 2))mtext("Pred. log odds", 4, 2.5)for (i inseq(nrow(fs_d))) {arrows(i, fs_d$lower[i], i, fs_d$upper[i], angle=90, code=3)}par(op)
Figure 21: The effect of final sibilants of the possessor
Interpretation: a clear articulatory effect: when the possessor ends in a sibilant, s-genitives are less likely.
1.2.3.3 Interaction 1: SPEAKER:POR_ANIMACY
Here’s the effects object for the first interaction, SPEAKER:POR_ANIMACY:
span_d <-data.frame(span <-effect("SPEAKER:POR_ANIMACY", m_final))# we add the predicted log odds as well:span_d$fit.lo <-qlogis(span_d$fit)# for later, we add the proportions of values in the data:span_d$PROP <- d$SPEAKER %>%table(d$POR_ANIMACY) %>% prop.table %>% as.numericspan_d
plot(span, type="response", ylim=c(0, 1), grid=TRUE,multiline=TRUE, confint=list(style="auto"))# the next plot is less useful, I think, but you should still have plotted it:# plot(span, type="response", ylim=c(0, 1), grid=TRUE,# x.var="SPEAKER", multiline=TRUE, confint=list(style="auto"))
Figure 22: The effect of Speaker interacting with Possessor animacy
Figure 23: The effect of Speaker interacting with Possessor animacy
Interpretation: Both NS and NNS are most likely to use s-genitives with animate possessors, followed by collective/temporal possessors, then locative, and finally inanimate ones (e.g., concrete objects), but the only significant difference in the interaction is that NS are significantly more likely to use s-genitives with collectives than NNS.
1.2.3.4 Interaction 2: SPEAKER:MODALITY
Here’s the effects object for the second interaction, SPEAKER:MODALITY:
spmo_d <-data.frame(spmo <-effect("SPEAKER:MODALITY", m_final))# we add the predicted log odds as well:spmo_d$fit.lo <-qlogis(spmo_d$fit)# for later, we add the proportions of values in the data:spmo_d$PROP <- d$SPEAKER %>%table(d$MODALITY) %>% prop.table %>% as.numericspmo_d
plot(spmo, type="response", ylim=c(0, 1), grid=TRUE,multiline=TRUE, confint=list(style="auto"))# the next plot is not any more useful, but you should still have plotted it:# plot(spmo, type="response", ylim=c(0, 1), grid=TRUE,# x.var="MODALITY", multiline=TRUE, confint=list(style="auto"))
Figure 24: The effect of Speaker interacting with Modality
Figure 25: The effect of Speaker interacting with Modality
Interpretation: NS are more likely to use s-genitives in writing than they are in speaking and than NNS are in writing. NNS are less likely to use s-genitives in writing than they are in speaking and than NS are in writing. The highest predicted probability of s-genitives is found for NNS in speaking and the lowest predicted probability of s-genitives is found for NS in speaking.
1.2.4 Write-up
To determine whether the choice of a genitive varies as a function of
whether the speaker is a native speaker or not;
the length difference between possessor and possessum and whether native and non-native speakers react to it in the same way;
the animacy of the possessor and whether native and non-native speakers react to it in the same way;
whether the possessor ends in a sibilant or not and whether native and non-native speakers react to it in the same way;
whether the modality is speaking or writing and whether native and non-native speakers react to it in the same way;
the interaction of POR_FINAL_SIB with MODALITY.
a generalized linear model selection process was undertaken (using backwards stepwise model selection of all predictors and controls based on significance testing). The final model resulting from the elimination of ns predictors contained the effects of LEN_PORmPUM_LOG:POR_ANIMACY and SPEAKER:POR_ANIMACY and was highly significant (LR=1596.2, df=14, p<0.0001) with a good amount of explanatory/predictive power (McFadden’s R2=0.397, Nagelkerke’s R2=0.534, C=0.899). The summary table of the model is shown here.
there is a clear main effect of LEN_PORmPUM_LOG amounting to a short-before-long effect;
there is a clear main effect of POR_FINAL_SIB: when the possessor ends in a sibilant, s-genitives are less likely;
the interaction SPEAKER:POR_ANIMACY only has the effect that native speakers are significantly more likely to use s genitives with collectives than learners are – all other speaker differences within each animacy level are ns. [Show plot(s)]
the interaction SPEAKER:MODALITY has the effect that NS are more likely to use s-genitives in writing than they are in speaking and than NNS are in writing. NNS are less likely to use s-genitives in writing than they are in speaking and than NS are in writing. [Show plot(s)]
---title: "Ling 204: session 01"author: - name: "[Stefan Th. Gries](https://www.stgries.info)" affiliation: - UC Santa Barbara - JLU Giessen orcid: 0000-0002-6497-3958date: "2026-01-07 12:34:56 PDT"date-format: "DD MMM YYYY HH-mm-ss"editor: sourceformat: html: page-layout: full code-fold: false code-link: true code-copy: true code-tools: true code-line-numbers: true code-overflow: scroll number-sections: true smooth-scroll: true toc: true toc-depth: 4 number-depth: 4 toc-location: left monofont: lucida console tbl-cap-location: top fig-cap-location: bottom fig-width: 4 fig-height: 4 fig-format: png fig-dpi: 300 fig-align: center embed-resources: true link-external-newwindow: trueexecute: cache: false echo: true eval: true warning: false---# Session 01: Recap of fixef regression modeling {#sec-s01}## Linear modeling {#sec-s01_lin}To do a quick recap of linear modeling, we will look at a data set in [_input/rts.csv](_input/rts.csv) involving reaction times as the response variable:```{r}rm(list=ls(all.names=TRUE))pkgs2load <-c("car", "effects", "emmeans", "magrittr", "multcomp")suppressMessages(sapply(pkgs2load, library,character.only=TRUE, logical.return=TRUE, quietly=TRUE)); rm(pkgs2load)summary(d <-read.delim( # summarize d, the result of loadingfile="_input/rts.csv", # this filestringsAsFactors=TRUE)) # change categorical variables into factors```You can see what the variables/columns contain in the file [_input/rts.r](_input/rts.r). Specifically, we are asking, does the (logged) reaction time to a word (in ms) vary as a function of* the Zipf frequency of that word (`ZIPFFREQ`);* the length of that word (`LENGTH`);* the language that word was presented in (`LANGUAGE`: *spanish* vs. *english*);* the speaker group that words was presented to (`GROUP`: *english* vs. *heritage*);* any pairwise interaction of these predictors? (Specifically, we are expecting that the effect of `LENGTH` will be stronger for when `LANGUAGE` is *spanish*.)### Exploration & preparation {#sec-s01_lin_explprep}As always, a decent amount of exploration is required to prepare the data for a hopefully unproblematic modeling process.#### The response variable {#sec-s01_lin_explprep_resp}We begin our exploration with the response variable `RT`. An initial histogram shows that we will probably want to transform the response given the long right tail. A log transformation seems to work well enough so that's what we'll go with:```{r}#| label: fig-RT#| fig-cap: The distribution of RT#| fig-width: 10#| fig-show: holdpar(mfrow=c(1, 2)) # make the plotting window have 1 rows & 2 columnshist(d$RT, main="")hist(log2(d$RT), main="")par(mfrow=c(1, 1)) # reset to defaultd$RT_log <-log2(d$RT)```#### Numeric predictors {#sec-s01_lin_explprep_numpred}Let's now look at the numeric predictors. We begin with `ZIPFFREQ`, which looks ok:```{r}#| label: fig-RT_n_ZIPFFREQ#| fig-cap: The distribution of ZIPFFREQ & RT#| fig-width: 10#| fig-show: holdpar(mfrow=c(1, 2)) # make the plotting window have 1 rows & 2 columnshist(d$ZIPFFREQ) # hist(d$ZIPFFREQ, breaks="FD") # usually, also check the bin widthplot(main="RT (in ms, logged) as a function of Zipf frequency",pch=16, col="#00000020",xlab="Zipf frequency" , x=d$ZIPFFREQ,ylab="RT (in ms, logged)", y=d$RT_log); grid()abline(lm(d$RT_log ~ d$ZIPFFREQ), col="green", lwd=2)par(mfrow=c(1, 1)) # reset to default```Let's now look at `LENGTH`, which looks ok:```{r}#| label: fig-RT_n_LENGTH#| fig-cap: The distribution of LENGTH & RT#| fig-width: 10#| fig-show: holdpar(mfrow=c(1, 2)) # make the plotting window have 1 rows & 2 columnshist(d$LENGTH) # hist(d$LENGTH, breaks="FD") # usually, also check the bin widthplot(main="RT (in ms, logged) as a function of Length",pch=16, col="#00000020",xlab="Length" , x=jitter(d$LENGTH),ylab="RT (in ms, logged)", y=d$RT_log); grid()abline(lm(d$RT_log ~ d$LENGTH), col="green", lwd=2)par(mfrow=c(1, 1)) # reset to default```#### Categorical predictors {#sec-s01_lin_explprep_catpred}Let's now look at the categorical predictors. We begin with `LANGUAGE`, which looks ok:```{r}#| label: fig-RT_n_LANGUAGE#| fig-cap: The distribution of LANGUAGE & RT#| fig-show: holdqwe <-boxplot(main="RT (in ms, logged) as a function of stimulus language", d$RT_log ~ d$LANGUAGE, notch=TRUE,xlab="Stimulus language",ylab="RT (in ms, logged)"); grid()```Let's now look at `GROUP`, which looks ok:```{r}#| label: fig-RT_n_GROUP#| fig-cap: The distribution of GROUP & RT#| fig-show: holdqwe <-boxplot(main="RT (in ms, logged) as a function of speaker group", d$RT_log ~ d$GROUP, notch=TRUE,xlab="Speaker group",ylab="RT (in ms, logged)"); grid()```Do we have a good distribution for the interaction of the categorical variables? Yes:```{r}table(d$LANGUAGE, d$GROUP)```#### Deviance {#sec-s01_lin_explprep_dev}First, we'll quickly compute the deviance of the null model:```{r}deviance(m_00 <-lm(RT_log ~1, data=d)) # 1581.951```### Model selection {#sec-s01_lin_modeling}The initial model involves all predictors and all their interactions:```{r}summary(m_01 <-lm( # make/summarize the linear model m_01: RT_log ~1+# RT_log ~ an overall intercept (1)# & these predictors & their pairwise interactions (ZIPFFREQ+LANGUAGE+GROUP+LENGTH)^2,data=d, # those vars are in dna.action=na.exclude)) # skip cases with NA/missing data (good habit)drop1(m_01, # drop each droppable predictor at a time from m_01 &test="F") %>%# test its significance w/ an F-test"["(order(.[,6]),) # order by p-values```Let's fit the 2nd model with that interaction `ZIPFFREQ:GROUP` deleted:```{r}summary(m_02 <-update( # make m_02 an update of m_01, .~. # m_01, namely all of it (.~.), but then- ZIPFFREQ:GROUP)) # remove this interactionanova(m_01, m_02, # compare m_01 to m_02test="F") # using an F-testdrop1(m_02, # drop each droppable predictor at a time from m_02 &test="F") %>%# test its significance w/ an F-test"["(order(.[,6]),) # order by p-values```Let's fit the 3rd model with the ns interaction `GROUP:LENGTH` deleted:```{r}summary(m_03 <-update( # make m_03 an update of m_02, .~. # m_02, namely all of it (.~.), but then- GROUP:LENGTH)) # remove this interactionanova(m_02, m_03, # compare m_02 to m_03test="F") # using an F-testdrop1(m_03, # drop each droppable predictor at a time from m_03 &test="F") %>%# test its significance w/ an F-test"["(order(.[,6]),) # order by p-values```Let's fit the 4th model with the ns interaction `ZIPFFREQ:LANGUAGE` deleted:```{r}summary(m_04 <-update( # make m_04 an update of m_03, .~. # m_03, namely all of it (.~.), but then- ZIPFFREQ:LANGUAGE)) # remove this interactionanova(m_03, m_04, # compare m_03 to m_04test="F") # using an F-testdrop1(m_04, # drop each droppable predictor at a time from m_04 &test="F") %>%# test its significance w/ an F-test"["(order(.[,6]),) # order by p-values```We're done now -- why? Because remember that we said "we are expecting that the effect of `LENGTH` will be stronger for when `LANGUAGE` is *spanish*", meaning if the coefficient for `LANGUAGEspanish:LENGTH` is positive -- which it is -- we can take half its *p*-value and then it is in fact significant! We call the final model `m_final` (for easy copying and pasting across analyses) and do some housekeeping:```{r}m_final <- m_04; rm(m_04); invisible(gc())```Normally, we would do some regression diagnostics **now** -- I have to skip those now to discuss how to interpret the model, but we will return to this later! For now, we compute the confidence intervals ...```{r}Confint(m_final)```... , add the predictions of the model to the original data frame:```{r}d$PRED <-# make d$PRED thepredict( # predictions for all cases m_final) # based on m_finalhead(d) # check the result```..., and generate an often super-useful kind of data frame that I always call `nd` (for *new data*) that contains all combinations of* all levels of all categorical predictors and* useful values of all numeric predictors (where 'useful' is usually 0, 1 and the mean):```{r}nd <-expand.grid(ZIPFFREQ=c(0, 1, mean(d$ZIPFFREQ)),LENGTH =c(0, 1, mean(d$LENGTH)),GROUP =levels(d$GROUP),LANGUAGE=levels(d$LANGUAGE))```To that, we add the predictions the model makes for all these combinations and give it some nice(r) names:```{r}nd <-cbind(nd,predict(m_final, newdata=nd, interval="confidence"))names(nd)[5:7] <-c("PRED", "PRED_LOW", "PRED_UPP")head(nd)```A maybe (!) nicer way to represent this might be this, but whether this is nicer or not is probably subjective:```{r}with(nd, # use the data in nd andtapply( # apply PRED, # to the predicted reaction times# a grouping like this: the numerics in the rows & columns, the factors in the layers/sliceslist(round(ZIPFFREQ, 2), round(LENGTH, 2), "GROUP"=GROUP, "LANGUAGE"=LANGUAGE),# just list the values, but rounded to 2 decimals c)) %>%round(2)```### Model evaluation & interpretation {#sec-s01_lin_modeleval}Let's look at the results.#### `LANGUAGE:GROUP` {#sec-s01_lin_modeleval_catcat}Let's visualize the nature of the effect of `LANGUAGE:GROUP` based on the predictions from `effects`:```{r}#| label: fig-LAGR#| fig-cap: The interaction LANGUAGE:GROUPlagr_d <-data.frame( # make lagr_d a data frame of lagr <-effect( # the effect"LANGUAGE:GROUP", # of LANGUAGE:GROUP m_final)) # in m_finalplot(lagr, # plot the effect of the interactionylim=range(d$RT_log, na.rm=TRUE), # w/ these y-axis limitsgrid=TRUE) # & w/ a gridplot(lagr, # plot the effect of the interactionylim=range(d$RT_log, na.rm=TRUE), # w/ these y-axis limitsx.var="GROUP", # w/ this predictor on the x-axisgrid=TRUE) # & w/ a grid```And here's a version (of the first of the two effects plots) with everything: observed values, predicted values, and the confidence intervals of the predictions:```{r}#| label: fig-LAGR_better#| fig-cap: The interaction LANGUAGE:GROUP#| fig-width: 10#| fig-show: hold# split data & predictions up by levels of GROUP:d_split <-split(d, d$GROUP)lagr_d_split <-split(lagr_d, lagr_d$GROUP, drop=TRUE)par(mfrow=c(1, 2))stripchart(main="RT ~ LANGUAGE for GROUP='english'", # stripchart w/ main headingpch=16, col="#00000020", ylim=c(8, 12), # filled semi-transparent grey circles, y limits d_split$english$RT_log ~ d_split$english$LANGUAGE, # data: RT_log ~ LANG (for eng)xlab="Language", ylab="RT (in ms, logged)", # axis labelsvertical=TRUE, method="jitter"); grid() # vertical plot & jittered observations, add grid# add predicted means:points(seq(lagr_d_split$english$fit), lagr_d_split$english$fit, pch="X", col=c("red", "blue"))# add confidence intervals:for (i inseq(lagr_d_split$english$fit)) { # as many arrows as neededarrows(i, lagr_d_split$english$lower[i], # each from here i, lagr_d_split$english$upper[i], # to herecode=3, angle=90, col=c("red", "blue")[i]) # both tips w/ 90 degrees}stripchart(main="RT ~ LANGUAGE for GROUP='heritage'", # stripchart w/ main headingpch=16, col="#00000020", ylim=c(8, 12), # filled semi-transparent grey circles, y limits d_split$heritage$RT_log ~ d_split$heritage$LANGUAGE, # data: RT_log ~ LANG (for her)xlab="Language", ylab="RT (in ms, logged)", # axis labelsvertical=TRUE, method="jitter"); grid() # vertical plot & jittered observations, add grid# add predicted means:points(seq(lagr_d_split$heritage$fit), lagr_d_split$heritage$fit, pch="X", col=c("red", "blue"))# add confidence intervals:for (i inseq(lagr_d_split$heritage$fit)) { # as many arrows as neededarrows(i, lagr_d_split$heritage$lower[i], # each from here i, lagr_d_split$heritage$upper[i], # to herecode=3, angle=90, col=c("red", "blue")[i]) # both tips w/ 90 degrees}par(mfrow=c(1, 1))```#### `LANGUAGE:LENGTH` {#sec-s01_lin_modeleval_catnum}Let's visualize the nature of the effect of `LANGUAGE:LENGTH` based on the predictions from `effects`:```{r}#| label: fig-LALE#| fig-cap: The interaction LANGUAGE:LENGTH#| fig-width: 10#| fig-show: holdlale_d <-data.frame( # make lale_d a data frame of lale <-effect( # the effect"LANGUAGE:LENGTH", # of LANGUAGE:LENGTH m_final, # in m_finalxlevels=list( # for whenLENGTH=3:9))) # LENGTH is these valueslale_d <- lale_d[order(lale_d$LANGUAGE, lale_d$LENGTH),]; lale_dplot(lale, # plot the effect of the interactionylim=range(d$RT_log, na.rm=TRUE), # w/ these y-axis limitsgrid=TRUE) # & w/ a grid```And here a version with everything: observed values, predicted values (red for `LANGUAGE`: *english*, blue for `LANGUAGE`: *spanish*), and the confidence intervals of the predictions (I'm building this up very stepwise, more so than I might do in actual research):```{r}#| label: fig-LALE_better#| fig-cap: The interaction LANGUAGE:LENGTHd_split <-split(d, d$LANGUAGE)lale_d_split <-split(lale_d, lale_d$LANGUAGE, drop=TRUE)plot(main="RT ~ LANGUAGE:LENGTH", type="n", # plot nothingxlab="Length" , x=d$LENGTH, # x-axisylab="RT (in ms, logged)", y=d$RT_log); grid() # y-axis & grid# plot observed data points for LANGUAGE being eng:points(jitter(d_split$eng$LENGTH), d_split$eng$RT_log, pch=16, col="#FF000010")lines(lale_d_split$eng$LENGTH, lale_d_split$eng$fit, col="red")# plot observed data points for LANGUAGE being spa:points(jitter(d_split$spa$LENGTH), d_split$spa$RT_log, pch=16, col="#0000FF10")lines(lale_d_split$spa$LENGTH, lale_d_split$spa$fit, col="blue")# plot confidence band for LANGUAGE being eng:polygon(col="#FF000030", border="red", # semi-transparent grey, for once w/ borderx=c( lale_d_split$eng$LENGTH, # x-coordinates left to right, thenrev(lale_d_split$eng$LENGTH)), # the same values back (in reverse order)y=c( lale_d_split$eng$lower, # y-coordinates lower boundaryrev(lale_d_split$eng$upper))) # y-coordinates upper boundary# plot confidence band for LANGUAGE being spa:polygon(col="#0000FF30", border="blue", # semi-transparent grey, for once w/ borderx=c( lale_d_split$spa$LENGTH, # x-coordinates left to right, thenrev(lale_d_split$spa$LENGTH)), # the same values back (in reverse order)y=c( lale_d_split$spa$lower, # y-coordinates lower boundaryrev(lale_d_split$spa$upper))) # y-coordinates upper boundary# plot legendlegend(title="Language", bty="n",x=6, xjust=0.5, y=11.5, yjust=0.5, ncol=2,legend=levels(d$LANGUAGE), fill=c("red", "blue"))``````{r}#| echo: false#| eval: false# These are the slopes of LENGTH one might report,# which differ by coef(m_final)[8]emmeans::emtrends( # of trend lines (i.e. slopes of regression lines)object=m_final, # for the model m_finalvar="LENGTH", # for this numeric predictorspecs=~ LANGUAGE) # for each level of this categorical predictor```#### `ZIPFFREQ:LENGTH` {#sec-s01_lin_modeleval_numnum}Let's visualize the nature of the effect of `ZIPFFREQ:LENGTH` based on the predictions from `effects`:```{r}#| label: fig-ZILE#| fig-cap: The interaction ZIPFFREQ:LENGTH(zile_d <-data.frame( # make zile_d a data frame of zile <-effect( # the effect"ZIPFFREQ:LENGTH", # of ZIPFFREQ:LENGTH m_final, # in m_finalxlevels=list( # for whenZIPFFREQ=seq(3, 6, 0.25), # ZIPFFREQ is these valuesLENGTH=seq(3, 9, 0.5))))) # LENGTH is these valuesplot(zile, # plot the effect of the interactionylim=range(d$RT_log, na.rm=TRUE), # w/ these y-axis limitsgrid=TRUE) # & w/ a grid```Not that great ... There are two main ways we could try and plot this better. One would be with an actual 3-D plot. The simplest version of this might use `rgl::plot3D` like this, which would open a new window with a rotatable 3-D plot (but this wouldn't be rendered into this HTML properly, which is why I am not executing it here):```{r}#| eval: falsergl::plot3d(zile_d$ZIPFFREQ, zile_d$LENGTH, zile_d$fit)```Another option would be to use `plotly::plot_ly`:```{r}library(plotly)fig <-plot_ly(zile_d,x=~ZIPFFREQ, y=~LENGTH, z=~fit, color=~fit, colors=c("#FFFFFF20", "#00000020"))fig <- fig %>%add_markers()fig <- fig %>%layout(scene =list(xaxis =list(title ='Zipf frequency'),yaxis =list(title ='Length'),zaxis =list(title ='Predicted RT')))fig```While this is nice, it's also hard to publish in print. Let's do this differently with a numeric heatmap kind of plot that I often use. For that, the numeric predictions are first cut into 9 ordinal bins with lower/higher bin numbers reflecting lower/higher numeric predictions ...```{r}zile_d$fitcat <-as.numeric( # create a numeric column fitcat in zile_dcut(zile_d$fit, # by cutting the predicted values upbreaks=seq( # at the values ranging frommin(zile_d$fit), # the smallest obs/pred RT tomax(zile_d$fit), # the largest obs/pred RTlength.out=10), # into 9 groupsinclude.lowest=TRUE, # include the lowest, andlabels=1:9)) # use the numbers 1:9 as labels```... and then we plot them into a coordinate system:```{r}#| label: fig-ZILE_better#| fig-cap: The interaction ZIPFFREQ:LENGTHplot(main="RT ~ ZIPFFREQ:LENGTH", type="n", # plot nothing, justxlab="Zipf frequency", xlim=range(zile_d$ZIPFFREQ), x=0, # get the coordinateylab="Length" , ylim=range(zile_d$LENGTH) , y=0) # system set uptext( # plot textx=zile_d$ZIPFFREQ, y=zile_d$LENGTH, # at these coordinates, namelylabels=zile_d$fitcat, # the 'categorical/ordinal' predictions# make the physical size of the numbers dependent on the valuescex=0.5+zile_d$fitcat/7,# make the darkness of the grey reflect the predicted RT_logcol=grey(level=zile_d$fitcat/10)) # see ?grey```However, there's one thing that's possibly problematic about this -- what is that? The problem is that we are binning with `cut` on the basis of `zile_d$fit` (check the above code), but the range of those predictions is much smaller than the range of the actually observed reaction times in `d$RT_log`:```{r}range(zile_d$fit)range(d$RT_log)```That in turn means the binning steps based on `zile_d$fit` will be based on much smaller differences between predicted values (of only 40, check `levels(cut(zile_d$fit, breaks=seq(min(zile_d$fit), max(zile_d$fit), length.out=10)))`), meaning a difference of predictions of only 81 ms could lead to a difference in 3 bins -- if we were to bin based on the actual value range, that would not happen, because there the binning steps are 429 ms wide (check `levels(cut(d$RT_log, breaks=seq(min(d$RT_log), max(d$RT_log), length.out=10)))`). Thus, the above plot exaggerates the effects somewhat.Accordingly, we determine the range of the values to be plotted in a potentially more representative way, by using* as a minimum, the minimum of the predictions *and * the actual values combined;* as a maximum, the maximum of the predictions *and * the actual values combined:```{r}# define the minimum of the range to be used for binningmin_of_obs_n_pred <-min(c( d$RT_log, zile_d$fit),na.rm=TRUE)# define the maximum of the range to be used for binningmax_of_obs_n_pred <-max(c( d$RT_log, zile_d$fit),na.rm=TRUE)```We can see that the initial plot indeed exaggerated the range: the upper plot showed prediction bins from 1 to 9, but we now see that the predicctions are actually only in a really small range of values.```{r}zile_d$fitcat2 <-as.numeric( # create a numeric column fitcat2 in zile_dcut(zile_d$fit, # by cutting the predicted values upbreaks=seq( # at the values ranging from min_of_obs_n_pred, # the smallest obs/pred RT_log to <--------- max_of_obs_n_pred, # the largest obs/pred RT_log <---------length.out=10), # into 9 groupsinclude.lowest=TRUE, # include the lowest, andlabels=1:9)) # use the numbers 1:9 as labelstable(zile_d$fitcat, zile_d$fitcat2)```Thus, we plot this again (with some additional info in the form of the decile-based grid), which here doesn't help much:```{r}#| label: fig-ZILE_evenbetter#| fig-cap: The interaction ZIPFFREQ:LENGTHplot(main="RT ~ ZIPFFREQ:LENGTH", type="n", # plot nothing, justxlab="Zipf frequency", xlim=range(zile_d$ZIPFFREQ), x=0, # get the coordinateylab="Length" , ylim=range(zile_d$LENGTH) , y=0) # system set uptext( # plot textx=zile_d$ZIPFFREQ, y=zile_d$LENGTH, # at these coordinates, namelylabels=zile_d$fitcat2, # the 'categorical/ordinal' predictions# make the physical size of the numbers dependent on the valuescex=0.5+zile_d$fitcat2/7,# make the darkness of the grey reflect the predicted RT_logcol=grey(level=zile_d$fitcat2/15)) # see ?greyabline(h=quantile(d$LENGTH , probs=seq(0, 1, 0.1)), lty=3, col="#00000020")abline(v=quantile(d$ZIPFFREQ, probs=seq(0, 1, 0.1)), lty=3, col="#00000020")``````{r}#| echo: false#| eval: false# LANGUAGE:GROUP# scenario 1a: what are the meansnd[nd$ZIPFFREQ==0& nd$LENGTH==0,]# scenario 1b: are they different from 0?c(1,0,0,0,0,0,0,0) %>%matrix(nrow=1) %>%glht(m_final, .) %>% summary # LANG:engl & GRP:englc(1,0,0,1,0,0,0,0) %>%matrix(nrow=1) %>%glht(m_final, .) %>% summary # LANG:engl & GRP:heric(1,0,1,0,0,0,0,0) %>%matrix(nrow=1) %>%glht(m_final, .) %>% summary # LANG:span & GRP:englc(1,0,1,1,0,0,1,0) %>%matrix(nrow=1) %>%glht(m_final, .) %>% summary # LANG:span & GRP:heri# scenario 1c: are they different from each other?pairs(emmeans::emmeans(object=m_final, ~ GROUP | LANGUAGE, at=list(ZIPFFREQ=0, LENGTH=0)), adjust="none")pairs(emmeans::emmeans(object=m_final, ~ LANGUAGE | GROUP, at=list(ZIPFFREQ=0, LENGTH=0)), adjust="none")pairs(emmeans::emmeans(object=m_final, ~ LANGUAGE:GROUP, at=list(ZIPFFREQ=0, LENGTH=0)), adjust="none")# without at: when ZIPFFREQ and LENGTH are at mean for the variable LENGTH interacting with them # or of course long glht comparisons# LANGUAGE:LENGTHnd[nd$ZIPFFREQ %in%0:1& nd$LENGTH %in%0:1,]# scenario 3a: what are the slopes of LENGTH for each level of LANGUAGE# scenario 3b: are they different from 0?c(0,0,0,0,1,0,0,0) %>%matrix(nrow=1) %>%glht(m_final, .) %>% summary # LANGUAGE: english & GROUP:englishc(0,0,0,0,1,0,0,1) %>%matrix(nrow=1) %>%glht(m_final, .) %>% summary # LANGUAGE: spanish & GROUP:english# scenario 3c: are they different from each other?c(0,0,0,0,0,0,0,1) %>%matrix(nrow=1) %>%glht(m_final, .) %>% summarypairs(emmeans::emtrends(object=m_final, ~ LANGUAGE, var="LENGTH"), adjust="none")(c(0,0,0,0,1,0,0,0) -c(0,0,0,0,1,0,0,1)) %>%matrix(nrow=1) %>%glht(m_final, .) %>% summary # LANGUAGE: english & GROUP:english# ZIPFFREQ:LENGTHnd[nd$GROUP=="english"& nd$LANGUAGE=="english"& nd$ZIPFFREQ %in%0:1& nd$LENGTH %in%0:1,]# scenario 3b: are all the following different from 0?# scenario 3a: what are the slopes of LENGTH when ZIPFFREQ=0?c(0,0,0,0,1,0,0,1) %>%matrix(nrow=1) %>%glht(m_final, .) %>% summary # LANGUAGE: english & GROUP:english# scenario 3a: what are the slopes of LENGTH when ZIPFFREQ=1?c(0,0,0,0,1,1,0,1) %>%matrix(nrow=1) %>%glht(m_final, .) %>% summary # LANGUAGE: english & GROUP:english# scenario 3a: what are the slopes of ZIPFFREQ when LENGTH=0?c(0,1,0,0,0,0,0,0) %>%matrix(nrow=1) %>%glht(m_final, .) %>% summary # LANGUAGE: english & GROUP:english# scenario 3a: what are the slopes of ZIPFFREQ when LENGTH=1?c(0,1,0,0,0,1,0,0) %>%matrix(nrow=1) %>%glht(m_final, .) %>% summary # LANGUAGE: english & GROUP:english# scenario 3b: are they different from each other?c(0,0,0,0,0,1,0,0) %>%matrix(nrow=1) %>%glht(m_final, .) %>% summary # LANGUAGE: english & GROUP:english```### Write-up {#sec-s01_lin_writeup}To determine whether the reaction time to a word (in ms) varies as a function of* the Zipf frequency of that word (`ZIPFFREQ`);* the length of that word (`LENGTH`);* the language that word was presented in (`LANGUAGE`: *english* vs. *spanish*);* the speaker group that words was presented to (`GROUP`: *english* vs. *heritage*);* any pairwise interaction of these predictors with the specific expectation that the effect of `LENGTH` will be stronger for when `LANGUAGE` is *spanish*a linear model selection process was undertaken. In order to deal with the very long right tail of the response variable the response variable was logged; then a backwards stepwise model selection of all predictors and controls (based on significance testing) was applied. The final model resulting from the elimination of ns predictors contained the effects of `ZIPFFREQ:LENGTH`, `LANGUAGE:GROUP`, and `LANGUAGE:LENGTH` and was highly significant (*F*=154.7, *df*~1~=7, *df*~2~=7744, *p*<0.0001), but it explains only little of the response variable (mult. *R*^2^=0.123, adj. *R*^2^=0.122). The summary table of the model is shown here.| | Estimate | 95%-CI | *se* | *t* | *p*~2-tailed~ ||:----------------------------------------------------------------------|:---------|:-----------------|:------|:--------|:--------------|| Intercept (ZIPFFREQ=0, LENGTH=0, LANGUAGE=*spanish*, GROUP=*english*) | 9 | [8.595, 9.409] | 0.208 | 43.33 | <0.001 || ZIPFFREQ~0→1~ | 0.007 | [-0.077, 0.092] | 0.043 | 0.171 | 0.864 || LANGUAGE~*spanish*→*english*~ | 0.258 | [0.149, 0.366] | 0.056 | 4.638 | <0.001 || GROUP~*english*→*heritage*~ | 0.103 | [0.076, 0.129] | 0.013 | 7.655 | <0.001 || LENGTH~0→1~ | 0.113 | [0.037, 0.19] | 0.039 | 2.913 | 0.004 || ZIPFFREQ~0→1~ + LENGTH~0→1~ | -0.021 | [-0.037 -0.006] | 0.008 | -2.648 | 0.008 || LANGUAGE~*spanish*→*english*~ + GROUP~*english*→*heritage*~ | -0.205 | [-0.243, -0.168] | 0.019 | -10.668 | <0.001 || LANGUAGE~*spanish*→*english*~ + LENGTH~0→1~ | 0.018 | [-0.002, 0.039] | 0.01 | 1.754 | 0.079 |The results indicate that* as hypothesized, the slowing-down effect of `LENGTH` is stronger for Spanish than for English words [show plot];* the interaction `LANGUAGE:GROUP` has the following effects [show at least one plot]: + reaction times are always higher when `LANGUAGE` is *spanish* than when `LANGUAGE` is *english*, but that difference is higher when `GROUP` is *english* than when `GROUP` is *spanish*; + when `LANGUAGE` is *english*, heritage speakers need more time than learners, but when `LANGUAGE` is *spanish*, heritage speakers need less time than learners; + all 4 means of `LANGUAGE:GROUP` are significantly different from each other (non-overlapping CIs);* the interaction `LANGUAGE:LENGTH` has the hypothesized effect that the effect of `LENGTH` is stronger -- in fact, nearly 2.25 times as strong -- when `LANGUAGE` is *spanish* than when `LANGUAGE` is *english* [show at least one plot];* the interaction of `ZIPFFREQ:LENGTH` is not strong: both variables behave as expected -- `LENGTH` slows down, `ZIPFFREQ` speeds up -- but + the effect of `LENGTH` is much stronger with rarer words than with frequent words; + the effect of `ZIPFFREQ` is much stronger with longer words than with shorter words.## Binary logistic regression modeling {#sec-s01_binlog}To do a quick recap of linear modeling, we will look at a data set in [_input/genitives.csv](_input/genitives.csv) involving genitive choices as the response variable:```{r}rm(list=ls(all.names=TRUE))pkgs2load <-c("car", "effects", "emmeans", "magrittr", "multcomp")suppressMessages(sapply(pkgs2load, library,character.only=TRUE, logical.return=TRUE, quietly=TRUE)); rm(pkgs2load)source("_helpers/C.score.r"); source("_helpers/R2s.r"); # helper functionssummary(d <-read.delim( # summarize d, the result of loadingfile="_input/genitives.csv", # this filestringsAsFactors=TRUE)) # change categorical variables into factors```You can see what the variables/columns contain in the file [_input/genitives.r](_input/genitives.r). Specifically, we want to see whether the choice of a genitive varies as a function of* whether the speaker is a native speaker or not (`SPEAKER`);* the length difference between possessor and possessum and whether native and non-native speakers react to it in the same way;* the animacy of the possessor and whether native and non-native speakers react to it in the same way (`POR_ANIMACY`);* whether the possessor ends in a sibilant or not and whether native and non-native speakers react to it in the same way (`POR_FINAL_SIB`);* whether the modality is speaking or writing and whether native and non-native speakers react to it in the same way (`MODALITY`);* the interaction of `POR_FINAL_SIB` with `MODALITY`.### Exploration & preparation {#sec-s01_binlog_explprep}#### All variables {#sec-s01_binlog_explprep_vars}As always, a decent amount of exploration is required to prepare the data for a hopefully unproblematic modeling process. We will only focus on the four variables of interest (and even there will take some didactically motivated shortcuts). We begin our exploration with the creation of the length difference predictor:Let's look at the numeric predictor first by exploring its 'components':```{r}#| label: fig-lengths#| fig-cap: The lengths of possessors/-umspar(mfrow=c(1, 2)) # make the plotting window have 1 rows & 2 columnshist(d$POR_LENGTH, main="")hist(d$PUM_LENGTH, main="")par(mfrow=c(1, 1)) # reset to default```This doesn't look all that great because it's not likely that the long tails will magically cancel each other out, and they don't:```{r}#| label: fig-lengthdiff1#| fig-cap: The difference of the logged lengthspar(mfrow=c(1, 2)) # make the plotting window have 1 rows & 2 columnshist(d$POR_LENGTH - d$PUM_LENGTH, main="")hist(d$POR_LENGTH - d$PUM_LENGTH, breaks="FD", main="")par(mfrow=c(1, 1)) # reset to default```What if we compute not the difference of the lengths, but the difference of the logged lengths?```{r}#| label: fig-lengthdiff#| fig-cap: The difference of the logged lengthshist(d$LEN_PORmPUM_LOG <-log2(d$POR_LENGTH)-log2(d$PUM_LENGTH), main="")```That looks awesome -- how this is related to the response? Very nicely:```{r}#| label: fig-lengthdiff_gen#| fig-cap: The difference of the logged lengths and genitivesspineplot(d$GENITIVE ~ d$LEN_PORmPUM_LOG)```Then we cross-tabulate for all interactions:```{r}cut(d$LEN_PORmPUM_LOG, 8) %>%table(d$SPEAKER) %>% addmarginsd$POR_ANIMACY %>%table(d$SPEAKER) %>% addmarginsd$POR_FINAL_SIB %>%table(d$SPEAKER) %>% addmarginsd$MODALITY %>%table(d$SPEAKER) %>% addmargins```All looking good, but what if we add the response?```{r}cut(d$LEN_PORmPUM_LOG, 8) %>%ftable(d$SPEAKER, d$GENITIVE) d$POR_ANIMACY %>%ftable(d$SPEAKER, d$GENITIVE) d$POR_FINAL_SIB %>%ftable(d$SPEAKER, d$GENITIVE) d$MODALITY %>%ftable(d$SPEAKER, d$GENITIVE)```Not too bad. There are some low-frequency cells and for the numeric predictor even some empty cells, but that problem will not surface in the model because there we won't use a factorized version of that predictor.#### Deviance {#sec-s01_binlog_explprep_dev}First, we'll quickly compute the deviance of the null model:```{r}deviance(m_00 <-glm(GENITIVE ~1, family=binomial, data=d)) # 4004.273```### Model selection {#sec-s01_binlog_modeling}Let's fit our initial regression model:```{r}summary(m_01 <-glm( # make/summarize the gen. linear model m_01: GENITIVE ~1+# GENITIVE ~ an overall intercept (1)# these predictors & their interaction, as specified above: SPEAKER * (LEN_PORmPUM_LOG + POR_ANIMACY + POR_FINAL_SIB + MODALITY) + POR_FINAL_SIB:MODALITY,family=binomial, # resp = binarydata=d, # those vars are in dna.action=na.exclude)) # skip cases with NA/missing datapchisq( # compute the area under the chi-squared curveq=m_01$null.deviance-m_01$deviance, # for this chi-squared value: null - res. dev.df=m_01$df.null-m_01$df.residual, # at this df: null - residual dflower.tail=FALSE) # only the right/upper tail/sidedrop1(m_01, # drop each droppable predictor at a time from m_01 &test="Chisq") %>%# test its significance w/ a LRT"["(order(.[,4]),) # order by p-values```Let's fit the 2nd model with the 2-way interaction `SPEAKER:LEN_PORmPUM_LOG` deleted:```{r}summary(m_02 <-update( # make m_02 an update of m_01, .~. # m_01, namely all of it (.~.), but then- SPEAKER:LEN_PORmPUM_LOG)) # remove this interactionanova(m_01, m_02, # compare m_01 to m_02test="Chisq") # using a LRTdrop1(m_02, # drop each droppable predictor at a time from m_02 &test="Chisq") %>%# test its significance w/ a LRT"["(order(.[,4]),) # order by p-values```Let's fit the 3rd model with the 2-way interaction `SPEAKER:POR_FINAL_SIB` deleted:```{r}summary(m_03 <-update( # make m_03 an update of m_02, .~. # m_02, namely all of it (.~.), but then- SPEAKER:POR_FINAL_SIB)) # remove this interactionanova(m_02, m_03, # compare m_02 to m_03test="Chisq") # using a LRTdrop1(m_03, # drop each droppable predictor at a time from m_03 &test="Chisq") %>%# test its significance w/ a LRT"["(order(.[,4]),) # order by p-values```Let's fit the 4th model with the 2-way interaction `POR_FINAL_SIB:MODALITY` deleted:```{r}summary(m_04 <-update( # make m_04 an update of m_03, .~. # m_03, namely all of it (.~.), but then- POR_FINAL_SIB:MODALITY)) # remove this interactionanova(m_03, m_04, # compare m_03 to m_04test="Chisq") # using a LRTdrop1(m_04, # drop each droppable predictor at a time from m_04 &test="Chisq") %>%# test its significance w/ a LRT"["(order(.[,4]),) # order by p-values```It seems like we're done:```{r}m_final <- m_04; rm(m_04); invisible(gc())```Here are the confidence intervals; in the interest of time, I am skipping `nd`:```{r}Confint(m_final)```### Model evaluation & interpretation: numerical & visual {#sec-s01_binlog_modeleval}We first assess the predictive power of the model in various ways. For that, we first compute all numeric and categorical predictions (as well as `PRED_PP_OBS`, which is the predicted probability of the level of the response variable that *should* have been predicted) ...```{r}d$PRED_PP_2 <-predict( # make d$PRED_PP_2 the result of predicting m_final, # from m_finaltype="response") # predicted probabilities of "s"d$PRED_LO_2 <-predict( # make d$PRED_LO_2 the result of predicting m_final) # from m_final -- this time the log oddsd$PRED_CAT <-factor( # d$PRED_CAT is determined by ifelseifelse(d$PRED_PP_2>=0.5, # if the predicted prob. of "s" is >=0.5levels(d$GENITIVE)[2], # predict "s"levels(d$GENITIVE)[1]), # otherwise, predict "of"levels=levels(d$GENITIVE)) # make sure both levels are registeredd$PRED_PP_OBS <-# d$PRED_PP_OBS is determined by ifelseifelse(d$GENITIVE==levels(d$GENITIVE)[2], # if the obs gen = 2nd level of GENITIVE d$PRED_PP_2, # take its predicted prob.1-d$PRED_PP_2) # otherwise take 1 minus its predicted probabilityhead(d)```... and then we evaluate the confusion matrix:```{r}(c_m <-table( # confusion matrix: cross-tabulateOBS =d$GENITIVE, # observed choices in the rowsPREDS=d$PRED_CAT)) # predicted choices in the columnsc( # evaluate the confusion matrix"Class. acc."=mean(d$GENITIVE==d$PRED_CAT, na.rm=TRUE), # (tp+tn) / (tp+tn+fp+fn)"Prec. for s"=c_m["s","s"] /sum(c_m[ ,"s" ]),"Acc./rec. for s"=c_m["s","s"] /sum(c_m["s" , ]),"Prec. for of"=c_m["of","of"] /sum(c_m[ ,"of"]),"Acc./rec. for of"=c_m["of","of"] /sum(c_m["of", ]))```Also, we compute *R*^2^-values, Cohen's κ, and *C*:```{r}c(R2s(m_final), "Cohen's kappa"=cohens.kappa(c_m)[[1]],"C"=C.score(cv.f=d$GENITIVE, d$PRED_PP_2))```Let's now look at the effects/interactions one by one.#### Main effect 1: `LEN_PORmPUM_LOG`Here's the effects object for the first interaction, `LEN_PORmPUM_LOG`:```{r}ld_d <-data.frame(ld <-effect("LEN_PORmPUM_LOG", m_final,xlevels=list(LEN_PORmPUM_LOG=seq(-5, 5, 1))))# we add the predicted log odds as well:ld_d$fit.lo <-qlogis(ld_d$fit)# for later, we add the proportions of values in the data:ld_d$PROP <- d$LEN_PORmPUM_LOG %>%cut(breaks=seq(-5.5, 5.5, 1), include.lowest=TRUE) %>% table %>% prop.table %>% as.numericld_d %>%round(4)```Let's visualize:```{r}#| label: fig-LD#| fig-cap: The effect of length differences# plot(ld, grid=TRUE)plot(ld, type="response", ylim=c(0, 1), grid=TRUE)```I might plot that effect like this:```{r}#| label: fig-LD_better#| fig-cap: The effect of length differencesop <-par(mar=c(4, 4, 4, 4)) # make the plotting margins big enoughplot(pch=16, col="#00000060", type="b", yaxp=c(0.1, 0.9, 4),cex=0.5+ (10* ld_d$PROP),xlab="Difference of logged lengths (or-um)" , x=ld_d$LEN_PORmPUM_LOG,ylab="Pred. prob. of s-genitives" , ylim=c(0, 1), y=ld_d$fit)abline(h=seq(0.1, 0.9, 0.2), lty=3, col="#00000030")abline(v=quantile(d$LEN_PORmPUM_LOG, prob=seq(0, 1, 0.1)), lty=3, col="#00000030")polygon(border=NA, col="#00000020",c(ld_d$LEN_PORmPUM_LOG, rev(ld_d$LEN_PORmPUM_LOG)),c(ld_d$lower , rev(ld_d$upper)))axis(4, at=seq(0.1, 0.9, 0.2), labels=round(qlogis(seq(0.1, 0.9, 0.2)), 2))mtext("Pred. log odds", 4, 2.5)par(op) # reset to defaults```Interpretation: a very clear short-before-long effect: when the possessor is longer than the possessum (on the right), *s*-genitives, which would have the long possessor in the front, are very unlikely and the other way round.#### Main effect 2: `POR_FINAL_SIB`Here's the effects object for the second main effect, `POR_FINAL_SIB`:```{r}fs_d <-data.frame(fs <-effect("POR_FINAL_SIB", m_final))# we add the predicted log odds as well:fs_d$fit.lo <-qlogis(fs_d$fit)# for later, we add the proportions of values in the data:fs_d$PROP <- d$POR_FINAL_SIB %>% table %>% prop.table %>% as.numericfs_d```Let's visualize:```{r}#| label: fig-FS#| fig-cap: The effect of final sibilants of the possessor# plot(fs, grid=TRUE)plot(fs, type="response", ylim=c(0, 1), grid=TRUE)```I would probably plot that effect like this:```{r}#| label: fig-FS_better#| fig-cap: The effect of final sibilants of the possessorop <-par(mar=c(4, 4, 4, 4)) # make the plotting margins big enoughplot(pch=16, col="#00000060", axes=FALSE,cex=4* fs_d$PROP,xlab="Sibilant at end of possessor" , x=1:2,ylab="Pred. prob. of s-genitives", ylim=c(0, 1), y=fs_d$fit)abline(h=seq(0.1, 0.9, 0.2), lty=3, col="#00000030")axis(1, at=1:2, labels=levels(d$POR_FINAL_SIB))axis(2, at=seq(0.1, 0.9, 0.2)); abline(h=0.5, lty=3)axis(4, at=seq(0.1, 0.9, 0.2), labels=round(qlogis(seq(0.1, 0.9, 0.2)), 2))mtext("Pred. log odds", 4, 2.5)for (i inseq(nrow(fs_d))) {arrows(i, fs_d$lower[i], i, fs_d$upper[i], angle=90, code=3)}par(op)```Interpretation: a clear articulatory effect: when the possessor ends in a sibilant, *s*-genitives are less likely.#### Interaction 1: `SPEAKER:POR_ANIMACY`Here's the effects object for the first interaction, `SPEAKER:POR_ANIMACY`:```{r}span_d <-data.frame(span <-effect("SPEAKER:POR_ANIMACY", m_final))# we add the predicted log odds as well:span_d$fit.lo <-qlogis(span_d$fit)# for later, we add the proportions of values in the data:span_d$PROP <- d$SPEAKER %>%table(d$POR_ANIMACY) %>% prop.table %>% as.numericspan_d```Let's visualize:```{r}#| label: fig-SPAN#| fig-cap: The effect of Speaker interacting with Possessor animacy#| fig-show: holdplot(span, type="response", ylim=c(0, 1), grid=TRUE,multiline=TRUE, confint=list(style="auto"))# the next plot is less useful, I think, but you should still have plotted it:# plot(span, type="response", ylim=c(0, 1), grid=TRUE,# x.var="SPEAKER", multiline=TRUE, confint=list(style="auto"))```I would probably plot that effect like this:```{r}#| label: fig-SPAN_better#| fig-cap: The effect of Speaker interacting with Possessor animacyx.coords <-rep(1:5, each=2)+c(-0.05, 0.05)y.cols <-rainbow(2, alpha=0.5)[span_d$SPEAKER]op <-par(mar=c(4, 4, 4, 4)) # make the plotting margins big enoughplot(pch=16, axes=FALSE,col=y.cols,cex=0.5+(5*span_d$PROP),xlab="Possessor animacy" , xlim=c(0.9, 5.1), x=x.coords,ylab="Pred. prob. of s-genitives", ylim=c(0 , 1) , y=span_d$fit)abline(h=seq(0.1, 0.9, 0.2), lty=3, col="#00000030")axis(1, at=1:5, labels=levels(d$POR_ANIMACY))axis(2, at=seq(0.1, 0.9, 0.2)); abline(h=0.5, lty=2)axis(4, at=seq(0.1, 0.9, 0.2), labels=round(qlogis(seq(0.1, 0.9, 0.2)), 2))mtext("Pred. log odds", 4, 2.5)for (i inseq(nrow(span_d))) {arrows(x.coords[i], span_d$lower[i], x.coords[i], span_d$upper[i],angle=90, code=3, col=y.cols[i])}legend(title="Speaker", bty="n",x=3, xjust=0.5, y=0.8, yjust=0.5, ncol=2,legend=levels(d$SPEAKER), fill=unique(y.cols))par(op)```Interpretation: Both NS and NNS are most likely to use *s*-genitives with animate possessors, followed by collective/temporal possessors, then locative, and finally inanimate ones (e.g., concrete objects), but the only significant difference in the interaction is that NS are significantly more likely to use *s*-genitives with collectives than NNS.#### Interaction 2: `SPEAKER:MODALITY`Here's the effects object for the second interaction, `SPEAKER:MODALITY`:```{r}spmo_d <-data.frame(spmo <-effect("SPEAKER:MODALITY", m_final))# we add the predicted log odds as well:spmo_d$fit.lo <-qlogis(spmo_d$fit)# for later, we add the proportions of values in the data:spmo_d$PROP <- d$SPEAKER %>%table(d$MODALITY) %>% prop.table %>% as.numericspmo_d```Let's visualize:```{r}#| label: fig-SPMO#| fig-cap: The effect of Speaker interacting with Modality#| fig-show: holdplot(spmo, type="response", ylim=c(0, 1), grid=TRUE,multiline=TRUE, confint=list(style="auto"))# the next plot is not any more useful, but you should still have plotted it:# plot(spmo, type="response", ylim=c(0, 1), grid=TRUE,# x.var="MODALITY", multiline=TRUE, confint=list(style="auto"))```I would probably plot that effect like this:```{r}#| label: fig-SPMO_better#| fig-cap: The effect of Speaker interacting with Modalityx.coords <-rep(1:2, each=2)+c(-0.05, 0.05)y.cols <-rainbow(2, alpha=0.5)[span_d$SPEAKER]op <-par(mar=c(4, 4, 4, 4)) # make the plotting margins big enoughplot(pch=16, axes=FALSE,cex=0.5+(5*spmo_d$PROP),col=y.cols,xlab="Modality" , xlim=c(0.9, 2.1), x=x.coords,ylab="Pred. prob. of s-genitives", ylim=c(0 , 1) , y=spmo_d$fit)abline(h=seq(0.1, 0.9, 0.2), lty=3, col="#00000030")axis(1, at=1:2, labels=levels(d$MODALITY))axis(2, at=seq(0.1, 0.9, 0.2)); abline(h=0.5, lty=2)axis(4, at=seq(0.1, 0.9, 0.2), labels=round(qlogis(seq(0.1, 0.9, 0.2)), 2))mtext("Pred. log odds", 4, 2.5)for (i inseq(nrow(spmo_d))) {arrows(x.coords[i], spmo_d$lower[i], x.coords[i], spmo_d$upper[i],angle=90, code=3, col=y.cols[i])}legend(title="Speaker", bty="n",x=1.5, xjust=0.5, y=0.8, yjust=0.5, ncol=2,legend=levels(d$SPEAKER), fill=unique(y.cols))par(op)```Interpretation: NS are more likely to use *s*-genitives in writing than they are in speaking and than NNS are in writing. NNS are less likely to use *s*-genitives in writing than they are in speaking and than NS are in writing. The highest predicted probability of *s*-genitives is found for NNS in speaking and the lowest predicted probability of *s*-genitives is found for NS in speaking.### Write-up {#sec-s01_binlog_writeup}To determine whether the choice of a genitive varies as a function of* whether the speaker is a native speaker or not;* the length difference between possessor and possessum and whether native and non-native speakers react to it in the same way;* the animacy of the possessor and whether native and non-native speakers react to it in the same way;* whether the possessor ends in a sibilant or not and whether native and non-native speakers react to it in the same way;* whether the modality is speaking or writing and whether native and non-native speakers react to it in the same way;* the interaction of `POR_FINAL_SIB` with `MODALITY`.a generalized linear model selection process was undertaken (using backwards stepwise model selection of all predictors and controls based on significance testing). The final model resulting from the elimination of ns predictors contained the effects of `LEN_PORmPUM_LOG:POR_ANIMACY` and `SPEAKER:POR_ANIMACY` and was highly significant (*LR*=1596.2, *df*=14, *p*<0.0001) with a good amount of explanatory/predictive power (McFadden's *R*^2^=0.397, Nagelkerke's *R*^2^=0.534, *C*=0.899). The summary table of the model is shown here.| | Estimate | 95%-CI | *se* | *z* | *p*~2-tailed~ ||:---------------------------------------------------------------------|:---------|:---------------|:------|:-------|:--------------|| (Intercept (LEN_PORmPUM_LOG=0, POR_ANIMACY=*animate*, SPEAKER=*nns*) | 1.22 | [0.99, 1.47] | 0.12 | 10.04 | <0.001 || SPEAKER~*nns*→*ns*~ | -0.54 | [-0.95, -0.12] | 0.21 | -2.53 | 0.011 || LEN_PORmPUM_LOG~0→1~ | -0.74 | [-0.84, -0.65] | 0.05 | -15.21 | <0.001 || POR_ANIMACY~*animate*→*collective*~ | -1.92 | [-2.23, -1.62] | 0.15 | -12.43 | <0.001 || POR_ANIMACY~*animate*→*inanimate*~ | -4.49 | [-4.95, -4.07] | 0.22 | -20.19 | <0.001 || POR_ANIMACY~*animate*→*locative*~ | -2.49 | [-2.96, -2.05] | 0.23 | -10.73 | <0.001 || POR_ANIMACY~*animate*→*temporal*~ | -1.65 | [-2.14, -1.19] | 0.24 | -6.83 | <0.001 || POR_FINAL_SIB~*absent*→*present*~ | -0.99 | [-1.26, -0.72] | 0.14 | -7.2 | <0.001 || MODALITY~*spoken*→*written*~ | -0.46 | [-0.71, -0.22] | 0.12 | -3.76 | 0.0002 || SPEAKER~*nns*→*ns*~:POR_ANIMACY~*animate*→*collective*~ | 0.87 | [0.33, 1.42] | 0.28 | 3.14 | 0.0017 || SPEAKER~*nns*→*ns*~:POR_ANIMACY~*animate*→*inanimate*~ | -0.27 | [-1.32, 0.63] | 0.49 | -0.55 | 0.581 || SPEAKER~*nns*→*ns*~:POR_ANIMACY~*animate*→*locative*~ | -0.5 | [-1.39, 0.34] | 0.44 | -1.13 | 0.286 || SPEAKER~*nns*→*ns*~:POR_ANIMACY~*animate*→*temporal*~ | 0.26 | [-0.59, 1.09] | 0.43 | 0.6 | 0.546 || SPEAKER~*nns*→*ns*~:MODALITY~*spoken*→*written*~ | 0.86 | [0.4, 1.33] | 0.24 | 3.65 | 0.0003 |The results indicate that* there is a clear main effect of `LEN_PORmPUM_LOG` amounting to a short-before-long effect;* there is a clear main effect of `POR_FINAL_SIB`: when the possessor ends in a sibilant, *s*-genitives are less likely;* the interaction `SPEAKER:POR_ANIMACY` only has the effect that native speakers are significantly more likely to use *s* genitives with collectives than learners are -- all other speaker differences within each animacy level are ns. [Show plot(s)]* the interaction `SPEAKER:MODALITY` has the effect that NS are more likely to use *s*-genitives in writing than they are in speaking and than NNS are in writing. NNS are less likely to use *s*-genitives in writing than they are in speaking and than NS are in writing. [Show plot(s)]# Session info```{r}sessionInfo()```
