0

I would like to fit a mixed model in R using the library nlme. I have three levels: election year, district, and candidates (unfortunately I cannot share the data). Districts are nested in election year. Now I would like one of my predictor variables from Level 1 to have a random slope at Level 3 (election year) but not at Level 2. I have tried different commands but only one ended up running without an error and I am quite sure it might not be right.

This code works:

m <- lme(y ~ x1 + x2 + xn,
          random = list(~x1-1|Level3, ~1|Level2), data=data, subset = x2 == 0,  method="REML", na.action="na.omit")
summary(m)

This seems to work, but I am worried that the "-1", which I implemented out of frustration, does something to the model I am not intending. In a further step, I introduce an interaction effect between my Level 3 variable and x1:

m <- lme(y ~ x1 + x2 + xn + x1:Level3,
          random = list(~x1-1|Level3, ~1|Level2), data=data, subset = x2 == 0,  method="REML", na.action="na.omit")
summary(m)

Any help would be appreciated! Thank you.

Edit: I would like to give a short example to illustrate my problem. I used the example from Richard Blissett and the data from Here.

library(haven)
library(nlme)
data <- read_dta("https://stats.idre.ucla.edu/stat/examples/imm/imm10.dta")
m1 <- lme (math ~ homework + ses,
           random = list(~ses-1|region, ~1|schid),
           data=data, subset = sex == 1,  
           method="REML", 
           na.action="na.omit")
summary(m1)

# Linear mixed-effects model fit by REML
#   Data: data 
#   Subset: sex == 1 
#        AIC      BIC    logLik
#   938.3002 955.4591 -463.1501

# Random effects:
#  Formula: ~ses - 1 | region
#                  ses
# StdDev: 0.0006488783

#  Formula: ~1 | schid %in% region
#         (Intercept) Residual
# StdDev:    2.761716 8.041486

# Fixed effects:  math ~ homework + ses 
#                Value Std.Error  DF   t-value p-value
# (Intercept) 46.65743 1.5362839 120 30.370316   0e+00
# homework     2.29734 0.5496755 120  4.179441   1e-04
# ses          3.83266 1.0350827 120  3.702753   3e-04
#  Correlation: 
#          (Intr) homwrk
# homework -0.645       
# ses       0.325 -0.264

# Standardized Within-Group Residuals:
#        Min          Q1         Med          Q3         Max 
# -2.43117427 -0.72287038 -0.08035428  0.56193617  3.29134643 

# Number of Observations: 132
# Number of Groups: 
#            region schid %in% region 
#                3                10

Now that I put together this example I realized that the -1 in the random part is taking away the (random) intercept of the Level 3 variable. Same happens with a +0 and works fine using the data above. Unfortunately, this is not what I want with my original data. As soon as I replace the -1 in my original code with a +1 (random Intercepts on level 3), which I want, I get an error message:

Error in logLik.lmeStructInt(lmeSt, lmePars) : 
  NA/NaN/Inf in foreign function call (arg 3)
In addition: There were 50 or more warnings (use warnings() to see the first 50)
# all of the warnings are of the same logic

If that helps: There are 5 categories on Level 3 (5 election years) and 299 units on level 2 (districts) or 1495 (299*5) districts within years, respectively. On Level 1 I have 5800 cases.

roalt102
  • 1
  • 1
  • I know you can't share your data, but can you share an *example*? This would be fairly easy to implement in `lmer`, if that's an option (`(1|district:year) + (1+x1|year)`) ... – Ben Bolker May 10 '23 at 17:03
  • Thank you! Yes, this is an option. All my other models are calculated with `nlme`. Therefore, I wanted to try this as well. However, I suspect even using `lme4` will result in error messages - I now suppose because of a data or fitting problem. I will update with the results of my attempt. – roalt102 May 10 '23 at 20:20
  • It works with `lme4`. However, I get the error message `boundary (singular) fit: see help('isSingular')`. I suppose this results from small random effects? The variance on level 2 is 0. There should not be much collinearity. – roalt102 May 11 '23 at 07:24

0 Answers0