2

I have categorized a patient's pattern of screening for a disease (annual, biennial, or else) and would now like to regress patient demographics on this pattern while adjusting for primary care provider (PCP) characteristics. I'm quite sure this requires a multinomial mixed effects model.

My response variable, "Pattern", is a character variable with 3 unordered factors and my grouping variable is "PCP", the PCP's ID. Here's a simplified reproducible example:

df<-data.frame("ID"=seq(1:20),
               "PCP"=rep((seq(1:10)*100),2),
               "Pattern"=rep(c("Annual","Biennial","Biennial","Annual","Else"),4),
               "Age"=runif(20,50,70))

df$PCP<-as.factor(as.character(df$PCP))

When I run what I believe the regression should be I get the following:

mod1<-glmer(Pattern~Age + (1|PCP), data=df)

Error in mkRespMod(fr, REML = REMLpass) : response must be numeric
In addition: Warning message:
In glmer(Pattern ~ Age + (1 | PCP), data = df) :
  calling glmer() with family=gaussian (identity link) as a shortcut to lmer() is deprecated; please call lmer() directly

Any help in troubleshooting would be most appreciated.

benschumacher12
  • 59
  • 1
  • 5
  • alternatives https://stackoverflow.com/questions/21082396/multinomial-logistic-multilevel-models-in-r . Could maybe also use https://stat.ethz.ch/R-manual/R-devel/library/mgcv/html/multinom.html for simple random effects structures – user20650 Jan 02 '20 at 00:31

1 Answers1

4

From the error message, I can see that your code was trying to fit Linear Mixed-Effects Model with family = gaussian. In this case, the dependent variable needs to be numeric but your Pattern variable is a factor. To fit binary (not multinomial) mixed effects models, you may need to define family: 

library(lme4)                               
mod1<-glmer(Pattern~Age + (1|PCP), data=df, family = binomial)
summary(mod1)

As pointed out by @user20650, glmer with family = binomial convert outcome variable into binary. I cannot work out a solution with glmer, but several other packages can be helpful, such as nnet::multinom() and mlogit::mlogit().

mod1 <-nnet::multinom(Pattern ~ Age, data=df)
broom::tidy(mod1)

> broom::tidy(mod1)
# A tibble: 4 x 6
  y.level  term         estimate std.error statistic p.value
  <chr>    <chr>           <dbl>     <dbl>     <dbl>   <dbl>
1 Biennial (Intercept) 0.0000124     6.70     -1.69   0.0916
2 Biennial Age         1.21          0.113     1.69   0.0901
3 Else     (Intercept) 0.000800      7.36     -0.968  0.333 
4 Else     Age         1.12          0.125     0.884  0.377
Zhiqiang Wang
  • 6,206
  • 2
  • 13
  • 27
  • after a bit of faffing around... your model seems to fit `glmer(I(Pattern != "Annual")~Age + (1|PCP), data=df, family = binomial)` so converts the outcome to binary instead of fitting a multinomial – user20650 Jan 02 '20 at 00:57
  • @user20650 Thanks for the information. You probably right. I am trying to test how `glmer` works with multinomial response variables. – Zhiqiang Wang Jan 02 '20 at 01:06