2
$\begingroup$

I am doing a linguistics project in which I am comparing the % of the time "um" follows a contraction vs. the % of time "um" follows what I call the full form (e.g. full form of don't is do not). I have 100 contractions I'm looking at. For each contraction, I treat "is followed by um or not" as a binomial random variable and run a z-test accordingly. I would like to then run a test on all these tests to test if contractions do have an effect on "um." The reason I can't just sum them all up and do one test is some contractions are way more frequent and therefore would dominate the data set. One issue is that it's possible the different contractions are not all distributed the same, such that the effect may only occur in some of them.

In general it seems like you should be able to interpret multiple hypothesis tests to see the significance of your testing. For instance, if I'm using an alpha of 0.05 and I only get a p-value < 0.05 in 5/100 of my tests, I could safely conclude those were just false positives and there isn't a general effect. Is there a more rigorous way to perform this test of tests?

$\endgroup$
  • $\begingroup$ Given "One issue is that it's possible the different contractions are not all distributed the same", why need to summarize them together? $\endgroup$ – user158565 Dec 12 '18 at 2:43
  • $\begingroup$ Because I am looking for a general trend. They might be distributed differently, they might not. I'd like to see if, in general, there is an effect. $\endgroup$ – jdude97 Dec 12 '18 at 4:23
2
$\begingroup$

Your "test of tests" sounds like a mixed model. In this formulation you aren't concerned with estimating the effect word by word, but rather want to know whether there's an overall effect, so you treat the variations in underlying rates by word as random.

My simulation experiments suggest you have a good chance to detect and quantify relatively subtle effects (such as a log odds ratio of $0.25$ or less) provided you have around 20 or more instances of each basic word on average. The simulation (detailed in code below) recreates all the features you have described: around 100 words, each possibly appearing in uncontracted or contracted form, with greatly varying frequencies.

Let's look at one simulation and what this model does with it. The data are represented in a table showing the base word form (a randomly created nonsense word), the form of the word that was observed, the number of times this form was followed by a pause ("um"), and the number times it was not followed by a pause.

> print(df[, .(base, form, contracted, pause, pause.not)])
              base          form contracted pause pause.not
  1:        abhbtf          abhb       TRUE     9        26
  2:        abhbtf        abhbtf      FALSE     1        40
  3:    agmkkgtrkv          agmk       TRUE     0        20
 ---                                                       
192: zprdpoobkstwu zprdpoobkstwu      FALSE     0         1
193:     zxpvvluon          zxpv       TRUE     5        35
194:     zxpvvluon     zxpvvluon      FALSE     4        34

Here is a plot of these data. Each dot reflects a single row of the table, located according to whether it is in contracted form (horizontally) and the observed frequency (vertically). Dots are sized in proportion to the total observed frequencies; for instance, the dot for the first data row (word abhb) has an area proportional to $9+26.$ The colored polygons are "violin plot" estimates of the data distributions (although they do not account for the variation in observation sizes).

Figure 1

There is a slight hint that pauses after contracted word forms (right violin) may be a little more frequent than pauses after base word forms (left violin). It's not perfectly clear that there is any real difference, though, because the many infrequent words (shown with tiny dots) have extremely variable frequencies, as high as 100% in several cases.

The mixed model supposes the true rate of response varies independently and randomly among all forms of all words (contracted or not), but that the mean rate for contractions differs from the mean rate for base forms of the words. This is a parsimonious model: it estimates just three parameters describing the rates for uncontracted word forms, the difference in mean log odds of rates (which is the parameter you wish to test), and the variance of those log odds.

Here is a typical model summary:

Random effects:
 Groups          Name        Variance Std.Dev.
 base:contracted (Intercept) 0.3397   0.5829  

Fixed effects:
               Estimate Std. Error z value Pr(>|z|)    
(Intercept)     -2.6114     0.1332 -19.607  < 2e-16 ***
contractedTRUE   0.4949     0.1724   2.871  0.00409 **

The value of $0.4949$ for contractedTRUE is the estimate of how much the "um" rate for contracted word forms exceeds the rate for uncontracted forms, expressed as a log odds ratio. Its small p-value of $0.00409$ would usually be considered significant evidence of a real underlying difference. Indeed, these data were simulated from a "reality" in which the base log odds is $-2.5$ and the difference is $+0.5.$ In this case, at least, the estimates are highly accurate.

The standard deviation of the random effect, $0.5829,$ indicates how variable the estimated rates are from word to word. The data were drawn from a population in which the true standard deviation is $0.5:$ the estimate is pretty good.

It is of interest to compare model-predicted rates with observed rates. As that SD of $0.58$ indicates, much of the variation in observed rates looks random. The mixed model "discounts" this inherent variation, essentially by "shrinking" the observed rates towards the overall average rate. Furthermore, words seen infrequently naturally experience more shrinkage towards the average. In the plot, this shows up in the small dots strung out to the right: they all had unusually high rates, but the model attributes that almost entirely to chance and thereby predicts their underlying rates to be near the average. The larger the dot, the more times that word was seen, and the closer the prediction comes to reproducing the observed frequencies.

Figure 2

Much more can be said. A good resource for learning about mixed models is Fitzmaurice et al., Applied Longitudinal Analysis (Wiley 2004). I found it informative to vary the parameters of the simulation below to see how well the model works under different circumstances and to understand better how to interpret the model output. Experimentation is fast and easy, because largish datasets (millions of observations) involving hundreds of words can be generated, modeled, summarized, and plotted within a second.

This simulation runs in R. It uses the data.table package for some data manipulation, lme4 for the modeling, and ggplot2 for plotting.

#
# Reality.
#
n.words <- 120               # Total *possible* word forms -- some might not appear
log.odds.base <- -2.5        # Log odds of base form followed by "um"
log.odds.contracted <- -2    # Log odds of contracted form followed by "um"
sigma <- 0.5                 # SD of log odds
kappa <- 0.5                 # Gamma parameter of relative frequency distribution
n <- 2.5e1 * n.words         # Mean observations per word form
#set.seed(17)                 # Set for reproducible results

#-- Create the base forms of the words
base <- replicate(n.words, paste0(sample(letters, sample.int(8)+5, replace=TRUE), collapse=""))

#-- Describe the actual word behavior
f <- function(y) 1 / (1 + exp(-y)) # Convert from log odds to probability
X <- data.frame(
  base=c(base, base),
  form=c(base, abbreviate(base)),
  contracted=c(rep(FALSE, n.words), rep(TRUE, n.words)),
  frequency=rep(rgamma(n.words, kappa), 2),
  prob.pause=f(c(rnorm(n.words, log.odds.base, sigma), 
                 rnorm(n.words, log.odds.contracted, sigma)))
)
#
# Data.
#
p <- X$frequency; p <- p / sum(p)
i <- sample.int(nrow(X), n, replace=TRUE, prob=p)
pause <- runif(n) < X$prob.pause[i]

#-- Create a data table to record the raw simulation results
library(data.table)
df.raw <- data.table(base = X$base[i],
                 form = X$form[i],
                     contracted = X$contracted[i],
                     pause = pause)

#-- Summarize the raw data in a form suitable for analysis
df <- df.raw[, .(pause=sum(pause), pause.not=sum(!pause)), 
             keyby=.(base, form, contracted)]
#
# Model.
#
library(lme4)
fit <- glmer(cbind(pause, pause.not) ~ contracted +  (1 | base:contracted), 
             data=df, family=binomial)
fit.base <- glmer(cbind(pause, pause.not) ~ contracted +  (1 | base), 
             data=df, family=binomial)
#
# Inspections.
#
# model.matrix(fit, type="random")
# model.matrix(fit.2, type="random")
# summary(fit.base)
summary(fit)
anova(fit.base, fit) # Compares competing models
#
# Plots.
#
library(ggplot2)

#-- Calculate predicted and observed rates for plotting
invisible({
  df[, rate:=predict(fit, type="response")]
  df[, observed := pause/(pause + pause.not)]
  df[, size := pause + pause.not]
})

#-- Establish a plotting range comprising all data
prange <- with(df, range(c(range(observed), range(rate))))

#-- Plot predicted vs. observed
ggplot(df, aes(observed, rate, group=contracted)) + 
  coord_fixed(ratio=1, xlim=prange, ylim=prange) + 
  geom_point(aes(fill=contracted, size=size), alpha=1/2, shape=21) + 
  geom_abline(intercept=0, slope=1, size=1) + 
  xlab("Observed pause frequency") + 
  ylab("Predicted pause frequency") + 
  ggtitle("Predicted vs. Observed", "Binomial Mixed Model")

#-- Plot the data
ggplot(df, aes(contracted, observed, fill=contracted)) + 
  geom_violin(alpha=1/3) + 
  geom_point(aes(size=size), alpha=1/2,
             position=position_jitter(height=0.01, width=0.05)) + 
  scale_size_continuous(range=c(0.1, 3)) + 
  ggtitle("Observed Frequencies", "Simulated Data")
$\endgroup$
  • $\begingroup$ Thank you for your thought out answer. I had to read it a few times to comprehend. Minor note: I think you have a typo right after the violin plots where you say "right violin" twice. I also have a question on the violin plots. What does the width of the plot for each row represent? What do you mean when you say 9+26? $\endgroup$ – jdude97 Dec 13 '18 at 17:54
  • $\begingroup$ Thanks for pointing out the typo--I'll fix it. The width is intended to represent an estimated density: these plots are kind of like vertical histograms. The value 9+26 is drawn from the first row of the data, indicating the "size" variable is the total number of observations of each word. $\endgroup$ – whuber Dec 13 '18 at 17:57
  • $\begingroup$ I don't follow what you mean by estimated density. Could you go through your example again of how the violin plot is generated? $\endgroup$ – jdude97 Dec 13 '18 at 18:08
  • $\begingroup$ Rather than trying to give you a full answer in such a limited space, let me refer you to the existing explanations at stats.stackexchange.com/search?q=violin+plot. $\endgroup$ – whuber Dec 13 '18 at 18:09
  • $\begingroup$ I reviewed those now. So I understand the violin plot is like a histogram with the number of counts for each row. Is the vertical axis of observed frequency the frequency with which the word is followed by a pause? I think I had been confused about what the axes represented. $\endgroup$ – jdude97 Dec 13 '18 at 18:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.