# Calculating the standard error of difference in proportions from paired multinomial variables

I have a data set where respondents are placed into one of four categories at two time points. I've worked up an example in R.

library(tidyverse)

set.seed(1234)

example_data <- tibble(id = seq_len(1000),
score_1 = sample(1:4, size = 1000, replace = TRUE),
score_2 = sample(1:4, size = 1000, replace = TRUE))
example_data
#> # A tibble: 1,000 x 3
#>       id score_1 score_2
#>    <int>   <int>   <int>
#>  1     1       4       4
#>  2     2       4       3
#>  3     3       2       4
#>  4     4       2       3
#>  5     5       1       2
#>  6     6       4       2
#>  7     7       3       3
#>  8     8       1       3
#>  9     9       1       2
#> 10    10       2       4
#> # … with 990 more rows


What I am hoping to do is calculate the proportion of respondents in each category at each time point, and then determine if the proportions changed "significantly" between the two time points. To this point, I have this by calculating the proportion of respondents in each category each time point, and then the standard error of that proportion as:

$$se = \frac{p(1-p)}{n}$$

time_1 <- example_data %>%
count(score = score_1, name = "n_1") %>%
mutate(prop_1 = n_1 / sum(n_1),
se_1 = sqrt((prop_1 * (1 - prop_1)) / n_1))

time_2 <- example_data %>%
count(score = score_2, name = "n_2") %>%
mutate(prop_2 = n_2 / sum(n_2),
se_2 = sqrt((prop_2 * (1 - prop_2)) / n_2))


I then calculate the difference between the two proportions, and the standard error of the difference as the square root of the sum of the two variances.

full_join(time_1, time_2, by = "score") %>%
mutate(change = prop_2 - prop_1,
change_se = sqrt(se_1^2 + se_2^2))
#> # A tibble: 4 x 9
#>   score   n_1 prop_1   se_1   n_2 prop_2   se_2   change change_se
#>   <int> <int>  <dbl>  <dbl> <int>  <dbl>  <dbl>    <dbl>     <dbl>
#> 1     1   233  0.233 0.0277   242  0.242 0.0275  0.00900    0.0391
#> 2     2   258  0.258 0.0272   243  0.243 0.0275 -0.015      0.0387
#> 3     3   281  0.281 0.0268   244  0.244 0.0275 -0.037      0.0384
#> 4     4   228  0.228 0.0278   271  0.271 0.027   0.043      0.0387


Created on 2019-12-04 by the reprex package (v0.3.0)

However, it is possible (and actually expected) that the two scores will be correlated between the first and second time point. The current method assumes the two proportions are independent, and thus this assumption is violated.

Is there a way to account for this dependency when looking at the difference in proportions?