# Which statistical model is suitable?

I have the results of a survey of $$n=132$$ patients with their socio-economic profile and their spending behavior on mobility-coins (my thesis topic).

In the survey, we asked people how they would spend $$100$$ coins for three categories (mobility, donation and selling).

My question is, which statistical model can help me to test out how the spending will change if the $$100$$ coins are replaced by $$100$$ million (suppose for a city like London where we want to target thousands of people).

I'm confused since I have input that varies and has three dependent variables.

• In your survey, did you ask people how they would spend these tokens, or did you observe how they did spent them? I'm alluding here to the intention-action gap. Commented May 2 at 20:31
• @Durden yes we asked them and i already have the result of how much they spenft on each category out of 100 coins Commented May 11 at 18:19
• @yaseen Does the answer below help you? If not, can you provide clarification on what else you need? Commented May 20 at 0:03

If you allow for the fractional spending of coins$$^\dagger$$, then perhaps the Dirichlet distribution would be an appropriate model. The Dirichlet distribution for a $$K$$-dimensional random vector requires that $$X_i \in [0,1]$$ for $$i=1,2,\ldots, K$$ and $$\sum_{i=1}^K X_i = 1$$.

Let $$X_i$$ ($$i=1,2,3)$$ be the proportion (i.e. between $$0$$ and $$1$$) of the coins which are spent on option $$i$$. Then we can assume that \begin{align*} (X_1, X_2, X_3) &\sim \text{Dirichlet}(\alpha_1, \alpha_2, \alpha_3) & \text{in the one-hundred coin case} \\ (X_1, X_2, X_3) &\sim \text{Dirichlet}(\beta_1, \beta_2, \beta_3) & \text{in the one-million coin case.} \end{align*}

Now the goal would be to perform statistical inference to determine if the parameters differ, e.g.,

\begin{align*} &H_0: \alpha_i = \beta_i & \text{for every } i \in \{1,2,3\} \\ &H_1: \alpha_i \neq \beta_i & \text{for some } i\in\{1,2,3\} \end{align*}

For ideas on how to actually perform this test, try Bhattacharya & Dunson (2012) or this dissertation.

## Likelihood Ratio Test in R

One option for conducting this hypothesis test is to use a likelihood ratio test (although Li (2015) suggests that an Energy test may be more powerful). To follow this approach, we need to find the maximum likelihood under both hypotheses.

1. Under the null hypothesis, we combine the two datasets into one and find maximum likelihood estimates of the dirichlet concentration parameters for the combined data. Now let $$\hat\ell_0$$ be the log-likelihood evaluated at these parameter estimates.
2. Under the alternative hypothesis, then the likelihood is maximized by finding MLE's of the two samples individually, denoted $$\hat\alpha$$ and $$\hat\beta$$. If we assume independence between the two samples then the maximum likelihood is given by $$\hat\ell_1 = \ell({\bf X}_1 |\hat\alpha) + \ell({\bf X}_2 |\hat\beta)$$ where $${\bf X}_1$$ and $${\bf X}_2$$ are the first sample (one-hundred coin case) and the second sample (one-million coin case) respectively.

Now the likelihood ratio test is given by $$L = -2\left(\hat\ell_0 - \hat\ell_1\right),$$ which asymptotically follows a Chi-square distribution with $$3$$ degrees of freedom (Wasserman 2013).

To see this in action, we can simulate data with $$n=132$$ and \begin{align*} {\bf\alpha} &= (1, 2, 3) \\ {\bf\beta} &= (1.05, 2, 2.95). \end{align*} For one simulated data-set, we find a test statistic of $$L \approx 6.2396$$ which yields a p-value of roughly $$1.005$$. R code is given below.

library(LaplacesDemon) # For dirichlet density and random generation
library(sirt)          # For dirichlet MLEs

# Make up some "true" parameters
set.seed(1111)
n <- 132
alpha <- c(1, 2, 3)
beta  <- c(1.05, 2, 2.95)

# Simulate data
X <- rdirichlet(n, alpha)
Y <- rdirichlet(n, beta)
Z <- rbind(X, Y)

# Get MLEs under null hypothesis
mle0 <- sirt::dirichlet.mle(Z)
# Get log likelihood under null hypothesis
ll0 <- sum(ddirichlet(Z, mle0\$alpha, log=TRUE))

# Get MLE's under alternative hypothesis
mleX <- sirt::dirichlet.mle(X)
mleY <- sirt::dirichlet.mle(Y)
# Get log likelihood under null hypothesis
llA <- sum(ddirichlet(X, mleX$$alpha, log=TRUE)) + sum(ddirichlet(Y, mleY$$alpha, log=TRUE))

# Calculate likelihood ratio test statistic
L <- -2*(ll0 - llA)

# Get p-value
pval <- pchisq(L, df=3, lower.tail=FALSE)

# Make plot
curve(dchisq(x, df=3), to=2*L,
ylab="density", main=paste0("P-value = ", round(pval, 3)))
abline(h=0)
xx = seq(L, 2*L, length.out=100)
yy = dchisq(xx, df=3)
polygon(c(xx, rev(xx)), c(rep(0, 100), rev(yy)), col='dodgerblue')


## References

Bhattacharya, Abhishek, and David Dunson. "Nonparametric Bayes classification and hypothesis testing on manifolds." Journal of multivariate analysis 111 (2012): 1-19.

Li, Yi. Goodness-of-fit tests for Dirichlet distributions with applications. Bowling Green State University, 2015.

Wasserman, Larry. All of statistics: a concise course in statistical inference. Springer Science & Business Media, 2013.

$$\dagger$$ It may not be that big of a deal if your data doesn't allow for fractional spending. It's hard to say for sure. You could consider checking the "sensitivity" of the results, by randomly adding small numbers to the data (don't forget to normalize so that they sum to one), and checking how much your results change. I would guess that the conclusions you draw will be quite robust to these perturbations.