Linear Mixed Model (NLME or LMER) - Putting bounds on coefficients/estimates

Question

I am trying to build put bounds on Mixed effect model. I usually use LMER function for my models but I am not able to find any way to put bounds on the coefficients. I tried using LME but even this was not helpful. Can anyone help here?

library(nlme)
library(lmerTest)


myDat = structure(list(Score = c(1.62, 2.18, 2.3, 3.46, 3.85, 4.7, 1.41, 
                                 2.21, 3.32, 2.73, 3.34, 3.27, 2.14, 2.73, 2.74, 3.39, 3.59, 4.01, 
                                 1.81, 1.83, 3.22, 3.64, 3.51, 4.26), Subject = structure(c(1L, 
                                                                                            1L, 1L, 2L, 2L, 2L, 3L, 3L, 3L, 4L, 4L, 4L, 5L, 5L, 5L, 6L, 6L, 
                                                                                            6L, 7L, 7L, 7L, 8L, 8L, 8L), .Label = c("A", "B", "C", "D", "E", 
                                                                                                                                    "F", "G", "H"), class = "factor"), Condition = structure(c(1L, 
                                                                                                                                                                                               1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 
                                                                                                                                                                                               2L, 1L, 1L, 1L, 2L, 2L, 2L), .Label = c("No", "Yes"), class = "factor"), 
                       Time = structure(c(1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 
                                          2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L), .Label = c("1PM", 
                                                                                                              "2PM", "3PM"), class = "factor")), .Names = c("Score", "Subject", 
                                                                                                                                                            "Condition", "Time"), class = "data.frame", row.names = c(NA, 
                                                                                                                                                                                                                      -24L))

m1 <- lmer(Score ~ Condition + Time + Condition*Time + (1 | Subject),
          data = myDat)
fixef(m1)
ranef(m1)


# This is a lame attempt to get positive coeffiecients but I am not sure what objective function, I need to spcify
m2 <- lme(Score ~ Condition + Time + Condition*Time,
          data = myDat, random=c(~1 | Subject)
          # ,
          # nlminb(start = c(0,0,0,0,0),lower = c(0,0,0,0,0))
          )

fixef(m2)
ranef(m2)

https://rstudio-pubs-static.s3.amazonaws.com/108675_02daa6544e584a928b296c3bca9a65d5.html (the comments at the end suggest I might not have worked out all of the bugs ... ??) — Ben Bolker, Dec 27 '21 at 15:27

score 2 · Accepted Answer · answered Dec 24 '21 at 20:30

If you only want to entertain the possibility of positive parameters, I do think Ben Bolker has written about implementing this in lme4. Although, in a Bayesian methodology you can place priors on the parameters with the support you desire. The example below uses the lower bound argument, lb, from the brms package (which has similar formula notation) to encode a folded normal prior. Of course, there are many possible priors with the correct support.

library(lme4)
library(brms)

myDat <-  structure(list(Score = c(1.62, 2.18, 2.3, 3.46, 3.85, 4.7, 1.41, 
                                 2.21, 3.32, 2.73, 3.34, 3.27, 2.14, 2.73, 2.74, 3.39, 3.59, 4.01, 
                                 1.81, 1.83, 3.22, 3.64, 3.51, 4.26), Subject = structure(c(1L, 
                                                                                            1L, 1L, 2L, 2L, 2L, 3L, 3L, 3L, 4L, 4L, 4L, 5L, 5L, 5L, 6L, 6L, 
                                                                                            6L, 7L, 7L, 7L, 8L, 8L, 8L), .Label = c("A", "B", "C", "D", "E", 
                                                                                                                                    "F", "G", "H"), class = "factor"), Condition = structure(c(1L, 
                                                                                                                                                                                               1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 
                                                                                                                                                                                               2L, 1L, 1L, 1L, 2L, 2L, 2L), .Label = c("No", "Yes"), class = "factor"), 
                       Time = structure(c(1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 
                                          2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L), .Label = c("1PM", 
                                                                                                              "2PM", "3PM"), class = "factor")), .Names = c("Score", "Subject", 
                                                                                                                                                            "Condition", "Time"), class = "data.frame", row.names = c(NA, 
## your first try                                                                                                                                                                                                                       -24L))

m1 <- lmer(Score ~ Condition + Time + Condition*Time + (1 | Subject),
           data = myDat)
summary(m1)
# 
# Linear mixed model fit by REML ['lmerMod']
# Formula: Score ~ Condition + Time + Condition * Time + (1 | Subject)
# Data: myDat
# 
# REML criterion at convergence: 25.7
# 
# Scaled residuals: 
#   Min      1Q  Median      3Q     Max 
# -1.4009 -0.5397  0.1392  0.4146  1.2662 

# Random effects:
#   Groups   Name        Variance Std.Dev.
# Subject  (Intercept) 0.05619  0.2370  
# Residual             0.11329  0.3366  
# Number of obs: 24, groups:  Subject, 8
# 
# Fixed effects:
#   Estimate Std. Error t value
# (Intercept)            1.7450     0.2058   8.477
# ConditionYes           1.5600     0.2911   5.359
# Time2PM                0.4925     0.2380   2.069
# Time3PM                1.1500     0.2380   4.832
# ConditionYes:Time2PM  -0.2250     0.3366  -0.668
# ConditionYes:Time3PM  -0.3950     0.3366  -1.174
# 
# Correlation of Fixed Effects:
#   (Intr) CndtnY Tim2PM Tim3PM CY:T2P
# ConditionYs -0.707                            
# Time2PM     -0.578  0.409                     
# Time3PM     -0.578  0.409  0.500              
# CndtnY:T2PM  0.409 -0.578 -0.707 -0.354       
# CndtnY:T3PM  0.409 -0.578 -0.354 -0.707  0.500

## brms - get priors

get_prior(Score ~ Condition + Time + Condition*Time + (1 | Subject),
          data = myDat, family= gaussian())

prior_positive <-  c(
    prior(normal(0, 1), class = Intercept),
    prior(normal(0, 1), class = b, lb = 0), # specify lower bound with "lb"
    #prior(cauchy(0,1), class = shape),
    prior(exponential(1), class = sd),
    prior(exponential(1), class = sigma))

  
mod1_positive <- brm(Score ~ Condition + Time + Condition*Time + (1 | Subject),
        data = myDat, family= gaussian(), prior = prior_positive)
  

summary(mod1_positive)

# Family: gaussian 
# Links: mu = identity; sigma = identity 
# Formula: Score ~ Condition + Time + Condition * Time + (1 | Subject) 
# Data: myDat (Number of observations: 24) 
# Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
# total post-warmup samples = 4000
# 
# Group-Level Effects: 
#   ~Subject (Number of levels: 8) 
# Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# sd(Intercept)     0.26      0.16     0.02     0.65 1.00      954     1653
# 
# Population-Level Effects: 
#   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# Intercept                1.91      0.21     1.48     2.32 1.00     1645     2280
# ConditionYes             1.11      0.27     0.55     1.65 1.00     1431     1742
# Time2PM                  0.27      0.17     0.02     0.65 1.00     1861     1176
# Time3PM                  0.82      0.22     0.37     1.23 1.00     1793     1129
# ConditionYes:Time2PM     0.28      0.21     0.01     0.76 1.00     2183     1490
# ConditionYes:Time3PM     0.24      0.20     0.01     0.75 1.00     1624     1475
# 
# Family Specific Parameters: 
#   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# sigma     0.40      0.08     0.27     0.58 1.00     1379     2100
# 
# Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
# and Tail_ESS are effective sample size measures, and Rhat is the potential
# scale reduction factor on split chains (at convergence, Rhat = 1).

Thanks a ton!! This is really helpful. Can you guide me to the document/link where Ben Bolkar has written/mentioned something for +ve coefficients from LME4? — kawsleo, Dec 27 '21 at 05:50

Linear Mixed Model (NLME or LMER) - Putting bounds on coefficients/estimates

1 Answers1