eels.r

## Write in the data
o <- rbind(
    c(264, 127, 99),
    c(161, 116, 67)
)
colnames(o) <- c("Border", "Grass", "Sand")
rownames(o) <- c("Moringa", "Vicinus")

## First let's get our chi value
## Expected values
e <- rowSums(o) %o% colSums(o) / sum(o)
## Chi Value
chi_obs <- sum((e-o)^2 / e)

## Estimate the proportions of eels at locations
## These are estimates of the population proportions
## Recall that we have assumed a multinomial distribution
## If that was true, this estimates the population parameters
n_m <- sum(o[1, ])
p_m <- o[1, ] / n_m

n_v <- sum(o[2, ])
p_v <- o[2, ] / n_v

## Simulate the Moringa Eels
moringa_sim <- rmultinom(
    n = 1, # Run 1 experiment
    size = n_m, # How many eels, roll dice 490 times
    prob = p_v
) # Proportion of Moringa Eels

## Simulate Vicinus Eels
vicinus_sim <- rmultinom(
    n = 1, # Run 1 experiment
    size = n_v, # How many eels, roll dice 490 times
    prob = p_m
) # Proportion of Moringa Eels

## Now we create the observed values
o <- rbind(c(moringa_sim), c(vicinus_sim))
## Now we get the expected values we would estimate given that
e <- rowSums(o) %o% colSums(o) / sum(o)
## Now we calculate chi
sum((o - e)^2 / e)



## _________________________________________________
## Below here we will replicate that to get a pvalue
## _________________________________________________


chi_sim <- replicate(10^5, {
    ## Simulate the Moringa Eels
    moringa_sim <- rmultinom(
        n = 1,      # Run 1 experiment
        size = n_m, # How many eels, roll dice 490 times
        prob = p_v) # Proportion of Moringa Eels

    ## Simulate Vicinus Eels
    vicinus_sim <- rmultinom(
        n = 1,      # Run 1 experiment
        size = n_v, # How many eels, roll dice 490 times
        prob = p_m) # Proportion of Moringa Eels

    ## Now we create the observed values
    o <- rbind(c(moringa_sim), c(vicinus_sim))
    ## Now we get the expected values we would estimate given that
    e <- rowSums(o) %o% colSums(o) / sum(o)
    ## Now we calculate chi
    sum((o - e)^2 / e)
})

## How often did we see values more extreme than our observation
## When the H_0 was true
pval <- mean(chi_sim > chi_obs)
round(pval, 2) ## 0.58