### Tire Example

library(tidyverse)

dat <- c(2,1,2,3)
expect <- rep(1.75, 4)

obs.stat <- sum( (dat - expect)^2 )

# Simulate one sample
dat.sim <- rmultinom(1, 7, rep(.25,4))
sim.stat <- sum( (dat.sim - expect)^2 )

# Repeat
sim.stats <- NULL
for(i in 1:10000){
    sim.sample <- rmultinom(1, 7, rep(.25,4))
    sim.stats[i] <- sum( (sim.sample - expect)^2 )
}

# p-value
mean(sim.stats >= obs.stat)
# or
sum(sim.stats >= obs.stat)/10000

# p-value = area to the right of our observed statistic:
hist(sim.stats)
abline(v = obs.stat, col = "red")