Why do I get pooled results with regression for a longitudinal dataset? I try to regress the outcomes of soccer matches (outcome number of goals scored) in a Poisson regression framework. My explanatory variables are soccer statistics.
My dataset contains 64 matches, so I've got 128 observations. There're 32 different teams in the database and the teams played different numbers of games (min:3, max:7). So the longitudinal dataset is pretty much unbalanced.
I'm trying to create a Three-level Mixed Effect Poisson model, where the levels are:


*

*individual observations (twice per matches)

*teams random effect

*opponent teams random effect


My problem is that my model doesn't find a solution. It iterates a lot, and displays "log likelihood is not-concave". Then the result are identical with the pooled Poisson regression. I'd like to see the Three-level Mixed Effect Poisson model, but I couldn't get the results of the Two-level model.
What can be the reason why my program doesn't find the Three-level or Two-level model, but the pooled one? Is my sample too small (degrees of freedom problem)? Or do my independent variables cause trouble? However, I could find the Three-level solution with nonsense variables!
Please give me some hints on what could be wrong. 
(I use mainly Stata, but I tried the same in other programs.)
 A: There's a number of possible issues. My best guess is that it is the first issue.
You did not say clearly what kind of algorithm you are using, but I assume some kind of algorithm that tries to maximize the likelihood of the hierarchical model. If you end up getting a (restricted/something else) maximum likelihood estimate of zero for one of the hierarchical scale parameters, the program may simply decide that there are no random effects (i.e. report the pooled model). This may be happening, because you have few observations on a very limited range (soccer matches are typically not very high scoring, so you get a lot of numbers of goals in 0, 1, 2) and you have a number of explanatory variables (how many?) so that these may end up giving you an apparent perfect fit already. You may in fact have so many variables in the model that you almost cannot help but have that situation (did you do a mental check of number of explanatory variables [counting categorical factors as (number of levels $-$ 1) variables, and random effects you want to add as a minimum as two] vs. observations?).
If this is the issue, then you could either get more data, switch to a penalized likelihood approach (see e.g. Chung, Y., Rabe‐Hesketh, S., & Choi, I. H. (2013). Avoiding zero between‐study variance estimates in random‐effects meta‐analysis. Statistics in Medicine 32(23), 4071-4089.) or switch to a Bayesian approach with at least some mild regularization on all model parameters via at least weakly-informative priors (e.g. see this vignette on count data in the R package rstanarm*: https://cran.r-project.org/web/packages/rstanarm/vignettes/count.html  in combination with the one on hierarchical models: https://cran.r-project.org/web/packages/rstanarm/vignettes/pooling.html). 
* I believe there is also Stata Stan (see http://mc-stan.org/interfaces/stata-stan), but I do not know how one uses it.
Another issue could be that the algorithm for fitting the model has trouble to converge and then decides to report on a simpler model (no idea whether that's what Stata would do). In that case, better guesses for initial parameter values may help, another fitting algorithm (not sure whether Stata allows user control there) or different options for the fitting algorithm (e.g. increasing the number of quadrature points) might help.
