# Mixed-effects log-linear regression for counts -

I've read other answers but couldn't find exactly what I was looking.

I have generated the following contingency table from my data:

'data.frame':   640 obs. of  6 variables:
$FN : Factor w/ 4 levels "0","1","2","3": 1 2 3 4 1 2 3 4 1 2 ...$ NSE         : Factor w/ 2 levels "nsb","nsm": 1 1 1 1 2 2 2 2 1 1 ...
$Hablante : Factor w/ 2 levels "ADU","CHI": 1 1 1 1 1 1 1 1 2 2 ...$ Destinatario: Factor w/ 2 levels "CDS","OHS": 1 1 1 1 1 1 1 1 1 1 ...
$Sujeto : Factor w/ 20 levels "alma","amadeo",..: 1 1 1 1 1 1 1 1 1 1 ...$ Freq        : int  8 11 7 0 0 0 0 0 0 0 ...


I've read that log-linear poisson regression can be used to model the relationship between categorical (and even ordinal, such as my "FN" column) variables.

It is sometimes point out that observations have to be independent, which doesn't happen in my data. As you can see, there are many observations that correspond to the same subject, so the data has hierarchical structure.

Thus, I created the following model:

a = glmer(Freq ~ FN * (NSE + Hablante + Destinatario) + (1|Sujeto), family= "poisson", data= test2)


and this is what summary returns:

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: poisson  ( log )
Formula: Freq ~ FN * (NSE + Hablante + Destinatario) + (1 | Sujeto)
Data: test2

AIC      BIC   logLik deviance df.resid
4157.3   4233.1  -2061.6   4123.3      623

Scaled residuals:
Min      1Q  Median      3Q     Max
-4.0761 -1.5003 -0.5607  0.3432 12.9388

Random effects:
Groups Name        Variance Std.Dev.
Sujeto (Intercept) 0.2411   0.491
Number of obs: 640, groups:  Sujeto, 20

Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept)          1.57052    0.13203  11.896  < 2e-16 ***
FN1                 -0.40935    0.10811  -3.786 0.000153 ***
FN2                 -1.67113    0.15942 -10.483  < 2e-16 ***
FN3                 -3.48744    0.41943  -8.315  < 2e-16 ***
NSEnsm              -0.23589    0.07742  -3.047 0.002313 **
HablanteCHI         -0.27812    0.07763  -3.583 0.000340 ***
DestinatarioOHS     -0.02077    0.07689  -0.270 0.787036
FN1:NSEnsm           0.52859    0.10882   4.857 1.19e-06 ***
FN2:NSEnsm           0.33486    0.14546   2.302 0.021333 *
FN3:NSEnsm           0.41820    0.43401   0.964 0.335251
FN1:HablanteCHI     -0.23042    0.11013  -2.092 0.036422 *
FN2:HablanteCHI      0.02583    0.14628   0.177 0.859848
FN3:HablanteCHI     -0.70261    0.48372  -1.453 0.146362
FN1:DestinatarioOHS  0.48045    0.10929   4.396 1.10e-05 ***
FN2:DestinatarioOHS  0.92113    0.15595   5.907 3.49e-09 ***
FN3:DestinatarioOHS  0.20314    0.43392   0.468 0.639677
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation matrix not shown by default, as p = 16 > 12.
Use print(x, correlation=TRUE)  or
vcov(x)        if you need it


I have several questions:

• Is it alright to proceed this way? I've found this article (http://users.stat.ufl.edu/~aa/articles/coull_agresti_2003.pdf) by Coull and Agresti that suggests it is possible to fit log-linear models with random effects, but the implementation is on SAS, so I don't know if I'm following the right procedure on R. I had fitted the model without the random effect and it looked really good, but given the variance of the random effect, it is importance to add it, right?

• I've run some diagnostics to assess the model, and I'm not sure the fit is good (but need help on interpreting the results).

qqnorm(resid(a))
qqline(resid(a))


• Perhaps I could get a better fit with backward/forward model selection but I get an error whenever trying to used a mixed-effects model with step(). I know model selection procedures are not recommended in general, but I was surprised to see they are widely used to get insights from log-linear models. How can I perform backward/forward selection on this model?

• Many Freq cells are 0 (449 out of 640), is this a problem ? Can it be solved somehow?

Thanks!