2
$\begingroup$

I could use some help Finding the make up of the intercept in a generalized linear mixed model.

FYI, I use a 2013 Macbook Pro with a 2.4 GHz dual-core intel chip, 8 GB of ram, macOS big sur 11.2.2, RStudio Version 1.4.1106, and the R Base Package 4.04.

This is the general model I used: price ~ cut + color + carat + (1 | clarity) + (1 | depth). I used the default prior for the bayesian approach. I used the first 300 rows from the ggplot2::diamonds dataset.

Please note that I took both frequentist (lme4) and bayesian (brms) approaches to analyzing these results.

How to I figure out what levels of the IVs were used to generate the basis of the intercept? Do I just have to track it down logically by looking at summary(mlm_bayes_proper)? Is there code that I can use to find this? Is it already present somewhere, and I just missed it? Is there some other method? Does it differ when using the frequentist v. bayesian approach or are the IV levels of the intercept the same?

The code I used for the analyses and some results are below.

Thanks.





Here is some information on the dataset used:

[1] "diamonds300 dataset info.txt"
[1] "# ---- NOTE: gives dataset info"
[1] " \n       "
[1] "head(diamonds300)"
  carat  cut color clarity depth table price    x    y    z
1  0.23 Good     E     VS1  56.9    65   327 4.05 4.07 2.31
2  0.86 Fair     E     SI2  55.1    69  2757 6.45 6.33 3.52
3  0.84 Fair     G     SI1  55.1    67  2782 6.39 6.20 3.47
4  0.70 Fair     G    VVS1  58.8    66  2797 5.81 5.90 3.44
5  0.76 Fair     G     VS1  59.0    70  2800 5.89 5.80 3.46
6  0.57 Fair     E    VVS1  58.7    66  2805 5.34 5.43 3.16
[1] " \n       "
[1] "str(diamonds300)"
'data.frame':   327 obs. of  10 variables:
 $ carat  : num  0.23 0.86 0.84 0.7 0.76 0.57 0.74 0.91 0.98 0.71 ...
 $ cut    : Ord.factor w/ 5 levels "Fair"<"Good"<..: 2 1 1 1 1 1 1 1 1 1 ...
 $ color  : Ord.factor w/ 7 levels "D"<"E"<"F"<"G"<..: 2 2 4 4 4 2 3 5 2 1 ...
 $ clarity: Ord.factor w/ 8 levels "I1"<"SI2"<"SI1"<..: 5 2 3 7 5 7 4 2 2 4 ...
 $ depth  : num  56.9 55.1 55.1 58.8 59 58.7 61.1 61.3 53.3 56.9 ...
 $ table  : num  65 69 67 66 70 66 68 67 67 65 ...
 $ price  : int  327 2757 2782 2797 2800 2805 2805 2825 2855 2858 ...
 $ x      : num  4.05 6.45 6.39 5.81 5.89 5.34 5.82 6.24 6.82 5.89 ...
 $ y      : num  4.07 6.33 6.2 5.9 5.8 5.43 5.75 6.19 6.74 5.84 ...
 $ z      : num  2.31 3.52 3.47 3.44 3.46 3.16 3.53 3.81 3.61 3.34 ...
NULL
[1] " \n       "
[1] "colnames(diamonds300)"
 [1] "carat"   "cut"     "color"   "clarity" "depth"   "table"   "price"   "x"       "y"       "z"      
[1] " \n       "
[1] "nrow(diamonds300)"
[1] 327
[1] " \n       "
[1] "# ---- NOTE: gives unique values of Fixed and Random effects"
[1] "unique(diamonds300$cut)"
[1] Good      Fair      Very Good
Levels: Fair < Good < Very Good < Premium < Ideal
[1] " \n       "
[1] "unique(diamonds300$color)"
[1] E G F H D I J
Levels: D < E < F < G < H < I < J
[1] " \n       "
[1] "unique(diamonds300$carat)"
 [1] 0.23 0.86 0.84 0.70 0.76 0.57 0.74 0.91 0.98 0.71 0.75 0.72 0.88 0.90 0.99 1.06 0.85 0.73 1.20 0.92 1.00 0.95
[23] 1.01 0.96 1.14 1.02 0.61 1.17 1.18 0.94 0.97 0.31 1.16 1.05 1.45 0.51 1.07 1.13 0.30 0.93 1.24 1.04 1.35 1.65
[45] 1.50 1.21 1.42 1.56 2.01 1.44 1.51 1.57 1.52 1.53 1.76 1.62 1.55 2.00 3.00 2.10 0.29 1.91 1.32 2.29 1.98 2.03
[67] 2.51 2.48 0.45 0.36 0.50 0.37 0.46 0.25 0.35 0.40 0.43 0.24 0.42 0.54 0.62 0.49 0.56 0.89 0.68 0.53 0.52 0.64
[89] 0.67 0.63 0.55 0.77 0.69 0.60 0.81 0.79 0.82 0.80 0.78
[1] " \n       "
[1] "unique(diamonds300$clarity)"
[1] VS1  SI2  SI1  VVS1 VS2  I1   VVS2 IF  
Levels: I1 < SI2 < SI1 < VS2 < VS1 < VVS2 < VVS1 < IF
[1] " \n       "
[1] "unique(diamonds300$depth)"
  [1] 56.9 55.1 58.8 59.0 58.7 61.1 61.3 53.3 55.8 56.6 64.1 56.0 58.0 61.6 59.1 61.0 57.7 58.6 62.2 62.5 62.4 57.4
 [23] 64.6 57.5 60.9 60.7 58.5 61.5 58.1 60.1 61.8 60.8 57.6 60.0 65.7 65.6 58.9 59.8 59.2 57.0 59.3 56.8 56.4 56.5
 [45] 55.2 64.8 62.6 60.2 64.9 59.5 62.0 60.6 55.9 61.7 63.1 61.9 61.4 56.3 60.4 68.5 57.8 59.6 55.6 58.4 67.3 59.9
 [67] 62.1 59.7 56.7 58.2 55.3 64.2 58.3 57.3 61.2 62.7 59.4 60.3 57.9 62.8 63.4 51.0 54.2 52.7 62.9 54.3 65.5 57.2
 [89] 56.1 62.3 55.0 52.2 63.6 53.4 63.9 68.8 68.2 63.0 65.3 56.2 65.8 55.5 79.0 54.7 64.5 64.3
[1] " \n       "
[1] "unique(diamonds300$table)"
 [1] 65.0 69.0 67.0 66.0 70.0 68.0 95.0 71.0 73.0 65.4 79.0 76.0
[1] " \n       "



Here is the R Code I used to create the model:



# generalized linear mixed models

## packages used
# ---- NOTE: for loading diamonds dataset
if(!require(ggplot2)){install.packages("ggplot2")}
# ---- NOTE: for interpreting mixed effect models
if(!require(jtools)){install.packages("jtools")}
# ---- NOTE: for bayes modeling
if(!require(rstan)){install.packages("rstan")}
# ---- NOTE: for bayes modeling
if(!require(brms)){install.packages("brms")}
# ---- NOTE: run mixed effects models
if(!require(lme4)){install.packages("lme4")}
# ---- NOTE: run mixed effects models comparisons
if(!require(lsmeans)){install.packages("lsmeans")}
# ---- NOTE: run mixed effects models comparisons
if(!require(emmeans)){install.packages("emmeans")}
# ---- NOTE: data wrangling
if(!require(tidyverse)){install.packages("tidyverse")}
# ---- NOTE: for mixed effect models
if(!require(car)){install.packages("car")}

### dataset
# ---- NOTE: selects only the top 300 rows of the dataset
diamonds300 <- data.frame(top_n(diamonds, 300, table))
# ---- NOTE: gives dataset info
head(diamonds300)
str(diamonds300)
colnames(diamonds300)
nrow(diamonds300)
# ---- NOTE: gives unique values of Fixed and Random effects
unique(diamonds300$cut)
unique(diamonds300$color)
unique(diamonds300$carat)
unique(diamonds300$clarity)
unique(diamonds300$depth)
unique(diamonds300$table)

## DV is price
# ---- NOTE: FIXED EFFECTS MAIN IV - cut
# ---- NOTE: FIXED EFFECTS CONTROLLED VARIABLES - color + carat
# ---- NOTE: RANDOM EFFECTS - (1 | clarity) + (1 | depth)
# ---- NOTE: full variable model - price ~ cut + color + carat + (1 | clarity) + (1 | depth)

### frequentist model

#### full model
# ---- NOTE: creates model
(mlm_freq = lme4::glmer(
  price ~ cut + color + carat + (1 | clarity) + (1 | depth),
  data = diamonds300,
  family = 'poisson'))
# ---- NOTE: gives model summary
summary(mlm_freq) 
# ---- NOTE: summary of model, with p values for F-statistic for fixed effects
Anova(mlm_freq)

##### Gives exponentialized table of estimates
summ(mlm_freq, exp = T)

#### prints frequentist results
sink("diamonds300 frequentist.txt")
print("diamonds300 frequentist.txt")
print(" 
       ")
print("mlm_freq")
print(mlm_freq)
print(" 
       ")
print("summary(mlm_freq) ")
print(summary(mlm_freq) )
print(" 
       ")
print("Anova(mlm_freq)")
print(Anova(mlm_freq))
print(" 
       ")
print("summ(mlm_freq, exp = T)")
print(summ(mlm_freq, exp = T))
print(" 
       ")
sink()

### bayesian approach

#### creates bayes model, with proper fixed and random effects
# ---- NOTE: used default priors
mlm_bayes_proper = brms::brm(
  price ~ cut + color + carat + (1 | clarity) + (1 | depth),
  data = diamonds300,
  family = 'poisson'
) 
# ---- NOTE: gives summary of model
summary(mlm_bayes_proper)
# ---- NOTE: gives 95% credible intervals, which can be used as a significance test for levels of fixed effect when compared to intercept, because it gives odds changes (see decimal points with 1 as the +/- point and https://www.rensvandeschoot.com/tutorials/generalised-linear-models-with-brms/)
exp(fixef(mlm_bayes_proper)[,-2])

#### prints bayesian results
sink("diamonds300 bayesian.txt")
print("diamonds300 bauesian.txt")
print(" 
       ")
print("mlm_bayes_proper")
print(mlm_bayes_proper)
print(" 
       ")
print("summary(mlm_bayes_proper)")
print(summary(mlm_bayes_proper))
print(" 
       ")
print("exp(fixef(mlm_bayes_proper)[,-2])")
print(exp(fixef(mlm_bayes_proper)[,-2]))
print(" 
       ")
sink()




### dataset
# ---- NOTE: selects only the top 300 rows of the dataset
diamonds300 <- data.frame(top_n(diamonds, 300, table))
# ---- NOTE: gives dataset info
head(diamonds300)
str(diamonds300)
colnames(diamonds300)
nrow(diamonds300)
# ---- NOTE: gives unique values of Fixed and Random effects
unique(diamonds300$cut)
unique(diamonds300$color)
unique(diamonds300$carat)
unique(diamonds300$clarity)
unique(diamonds300$depth)
unique(diamonds300$table)

#### prints dataset info
sink("diamonds300 dataset.txt")
print("diamonds300 dataset info.txt")
print("# ---- NOTE: gives dataset info")
print(" 
       ")
print("head(diamonds300)")
print(head(diamonds300))
print(" 
       ")
print("str(diamonds300)")
print(str(diamonds300))
print(" 
       ")
print("colnames(diamonds300)")
print(colnames(diamonds300))
print(" 
       ")
print("nrow(diamonds300)")
print(nrow(diamonds300))
print(" 
       ")
print("# ---- NOTE: gives unique values of Fixed and Random effects")
print("unique(diamonds300$cut)")
print(unique(diamonds300$cut))
print(" 
       ")
print("unique(diamonds300$color)")
print(unique(diamonds300$color))
print(" 
       ")
print("unique(diamonds300$carat)")
print(unique(diamonds300$carat))
print(" 
       ")
print("unique(diamonds300$clarity)")
print(unique(diamonds300$clarity))
print(" 
       ")
print("unique(diamonds300$depth)")
print(unique(diamonds300$depth))
print(" 
       ")
print("unique(diamonds300$table)")
print(unique(diamonds300$table))
print(" 
       ")
sink()

```
$\endgroup$
2
  • $\begingroup$ Did you just re-ask your question that was closed without modifications? $\endgroup$ – Arya McCarthy Mar 4 at 5:31
  • $\begingroup$ @AryaMcCarthy I'm not too sure why the other question was closed as it wasn't really a duplicate. $\endgroup$ – Robert Long Mar 4 at 11:02
1
$\begingroup$

Some classes of objects will provide ways to do this, so it depends a lot on which package/function you used, but such a question would be off-topic here.

So without it being about software, if you have default contrast coding it is very easy to know the reference level, because it is simply the one which doesn't have it's own estimate.

$\endgroup$
2
  • $\begingroup$ Cool. Thanks for the clarification. $\endgroup$ – Mel Mar 4 at 19:36
  • $\begingroup$ Thanks. I'm still new. $\endgroup$ – Mel Mar 4 at 20:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.