Using standardization with small noisy values

Question

Let's say we have a dataset of how often people say certain words that looks like

+-------+------+-----+---------+
|       | this | cat | uranium |
+-------+------+-----+---------+
| Alice |   32 |   8 |    0.05 |
| Bob   |   22 |   3 |    0.01 |
| Carol |   25 |   4 |   0.005 |
+-------+------+-----+---------+

If we want to compare the relative usage of the words we can standardize the data's mean and standard deviation to get

+-------+-------+-------+---------+
|       | this  |  cat  | uranium |
+-------+-------+-------+---------+
| Alice |  1.35 |  1.38 |    1.40 |
| Bob   | -1.03 | -0.92 |   -0.57 |
| Carol | -0.31 | -0.46 |   -0.82 |
+-------+-------+-------+---------+

Now we can see that Alice says the word "cat" more often than the word "this" compared to other people, even though the original value of "cat" was much lower than "this". However it also indicates that she says "uranium" significantly more than others, even though the original value for "uranium" was very small. She probably just happened to say it a few times more by chance, and doesn't actually talk about uranium much more than others.

Is there a way to use standardization that takes into account the fact that small values are more susceptible to random noise?

Answer 1

Log linear model involves assuming that the expected counts are a linear function of the person, the word. You would do something like: $E[Y] = \beta_0 + \beta_1 x_1 + \beta_2 x_2$. You assume that the $Y$ are distributed as a multinomial, but it is possible to fit a Poisson model so that you have.

y <- c(3200, 2200, 2500, 800, 300, 400, 5, 1, 1)
x1 <- gl(3, 3, labels=c("Alice", "Bob", "Carol"))
x2 <- factor(rep(c("this", "cat", "uranium"), 3))
m1 <- glm(y ~ x1 + x2, family = poisson())
summary(m1)

You get output like this:

Coefficients:
            Estimate Std. Error z value   
(Intercept)  7.64985    0.02050 373.222
x1Bob       -1.66140    0.02816 -58.989
x1Carol     -7.02871    0.37812 -18.589
x2this       0.47085    0.02549  18.475
x2uranium    0.14836    0.02729   5.437

So the intercept is the log of the expected value of the number of times Alice says cat, i.e. $exp(7.6)$ about 2900. The other coefficients, because they are on a log scale are telling you how much more likely Bob is to say cat than Alice (exp(-1.66) times as much) or how much more likely Alice is to say this than cat (exp(0.47) times as much). The longest one is that Carol is exp(-7.02871 + 0.14836) times as likely to say uranium than Alice is to say cat.

However, I think for your question you want to fit a model with interactions.

$E[Y] = \beta_0 + \beta_1 x_1 + \beta_2 x_2$ + \beta_{12} x_1 * x_2$

With those made up numbers I get:

(Intercept)        7.69621    0.02132 360.984
x1Bob             -1.99243    0.06155 -32.373
x1Carol           -7.69621    1.00022  -7.695
x2this             0.37469    0.02770  13.529
x2uranium          0.12783    0.02923   4.373
x1Bob:x2this       0.60614    0.07315   8.287
x1Carol:x2this     1.23474    1.09579   1.127    
x1Bob:x2uranium    0.15985    0.08178   1.955  
x1Carol:x2uranium -0.12783    1.41451  -0.090    

Now this is where life gets interesting. You're not really interested (I think??????) in whether Bob talks more than Alice, or whether people say cat more often than uranium. So these interactions show that given Bob being more/less talkative than Alice and this being more/less popular than cat, Bob says this exp(0.6) 1.8 times more often than you can account for. Carol says uranium exp(-0.12) 0.88 times as often, again conditioned on uranium already being a less popular word, and Carol being more/less talkative.

Another way of examining this is through graphical models, but Alan Agresti Analysis of Categorical Data is my goto guide for these kinds of questions. There's a lot more detail, but I hope that at least lets you decide whether log-linear model is the right way to go for your problem.

Answer 2

I suppose I should give an update with what I've done. After reading through the textbook Paul suggested, I've settled on using the standardized residuals of the chi squared test.

With the dataset

       this  cat  uranium
Alice  3200  800        5
Bob    2200  300        1
Carol  2500  400        1

using regular standardization gives

          this     cat    uranium  
 ------- ------- ------- --------- 
  Alice    1.35    1.38      1.41  
  Bob     -1.03   -0.92     -0.70  
  Carol   -0.31   -0.46     -0.70  

putting a lot of unwanted emphasis on the smaller values. Using the standardized residuals of the chi squared test gives

           this       cat   uranium
Alice -9.289214  9.191701  1.544499
Bob    6.340982 -6.298057 -0.736932
Carol  3.879593 -3.816253 -0.948655

which treats the variation between bigger numbers as more than the variation between smaller numbers. The scores seem pretty high, but I was just comparing the values between categories so for my use case it works well.

The R code to do this is

y <- c(3200, 2200, 2500, 800, 300, 400, 5, 1, 1)
x1 <- gl(3, 3, labels=c("this", "cat", "uranium"))
x2 <- factor(rep(c("Alice", "Bob", "Carol"), 3))
chisq.test(xtabs(y ~ x2 + x1))$stdres

and python code is

import pandas as pd
import statsmodels.api as sm

df = pd.DataFrame(
    {"this": [3200, 2200, 2500], "cat": [800, 300, 400], "uranium": [5, 1, 1]},
    index=["Alice", "Bob", "Carol"],
)
table = sm.stats.Table(df)
print(table.standardized_resids)

Side note: I also tried the standardized residuals of a log linear model, but that gave the exact same results as the chi squared with much more computation time.