It is often the case that a factor only makes sense in a subset of a dataset (i.e. for some observations a factor may simply not be meaningful or contrastive), or that with observational datasets there are no observations in some levels of an interaction term. There are also cases where a random effects grouping factor is only applicable in a subset the data, and it is desireable to model the noise introduced by the repeated measures on the group members within the subset of the data where the repeated measures exist. The nauf package allows unordered factors and random effects grouping factors to be coded as NA in the subsets of the data where they are not applicable or otherwise not contrastive. Sum contrasts are used for all unordered factors (using named_contr_sum in the standardize package), and then NA values are set to 0. This allows all of the data to be modeled together without creating collinearity or making the output difficult to interpret. It is highly recommended that regression variables be put on the same scale with the standardize function in the standardize package prior to using nauf functions (though this is not required for the functions to work). The “Using the standardize package” vignette also provides useful information on the difference between unordered and ordered factors.
The plosives dataset (Eager, 2017) included in the nauf package contains measures of plosive strength for instances of intervocalic Spanish /p/, /t/, /k/, /b/, /d/ and /g/ from speakers from three dialects. For this first example we will focus on the voiceless plosives /ptk/ from the 30 speakers in the Cuzco dialect. The dataset contains a combination of experimental data (with exactly repeated measures on items from a read speech task) and observational data (without exactly repeated measures; from interviews with the same speakers). For the spontaneous speech (as indicated by the logical variable spont) the item factor is coded as NA.
library(nauf)
summary(plosives)
#> cdur vdur vpct intdiff
#> Min. : 0.0 Min. : 0.00 Min. :0.0000 Min. : 0.000
#> 1st Qu.: 93.0 1st Qu.: 0.00 1st Qu.:0.0000 1st Qu.: 5.951
#> Median :127.0 Median : 0.00 Median :0.0000 Median :22.452
#> Mean :123.9 Mean : 39.21 Mean :0.2437 Mean :22.201
#> 3rd Qu.:161.0 3rd Qu.: 81.00 3rd Qu.:0.5213 3rd Qu.:36.277
#> Max. :298.0 Max. :212.00 Max. :0.9079 Max. :73.952
#>
#> intvel voicing place stress
#> Min. :0.0000 Voiced :2694 Bilabial:1752 Post-Tonic:1629
#> 1st Qu.:0.1974 Voiceless:2587 Dental :1862 Tonic :2007
#> Median :0.7440 Velar :1667 Unstressed:1645
#> Mean :0.8341
#> 3rd Qu.:1.2870
#> Max. :5.5594
#>
#> prevowel posvowel wordpos wordfreq speechrate
#> a:1675 a:1942 Initial:1166 Min. : 1 Min. : 1.000
#> e:1237 e: 900 Medial :4115 1st Qu.: 1057 1st Qu.: 4.115
#> i:1280 i: 782 Median : 4579 Median : 5.000
#> o: 708 o:1217 Mean : 20189 Mean : 5.236
#> u: 381 u: 440 3rd Qu.: 20368 3rd Qu.: 6.000
#> Max. :247340 Max. :10.000
#>
#> spont dialect sex age
#> Mode :logical Cuzco :2872 Female:2599 Older :1362
#> FALSE:2019 Lima : 862 Male :2682 Younger:3919
#> TRUE :3262 Valladolid:1547
#> NA's :0
#>
#>
#>
#> ed ling speaker item
#> Secondary :1445 Bilingual :1379 s34 : 133 i01 : 38
#> University:3836 Monolingual:3902 s09 : 124 i02 : 38
#> s32 : 124 i03 : 38
#> s22 : 120 i04 : 38
#> s38 : 116 i05 : 38
#> s25 : 114 (Other):1829
#> (Other):4550 NA's :3262
dat <- droplevels(subset(plosives, dialect == "Cuzco" & voicing == "Voiceless"))
xtabs(~ is.na(item) + spont, dat)
#> spont
#> is.na(item) FALSE TRUE
#> FALSE 795 0
#> TRUE 0 620For our example, we want to model the voiceless duration of the plosives as a function of place of articulation, stress, and whether or not the speech was spontaneous. When modeling the read speech data (spont = FALSE), we want to account for noise introduced by the repeated measures on item, but we can’t apply this random effects structure to the interview data (spont = TRUE). In addition to this, we want to account for the noise introduced by the individual speakers. Rather than leaving out the item effects to analyze all the data together, or keeping the item effects and analyzing the read and spontaneous speech separately, we can model both subsets together and have the item effects apply only within the read speech data using nauf. We just need to have item coded as NA when it is not applicable, which it already is. Then the nauf_lmer function takes care of the rest.
sobj <- standardize(vdur ~ spont + place + stress +
(1 + spont + place + stress | speaker) + (1 | item), dat)
mod <- nauf_lmer(sobj$formula, sobj$data)
summary(mod)
#> Linear mixed model fit by REML ['nauf.lmerMod']
#> Formula: vdur ~ spont + place + stress + (1 + spont + place + stress |
#> speaker) + (1 | item)
#> Data: sobj$data
#>
#> REML criterion at convergence: 3543
#>
#> Scaled residuals:
#> Min 1Q Median 3Q Max
#> -5.1474 -0.5000 0.0235 0.5397 4.6395
#>
#> Random effects:
#> Groups Name Variance Std.Dev. Corr
#> speaker (Intercept) 0.215923 0.46468
#> spontTRUE 0.002052 0.04530 -0.09
#> placeBilabial 0.016284 0.12761 -0.44 0.49
#> placeDental 0.008709 0.09332 0.01 -0.93 -0.14
#> stressPost-Tonic 0.024526 0.15661 0.35 -0.92 -0.78 0.73
#> stressTonic 0.020310 0.14251 -0.41 0.51 1.00 -0.16 -0.79
#> item (Intercept) 0.109717 0.33123
#> Residual 0.621677 0.78846
#> Number of obs: 1415, groups: speaker, 30; item, 27
#>
#> Fixed effects:
#> Estimate Std. Error t value
#> (Intercept) -0.011606 0.093282 -0.124
#> spontTRUE -0.182422 0.039642 -4.602
#> placeBilabial -0.014944 0.048795 -0.306
#> placeDental 0.001387 0.044381 0.031
#> stressPost-Tonic 0.098136 0.054977 1.785
#> stressTonic 0.143076 0.047996 2.981
#>
#> Correlation of Fixed Effects:
#> (Intr) spTRUE plcBlb plcDnt strP-T
#> spontTRUE -0.266
#> placeBilabl -0.176 0.080
#> placeDental -0.008 -0.102 -0.450
#> strssPst-Tn 0.206 -0.005 -0.162 0.058
#> stressTonic -0.235 -0.021 0.292 -0.021 -0.613This way, we are making use of all of the information in the dataset. We can obtain a more principled statistical test for the spont factor, and also get better information about the other fixed effects and the individual speakers, since the same random effects for speaker apply in both the read speech and spontaneous subsets. We can obtain Type III tests using anova (here with method = “S” to indicate Satterthwaite approximation for the denominator degrees of freedom in the F tests).
anova(mod, method = "S")
#> Mixed Model Anova Table (Type III tests, S-method)
#>
#> Model: vdur ~ spont + place + stress + (1 + spont + place + stress |
#> Model: speaker) + (1 | item)
#> Data: $
#> Data: sobj
#> Data: data
#> num Df den Df F Pr(>F)
#> spont 1 36.900 21.176 4.817e-05 ***
#> place 2 80.249 0.054 0.9475
#> stress 2 74.356 14.893 3.635e-06 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1We see that stress is significant, and can then get predicted marginal means (often called least-squares means) and pairwise comparisons for its levels using nauf_ref.grid to create a reference grid, and then calling nauf_pmmeans.
rg <- nauf_ref.grid(mod)
nauf_pmmeans(rg, "stress", pairwise = TRUE)
#>
#> Predicted marginal means for 'stress'
#>
#> Factors averaged over: 'spont' 'place'
#>
#> $pmmeans
#> stress pmmean SE df lower.CL upper.CL
#> Post-Tonic 0.08653011 0.11764046 42.24 -0.15083790 0.32389813
#> Tonic 0.13146972 0.09434706 41.10 -0.05905417 0.32199361
#> Unstressed -0.25281748 0.10383491 41.03 -0.46251171 -0.04312325
#>
#> Confidence level used: 0.95
#>
#> $contrasts
#> contrast estimate SE df t.ratio p.value
#> Post-Tonic - Tonic -0.04493961 0.09252431 59.26 -0.486 0.8783
#> Post-Tonic - Unstressed 0.33934759 0.08901834 87.70 3.812 0.0007
#> Tonic - Unstressed 0.38428720 0.07594482 82.73 5.060 <.0001
#>
#> P value adjustment: tukey method for comparing a family of 3 estimatesThe fricatives dataset (Hualde and Prieto, 2014) included in the nauf package contains measures of duration and voicing for intervocalic alveolar fricatives in Spanish and Catalan. Spanish has only one such fricative, /s/ (underlyingly voiceless), which can occur in any position in the word (initial, medial, or final). In Catalan, the situation is much more complicated. At the beginning of a word and in the middle of a word, both /s/ (underlyingly voiceless) and /z/ (underlyingly voiced) can occur (though word-initial /z/ is rare, and does not occur in the dataset). In word-final position, there is no contrast between /s/ and /z/ (labeled /S/, underlyingly neutralized), with the voicing of the fricative determined by the following sound. Because all of the fricatives in the dataset are intervocalic, /S/ ought to be like /z/ (according to traditional descriptions of Catalan), but may be shorter and possibly less voiced. That is, in the dataset, we have the following set of unique values.
summary(fricatives)
#> dur pvoi lang wordpos
#> Min. : 29.71 Min. :0.0625 Catalan:974 Final :296
#> 1st Qu.: 71.14 1st Qu.:0.5263 Spanish:648 Initial:397
#> Median : 86.29 Median :0.6471 Medial :929
#> Mean : 90.35 Mean :0.6937
#> 3rd Qu.:105.53 3rd Qu.:0.9231
#> Max. :250.78 Max. :1.0000
#>
#> uvoi speaker
#> Neutralized: 160 s12 : 142
#> Voiced : 327 s09 : 141
#> Voiceless :1135 s10 : 102
#> s20 : 71
#> s36 : 70
#> s19 : 64
#> (Other):1032
dat <- fricatives
u <- unique(dat[, c("lang", "wordpos", "uvoi")])
u <- u[order(u$lang, u$wordpos, u$uvoi), ]
rownames(u) <- NULL
u
#> lang wordpos uvoi
#> 1 Catalan Final Neutralized
#> 2 Catalan Initial Voiceless
#> 3 Catalan Medial Voiced
#> 4 Catalan Medial Voiceless
#> 5 Spanish Final Voiceless
#> 6 Spanish Initial Voiceless
#> 7 Spanish Medial VoicelessThis raises a problem for the regression analysis, since underlying voicing uvoi is only contrastive within a subset of the data (specifically, when lang = “Catalan” and wordpos = “Medial”), and everywhere else is uniquely determined by lang and wordpos. Ideally, we would like to be able to have uvoi apply only where it is contrastive, and to be able to include random slopes for uvoi for the Catalan speakers, but not the Spanish speakers. Using traditional methods won’t help us here. If we were to take the full interaction of the three factors as a new factor with 7 levels, we can’t include slopes. If we run separate regressions on different subsets of the data, then we will have to limit the comparisons we can make, and won’t be making use of all of the information the data is providing us. Note also that this issue has nothing to do with the way the data were collected; the nature of the languages and the research questions creates these imbalances. With nauf, we can solve this problem by coding uvoi as NA in the subsets of the data where it is not contrastive. That is we can code the unique possible combinations as:
u$uvoi[!(u$lang == "Catalan" & u$wordpos == "Medial")] <- NA
u
#> lang wordpos uvoi
#> 1 Catalan Final <NA>
#> 2 Catalan Initial <NA>
#> 3 Catalan Medial Voiced
#> 4 Catalan Medial Voiceless
#> 5 Spanish Final <NA>
#> 6 Spanish Initial <NA>
#> 7 Spanish Medial <NA>
dat$uvoi[!(dat$lang == "Catalan" & dat$wordpos == "Medial")] <- NAIn this way, we are recognizing that word-final Catalan fricatives are limited to one value for uvoi, as are all Spanish fricatives, etc. The meaning of NA in this table is not always the same, so we have to keep track of it, and we will definitely want to include an interaction between lang and wordpos since, for example, “Catalan:Final:NA” means neutralized /S/ while “Spanish:Final:NA” means voiceless /s/. As for being able to include uvoi slopes for Catalan speakers and not Spanish speakers, we can create two new random effects grouping factors, one for each language, where the factor is the speaker identifier or NA based on the language the observation comes from:
dat$ca_speaker <- dat$sp_speaker <- dat$speaker
dat$ca_speaker[dat$lang != "Catalan"] <- NA
dat$sp_speaker[dat$lang != "Spanish"] <- NAWith these NA values, we can then fit a model with nauf_lmer to predict the duration of the fricatives based on lang, wordpos, and uvoi:
sobj <- standardize(dur ~ lang * wordpos + uvoi +
(1 + wordpos + uvoi | ca_speaker) + (1 + wordpos | sp_speaker), dat)
mod <- nauf_lmer(sobj$formula, sobj$data)
summary(mod)
#> Linear mixed model fit by REML ['nauf.lmerMod']
#> Formula:
#> dur ~ lang * wordpos + uvoi + (1 + wordpos + uvoi | ca_speaker) +
#> (1 + wordpos | sp_speaker)
#> Data: sobj$data
#>
#> REML criterion at convergence: 4003.3
#>
#> Scaled residuals:
#> Min 1Q Median 3Q Max
#> -3.6998 -0.6383 -0.0850 0.5507 6.6149
#>
#> Random effects:
#> Groups Name Variance Std.Dev. Corr
#> ca_speaker (Intercept) 0.057124 0.2390
#> wordposFinal 0.076718 0.2770 -0.75
#> wordposInitial 0.040187 0.2005 0.59 -0.98
#> uvoiVoiced 0.013167 0.1147 -1.00 0.69 -0.52
#> sp_speaker (Intercept) 0.068574 0.2619
#> wordposFinal 0.004475 0.0669 0.87
#> wordposInitial 0.021641 0.1471 -0.72 -0.97
#> Residual 0.645343 0.8033
#> Number of obs: 1622, groups: ca_speaker, 26; sp_speaker, 16
#>
#> Fixed effects:
#> Estimate Std. Error t value
#> (Intercept) -0.007564 0.048331 -0.157
#> langCatalan 0.029407 0.048331 0.608
#> wordposFinal -0.505660 0.048465 -10.433
#> wordposInitial 0.408757 0.045600 8.964
#> uvoiVoiced -0.514556 0.042640 -12.067
#> langCatalan:wordposFinal -0.231350 0.048465 -4.774
#> langCatalan:wordposInitial 0.338058 0.045600 7.414
#>
#> Correlation of Fixed Effects:
#> (Intr) lngCtl wrdpsF wrdpsI uvoVcd lngC:F
#> langCatalan -0.285
#> wordposFinl 0.028 -0.297
#> wordposIntl -0.067 0.290 -0.768
#> uvoiVoiced -0.292 -0.292 0.216 -0.083
#> lngCtln:wrF -0.297 0.028 0.289 -0.154 0.216
#> lngCtln:wrI 0.290 -0.067 -0.154 -0.077 -0.083 -0.768To understand how nauf treats the NA values in uvoi, we call nauf_contrasts:
nauf_contrasts(mod)
#> $lang
#> Catalan
#> Catalan 1
#> Spanish -1
#>
#> $wordpos
#> Final Initial
#> Final 1 0
#> Initial 0 1
#> Medial -1 -1
#>
#> $uvoi
#> Voiced
#> Voiced 1
#> Voiceless -1
#> <NA> 0The function returns a list with an element for each of the three factors. Sum contrasts have been applied for all three with named_contr_sum from the standardize package, but for uvoi there is an additional row indicating that a value of 0 is used when it is NA. In this way, the uvoiVoiced coefficient in the summary never contributes to the fitted value for any observation that belongs to a subset where underlying voicing is not contrastive, and its estimate represents half the difference between underlyingly voiced and voiceless Catalan word-medial observations. We can run an anova just as we did in the Cuzco example (though in this case it is not as useful since we are interested in making specific comparisons), and also create a reference grid:
anova(mod, method = "S")
#> Mixed Model Anova Table (Type III tests, S-method)
#>
#> Model: dur ~ lang * wordpos + uvoi + (1 + wordpos + uvoi | ca_speaker) +
#> Model: (1 + wordpos | sp_speaker)
#> Data: $
#> Data: sobj
#> Data: data
#> num Df den Df F Pr(>F)
#> lang 1 27.374 0.3702 0.5479
#> wordpos 2 26.862 55.5357 2.863e-10 ***
#> uvoi 1 44.259 145.6220 1.332e-15 ***
#> lang:wordpos 2 26.862 28.5074 2.282e-07 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
rg <- nauf_ref.grid(mod)We can now use the nauf_pmmeans function’s subset, by, and na_as_level arguments to test hypotheses in the subset of the reference grid where comparisons are valid. For example, if we want to test whether Spanish fricatives are longer than Catalan fricatives, the only just comparison we can make is for word-initial and word-medial /s/. To do this, we provide nauf_pmmeans with a list of valid groups, with each group specified as a named list of unordered factor levels:
nauf_pmmeans(rg, "lang", pairwise = TRUE,
subset = list(
list(lang = "Catalan", wordpos = "Initial", uvoi = NA),
list(lang = "Catalan", wordpos = "Medial", uvoi = "Voiceless"),
list(lang = "Spanish", wordpos = c("Initial", "Medial"), uvoi = NA)
)
)
#>
#> Predicted marginal means for 'lang'
#> Note: 'lang' is included in higher order interaction(s) 'lang:wordpos'
#>
#> Factors conditioned on: 'lang' 'wordpos' 'uvoi'
#>
#> See the 'subset' element of the 'specs' attribute for subsetted groups
#>
#> $pmmeans
#> lang pmmean SE df lower.CL upper.CL
#> Catalan 0.6476254 0.09069385 26.54 0.46138613 0.8338646
#> Spanish 0.1001844 0.07289153 13.23 -0.05700457 0.2573734
#>
#> Confidence level used: 0.95
#>
#> $contrasts
#> contrast estimate SE df t.ratio p.value
#> Catalan - Spanish 0.547441 0.1163553 39.15 4.705 <.0001Each of the lists specified in the subset list represents a set of valid combinations of factor levels where language is truly contrastive. Any row in the reference grid which belongs to at least one of the specified groups is kept and the others are ignored. So the Catalan estimate is an average of Catalan:Initial:NA and Catalan:Medial:Voiceless, and the Spanish estimate is an average of Spanish:Initial:NA and Spanish:Medial:NA. To test the effect of word position in the two languages separately, we can make use of the by argument, and specify the groups for each language where word position is truly contrastive:
nauf_pmmeans(rg, c("lang", "wordpos"), pairwise = TRUE, by = "lang",
subset = list(
list(lang = "Catalan", wordpos = "Initial", uvoi = NA),
list(lang = "Catalan", wordpos = "Medial", uvoi = "Voiceless"),
list(lang = "Spanish", uvoi = NA)
)
)
#>
#> Predicted marginal means for 'lang:wordpos'
#> Pairwise comparisons by 'lang'
#>
#> Factors conditioned on: 'lang' 'wordpos' 'uvoi'
#>
#> See the 'subset' element of the 'specs' attribute for subsetted groups
#>
#> $pmmeans
#> lang wordpos pmmean SE df lower.CL upper.CL
#> Catalan Initial 0.76865796 0.09560913 26.21 0.57220781 0.96510810
#> Catalan Medial 0.52659283 0.10499041 26.58 0.31100913 0.74217652
#> Spanish Final -0.31128082 0.11129830 13.47 -0.55087547 -0.07168617
#> Spanish Initial 0.03372864 0.08563586 10.03 -0.15700539 0.22446267
#> Spanish Medial 0.16664021 0.09218783 12.00 -0.03421486 0.36749527
#>
#> Confidence level used: 0.95
#>
#> $contrasts_1
#> contrast estimate SE df t.ratio
#> Catalan,Initial - Catalan,Medial 0.2420651 0.08617821 29.03 2.809
#> p.value
#> 0.0088
#>
#>
#> $contrasts_2
#> contrast estimate SE df t.ratio
#> Spanish,Final - Spanish,Initial -0.3450095 0.11508620 12.33 -2.998
#> Spanish,Final - Spanish,Medial -0.4779210 0.08357926 120.38 -5.718
#> Spanish,Initial - Spanish,Medial -0.1329116 0.10203672 5.96 -1.303
#> p.value
#> 0.0272
#> <.0001
#> 0.4444
#>
#> P value adjustment: tukey method for comparing a family of 3 estimatesIn this way we obtain 5 predicted marginal means, and pairwise comparisons within-language. Underlying voicing is only contrastive for Catalan medial fricatives, and so we can call:
nauf_pmmeans(rg, "uvoi", pairwise = TRUE,
subset = list(lang = "Catalan", wordpos = "Medial")
)
#>
#> Predicted marginal means for 'uvoi'
#> NA not considered a level for: 'uvoi'
#>
#> Factors conditioned on: 'lang' 'wordpos'
#>
#> See the 'subset' element of the 'specs' attribute for subsetted groups
#>
#> $pmmeans
#> uvoi pmmean SE df lower.CL upper.CL
#> Voiced -0.5025186 0.06253773 31.46 -0.6299897 -0.3750475
#> Voiceless 0.5265928 0.10499041 26.58 0.3110091 0.7421765
#>
#> Confidence level used: 0.95
#>
#> $contrasts
#> contrast estimate SE df t.ratio p.value
#> Voiced - Voiceless -1.029111 0.08528034 44.26 -12.067 <.0001In this call, we are only specifying one group in subset, so we don’t need to double-list it (i.e. nauf will understnd list(lang = “Catalan”, wordpos = “Medial”) as list(list(lang = “Catalan”, wordpos = “Medial”))). For Catalan word-final neutralized /S/, we may be interested in making a comparison to see whether it is somewhere in between /s/ and /z/ in terms of duration, or not significantly different from /z/. This case requires the additional argument na_as_level to be specified since we want an estimate for a group where uvoi is NA (the default is that estiamtes are not generated for any group that has NA’s in the estimate table):
nauf_pmmeans(rg, c("wordpos", "uvoi"), pairwise = TRUE, na_as_level = "uvoi",
subset = list(
list(lang = "Catalan", wordpos = "Medial", uvoi = c("Voiced", "Voiceless")),
list(lang = "Catalan", wordpos = "Final", uvoi = NA)
)
)
#>
#> Predicted marginal means for 'wordpos:uvoi'
#> NA considered a level for: 'uvoi'
#> Note: The interaction term 'wordpos:uvoi' is not in the model.
#>
#> Factors conditioned on: 'lang' 'wordpos' 'uvoi'
#>
#> See the 'subset' element of the 'specs' attribute for subsetted groups
#>
#> $pmmeans
#> wordpos uvoi pmmean SE df lower.CL upper.CL
#> Final NA -0.7151669 0.08290231 14.38 -0.8925304 -0.5378035
#> Medial Voiced -0.5025186 0.06253773 31.46 -0.6299897 -0.3750475
#> Medial Voiceless 0.5265928 0.10499041 26.58 0.3110091 0.7421765
#>
#> Confidence level used: 0.95
#>
#> $contrasts
#> contrast estimate SE df t.ratio
#> Final,NA - Medial,Voiced -0.2126484 0.10320044 19.49 -2.061
#> Final,NA - Medial,Voiceless -1.2417598 0.12851501 24.82 -9.662
#> Medial,Voiced - Medial,Voiceless -1.0291114 0.08528034 44.26 -12.067
#> p.value
#> 0.1245
#> <.0001
#> <.0001
#>
#> P value adjustment: tukey method for comparing a family of 3 estimatesAs a final example, we could test if Spanish word-final /s/ (being shorter than Spanish word-medial /s/ above) is as short as Catalan word-final /S/:
nauf_pmmeans(rg, "lang", pairwise = TRUE,
subset = list(wordpos = "Final", uvoi = NA)
)
#>
#> Predicted marginal means for 'lang'
#> Note: 'lang' is included in higher order interaction(s) 'lang:wordpos'
#>
#> Factors conditioned on: 'wordpos' 'uvoi'
#>
#> See the 'subset' element of the 'specs' attribute for subsetted groups
#>
#> $pmmeans
#> lang pmmean SE df lower.CL upper.CL
#> Catalan -0.7151669 0.08290231 14.38 -0.8925304 -0.53780347
#> Spanish -0.3112808 0.11129830 13.47 -0.5508755 -0.07168617
#>
#> Confidence level used: 0.95
#>
#> $contrasts
#> contrast estimate SE df t.ratio p.value
#> Catalan - Spanish -0.4038861 0.1387808 25.26 -2.91 0.0074In this way, we are able to use all of the information in the dataset when fitting the model, account for the uncertainty introduced by repeated measures on the subjects, assigning different random effects to the two langauges, and then test any variety of hypotheses we might have about the effects of the factors within the subsets of the data where the comparisons make sense.
There are two situations in which unordered factors will need more than one set of contrasts: (1) when an unordered factor with NA values interacts with another unordered factor, and some levels are collinear with NA; and (2) when an unordered factor is included as a slope for a random effects grouping factor that has NA values, but only a subset of the levels for the slope factor occur when the grouping factor is applicable. Both of these situations occur when we consider all three dialects in the plosives dataset jointly (here we will look at the voiceless subset):
dat <- subset(plosives, voicing == "Voiceless")
xtabs(~ dialect + spont, dat)
#> spont
#> dialect FALSE TRUE
#> Cuzco 795 620
#> Lima 215 206
#> Valladolid 0 751The data for Cuzco and Lima consist of two subsets: read speech and spontaneous speech (with the read speech task being the same for both dialects). The Valladolid data, however, come from a different corpus consisting only of spontaneous speech, with the measurements taken in the same way as for Cuzco and Lima. While it would of course be ideal to have read speech for the Valladolid speakers as well, this doesn’t mean that we need to split up the data and run multiple regressions. We can simply code spont as NA for Valladolid, split the speaker random effects by dialect, and only include a spont slope for Cuzco and Lima:
dat$spont[dat$dialect == "Valladolid"] <- NA
dat$c_speaker <- dat$l_speaker <- dat$v_speaker <- dat$speaker
dat$c_speaker[dat$dialect != "Cuzco"] <- NA
dat$l_speaker[dat$dialect != "Lima"] <- NA
dat$v_speaker[dat$dialect != "Valladolid"] <- NA
sobj <- standardize(cdur ~ spont * dialect +
(1 + spont | c_speaker) + (1 + spont | l_speaker) + (1 | v_speaker) +
(1 + dialect | item),
dat)
mod <- nauf_lmer(sobj$formula, sobj$data)
summary(mod)
#> Linear mixed model fit by REML ['nauf.lmerMod']
#> Formula: cdur ~ spont * dialect + (1 + spont | c_speaker) + (1 + spont |
#> l_speaker) + (1 | v_speaker) + (1 + dialect | item)
#> Data: sobj$data
#>
#> REML criterion at convergence: 6242.3
#>
#> Scaled residuals:
#> Min 1Q Median 3Q Max
#> -3.0413 -0.6728 -0.1214 0.5382 4.5849
#>
#> Random effects:
#> Groups Name Variance Std.Dev. Corr
#> c_speaker (Intercept) 0.0846787 0.29100
#> spontTRUE 0.0065212 0.08075 -0.53
#> item (Intercept) 0.2805902 0.52971
#> dialect.c2.Cuzco 0.0009203 0.03034 -1.00
#> v_speaker (Intercept) 0.0633240 0.25164
#> l_speaker (Intercept) 0.0428048 0.20689
#> spontTRUE 0.0158277 0.12581 -1.00
#> Residual 0.6030539 0.77657
#> Number of obs: 2587, groups:
#> c_speaker, 30; item, 27; v_speaker, 18; l_speaker, 8
#>
#> Fixed effects:
#> Estimate Std. Error t value
#> (Intercept) -0.180515 0.052463 -3.441
#> spontTRUE -0.152609 0.060147 -2.537
#> dialectCuzco 0.574724 0.053716 10.699
#> dialectLima -0.203724 0.065241 -3.123
#> spontTRUE:dialect.c2.Cuzco -0.004251 0.032064 -0.133
#>
#> Correlation of Fixed Effects:
#> (Intr) spTRUE dlctCz dlctLm
#> spontTRUE -0.738
#> dialectCuzc -0.011 -0.092
#> dialectLima 0.392 -0.517 -0.421
#> spTRUE:.2.C 0.343 -0.434 -0.356 0.569In this way, speaker effects are accounted for in each dialect, and within the spont = FALSE subset, item effects are additionally accounted for. However, note that for the interaction term spont:dialect, .c2. appears before Cuzco, and there is only one coefficient for the interaction term rather than two as we would normally expect. The same goes for the item slope for dialect. This is because using the main effect contrasts for dialect in the spont:dialect interaction term and item dialect slope would lead to collinearity issues (as explained in detail in the help page for the nauf_contrasts function). The nauf package automatically recognizes when this happens, and creates additional sets of contrasts which it uses only when it needs to. To see how the contrasts are coded, we can call nauf_contrasts:
nauf_contrasts(mod)
#> $spont
#> TRUE
#> TRUE 1
#> FALSE -1
#> <NA> 0
#>
#> $dialect
#> Cuzco Lima .c2.Cuzco
#> Cuzco 1 0 1
#> Lima 0 1 -1
#> Valladolid -1 -1 0There are nauf regression functions for fixed effects regressions which would normally be fit with lm, glm, and glm.nb (from the stats and MASS packages), mixed effects regressions which would normally be fit with lmer, glmer, and glmer.nb (from the lme4 package), and Bayesian versions of these six regression functions which would normally be fit with stan_lm, stan_glm, stan_glm.nb, stan_lmer, stan_glmer, and stan_glmer.nb (from the rstanarm package). In each case, the nauf regression function has the same name but preceded by nauf_ (e.g. nauf_lm, nauf_glmer, nauf_stan_glmer.nb, etc.). For Bayesian regressions, there are methods for the generic functions in the rstantools package, as well as functions nauf_kfold and nauf_launch_shinystan, and nauf_ref.grid and nauf_pmmeans work in the same way as illustrated above, but returning summaries based on computing posterior marginal means at each iteration of the model.
The nauf package allows unordered factors and random effects grouping factors to be used even when they are only applicable within a subset of a dataset. The user only needs to code them as NA when they are not applicable/contrastive, and nauf regression fitting functions take care of the rest. Different random effects structures for the same grouping variable can be fit in different subsets by creating new factors from the grouping factor, and setting them to NA in the appropriate subsets. The nauf_pmmeans function can then be used to test hypotheses conditioning on subsets of the data where a just comparison can be made.
If you use the nauf package in a publication, please cite:
Eager, Christopher D. (2017). nauf: Regression with NA Values in Unordered Factors. R package version 1.1.0. https://CRAN.R-project.org/package=naufIf you analyze the plosives dataset in a publication, please cite:
Eager, Christopher D. (2017). Contrast preservation and constraints on individual phonetic variation. Doctoral thesis. University of Illinois at Urbana-Champaign.If you analyze the fricatives dataset in a publication, please cite:
Hualde, J. I., & Prieto, P. (2014). Lenition of intervocalic alveolar fricatives in Catalan and Spanish. Phonetica, 71(2), 109-127.