Intro to Multilevel Model: Lecture 7 Example Code
Sunthud Pornprasertmanit
26/01/2024
No-Centering vs. Group-Mean-Centering
Example 1: Customers Satisfaction within Each Table in Dining Services
Read the data set from Lecture 4.
dat3 <-read.table("lecture4ex3.csv", sep=",", header=TRUE)Calculate group means by the ave function. The result
has the same length as the original score and provides the group means
that each data point belongs to. The group-mean centered variable is the
difference between the original score and the group mean.
dat3$aveage <- ave(dat3$age, dat3$tableid)
dat3$diffage <- dat3$age - dat3$aveageCalculate the output when age is not centered. The
regression coefficient of age equals both the between-table
and within-table effects.
library(lme4)
out3m1 <- lmer(sat ~ 1 + age + (1|tableid), data=dat3, REML=FALSE)
summary(out3m1)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: sat ~ 1 + age + (1 | tableid)
## Data: dat3
##
## AIC BIC logLik deviance df.resid
## 16585.8 16608.6 -8288.9 16577.8 2258
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.05873 -0.59420 -0.00433 0.60538 3.07580
##
## Random effects:
## Groups Name Variance Std.Dev.
## tableid (Intercept) 68.44 8.273
## Residual 60.82 7.799
## Number of obs: 2262, groups: tableid, 500
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 55.1116 0.6093 90.45
## age 0.2553 0.0125 20.42
##
## Correlation of Fixed Effects:
## (Intr)
## age -0.738Check when age is group-mean centered. The regression
coefficient represents the within-table effect.
out3m2 <- lmer(sat ~ 1 + diffage + (1|tableid), data=dat3, REML=FALSE)
summary(out3m2)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: sat ~ 1 + diffage + (1 | tableid)
## Data: dat3
##
## AIC BIC logLik deviance df.resid
## 16561.4 16584.3 -8276.7 16553.4 2258
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.14691 -0.59933 -0.00696 0.60899 3.07337
##
## Random effects:
## Groups Name Variance Std.Dev.
## tableid (Intercept) 64.45 8.028
## Residual 60.82 7.798
## Number of obs: 2262, groups: tableid, 500
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 64.28983 0.40129 160.21
## diffage 0.26783 0.01276 20.99
##
## Correlation of Fixed Effects:
## (Intr)
## diffage 0.000Let’s predict the customer satisfaction with group means of
age too. Check the result when both the original scale of
age and the group means of age are the
predictors. The regression coefficient of the original scale of
age is the within-table effect. However, the regression
coefficient of the age table averages is the difference
between within- and between-table effects.
out3m3 <- lmer(sat ~ 1 + age + aveage + (1|tableid), data=dat3, REML=FALSE)
summary(out3m3)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: sat ~ 1 + age + aveage + (1 | tableid)
## Data: dat3
##
## AIC BIC logLik deviance df.resid
## 16562.7 16591.3 -8276.4 16552.7 2257
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1638 -0.6006 -0.0048 0.6144 3.0681
##
## Random effects:
## Groups Name Variance Std.Dev.
## tableid (Intercept) 64.30 8.019
## Residual 60.82 7.799
## Number of obs: 2262, groups: tableid, 500
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 66.10005 2.24652 29.423
## age 0.26783 0.01276 20.992
## aveage -0.31818 0.06279 -5.068
##
## Correlation of Fixed Effects:
## (Intr) age
## age 0.000
## aveage -0.963 -0.203Check the result when both the group-centered age and
the group means of age are the predictors. The regression
coefficient of the group-centered age is the within-table
effect and the regression coefficient of the age table
averages is the between-table effect.
out3m4 <- lmer(sat ~ 1 + diffage + aveage + (1|tableid), data=dat3, REML=FALSE)
summary(out3m4)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: sat ~ 1 + diffage + aveage + (1 | tableid)
## Data: dat3
##
## AIC BIC logLik deviance df.resid
## 16562.7 16591.3 -8276.4 16552.7 2257
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1638 -0.6006 -0.0048 0.6144 3.0681
##
## Random effects:
## Groups Name Variance Std.Dev.
## tableid (Intercept) 64.30 8.019
## Residual 60.82 7.799
## Number of obs: 2262, groups: tableid, 500
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 66.10005 2.24652 29.423
## diffage 0.26783 0.01276 20.992
## aveage -0.05035 0.06148 -0.819
##
## Correlation of Fixed Effects:
## (Intr) diffag
## diffage 0.000
## aveage -0.984 0.000Before trying to find the cross-level interaction of
age, the age needs to be rescaled to make the
multilevel analysis convergent. The scaled age, sage, is
the original scale of age divided by 10.
dat3$sage <- dat3$age / 10
dat3$avesage <- ave(dat3$sage, dat3$tableid)
dat3$diffsage <- dat3$sage - dat3$avesageLet’s rerun the previous analyses.
out3m1 <- lmer(sat ~ 1 + sage + (1|tableid), data=dat3, REML=FALSE)
summary(out3m1)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: sat ~ 1 + sage + (1 | tableid)
## Data: dat3
##
## AIC BIC logLik deviance df.resid
## 16585.8 16608.6 -8288.9 16577.8 2258
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.05873 -0.59420 -0.00433 0.60538 3.07580
##
## Random effects:
## Groups Name Variance Std.Dev.
## tableid (Intercept) 68.44 8.273
## Residual 60.82 7.799
## Number of obs: 2262, groups: tableid, 500
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 55.1116 0.6093 90.45
## sage 2.5532 0.1250 20.42
##
## Correlation of Fixed Effects:
## (Intr)
## sage -0.738out3m2 <- lmer(sat ~ 1 + diffsage + (1|tableid), data=dat3, REML=FALSE)
summary(out3m2)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: sat ~ 1 + diffsage + (1 | tableid)
## Data: dat3
##
## AIC BIC logLik deviance df.resid
## 16561.4 16584.3 -8276.7 16553.4 2258
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.14691 -0.59933 -0.00696 0.60899 3.07337
##
## Random effects:
## Groups Name Variance Std.Dev.
## tableid (Intercept) 64.45 8.028
## Residual 60.82 7.798
## Number of obs: 2262, groups: tableid, 500
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 64.2898 0.4013 160.21
## diffsage 2.6783 0.1276 20.99
##
## Correlation of Fixed Effects:
## (Intr)
## diffsage 0.000out3m3 <- lmer(sat ~ 1 + sage + avesage + (1|tableid),
data=dat3, REML=FALSE)
summary(out3m3)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: sat ~ 1 + sage + avesage + (1 | tableid)
## Data: dat3
##
## AIC BIC logLik deviance df.resid
## 16562.7 16591.3 -8276.4 16552.7 2257
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1638 -0.6006 -0.0048 0.6144 3.0681
##
## Random effects:
## Groups Name Variance Std.Dev.
## tableid (Intercept) 64.30 8.019
## Residual 60.82 7.799
## Number of obs: 2262, groups: tableid, 500
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 66.1000 2.2465 29.423
## sage 2.6783 0.1276 20.992
## avesage -3.1818 0.6279 -5.068
##
## Correlation of Fixed Effects:
## (Intr) sage
## sage 0.000
## avesage -0.963 -0.203out3m4 <- lmer(sat ~ 1 + diffsage + avesage + (1|tableid),
data=dat3, REML=FALSE)
summary(out3m4)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: sat ~ 1 + diffsage + avesage + (1 | tableid)
## Data: dat3
##
## AIC BIC logLik deviance df.resid
## 16562.7 16591.3 -8276.4 16552.7 2257
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1638 -0.6006 -0.0048 0.6144 3.0681
##
## Random effects:
## Groups Name Variance Std.Dev.
## tableid (Intercept) 64.30 8.019
## Residual 60.82 7.799
## Number of obs: 2262, groups: tableid, 500
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 66.1000 2.2465 29.423
## diffsage 2.6783 0.1276 20.992
## avesage -0.5035 0.6148 -0.819
##
## Correlation of Fixed Effects:
## (Intr) diffsg
## diffsage 0.000
## avesage -0.984 0.000The fifth model has sage (the original scale divided by
10), avesage (the group means), and their interaction. As
mentioned in the lecture, the regression coefficients are not intuitive.
The effect of sage represents the within-table effect at
avesage of 0. The effect of avesage represents
the difference between the within-table and between-table effect at
avesage of 0. The interaction effect represents the
increase in both within-table and between-table effects.
out3m5 <- lmer(sat ~ 1 + sage + avesage + sage:avesage
+ (1 + sage|tableid), data=dat3, REML=FALSE)
summary(out3m5)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: sat ~ 1 + sage + avesage + sage:avesage + (1 + sage | tableid)
## Data: dat3
##
## AIC BIC logLik deviance df.resid
## 16566.8 16612.6 -8275.4 16550.8 2254
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1645 -0.5896 -0.0027 0.5992 3.0738
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## tableid (Intercept) 68.4290 8.2722
## sage 0.4891 0.6994 -0.26
## Residual 59.8437 7.7359
## Number of obs: 2262, groups: tableid, 500
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 63.6714 3.7641 16.915
## sage 3.3543 0.8770 3.825
## avesage -2.4887 1.0657 -2.335
## sage:avesage -0.1867 0.2393 -0.781
##
## Correlation of Fixed Effects:
## (Intr) sage avesag
## sage -0.794
## avesage -0.987 0.780
## sage:avesag 0.803 -0.988 -0.808The sixth model has diffsage (the group-centered age),
avesage (the group means), and their interaction. As
mentioned in the lecture, the regression coefficients are easier to
understand. The effect of diffsage represents the
within-table effect at avesage of 0. The effect of
avesage represents the between-table effect. The
interaction effect represents the increase in within-table effects as
avesage increased by 1.
out3m6 <- lmer(sat ~ 1 + diffsage + avesage + diffsage:avesage
+ (1 + diffsage|tableid), data=dat3, REML=FALSE)
summary(out3m6)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: sat ~ 1 + diffsage + avesage + diffsage:avesage + (1 + diffsage |
## tableid)
## Data: dat3
##
## AIC BIC logLik deviance df.resid
## 16568.2 16614.0 -8276.1 16552.2 2254
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1662 -0.6010 -0.0038 0.6128 3.0589
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## tableid (Intercept) 64.4449 8.0278
## diffsage 0.2536 0.5036 -0.11
## Residual 60.2845 7.7643
## Number of obs: 2262, groups: tableid, 500
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 66.09848 2.24635 29.425
## diffsage 2.81671 0.91338 3.084
## avesage -0.50294 0.61473 -0.818
## diffsage:avesage -0.03791 0.24956 -0.152
##
## Correlation of Fixed Effects:
## (Intr) diffsg avesag
## diffsage -0.014
## avesage -0.984 0.014
## diffsag:vsg 0.014 -0.990 -0.014Probing Cross-Level Interactions between the Same Predictor
Example 2: The Effect of Math Achievement on Perceived Self-Efficacy
Read the new data set and check the variables within data.
dat4 <- read.table("lecture7ex1.csv", sep=",", header=TRUE)
library(psych)
describe(dat4)## vars n mean sd median trimmed mad min max range skew
## id 1 1991 996.00 574.90 996 996.00 738.33 1 1991 1990 0.00
## schoolid 2 1991 25.40 14.66 26 25.37 19.27 1 50 49 0.01
## efficacy 3 1991 35.76 3.26 36 35.76 2.97 20 50 30 0.04
## ach 4 1991 54.32 13.77 54 54.18 14.83 17 96 79 0.10
## private 5 1991 0.50 0.50 0 0.50 0.00 0 1 1 0.00
## kurtosis se
## id -1.20 12.88
## schoolid -1.22 0.33
## efficacy 0.96 0.07
## ach -0.35 0.31
## private -2.00 0.01Calculate the group means and group-mean-centered math achievement scores.
dat4$aveach <- ave(dat4$ach, dat4$schoolid)
dat4$diffach <- dat4$ach - dat4$aveachPredict self-efficacy with the group-mean-centered achievement,
diffach, group means of achievement, aveach,
its interaction, and whether the school is private,
private.
out4m1a <- lmer(efficacy ~ 1 + diffach + aveach + private
+ diffach:aveach
+ (1 + diffach|schoolid), data=dat4, REML=FALSE)
summary(out4m1a)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: efficacy ~ 1 + diffach + aveach + private + diffach:aveach +
## (1 + diffach | schoolid)
## Data: dat4
##
## AIC BIC logLik deviance df.resid
## 8888.7 8939.0 -4435.3 8870.7 1982
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1594 -0.6500 0.0071 0.6654 3.2921
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## schoolid (Intercept) 1.53574 1.2392
## diffach 0.02968 0.1723 -0.22
## Residual 4.34903 2.0854
## Number of obs: 1991, groups: schoolid, 50
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 33.168470 1.320101 25.126
## diffach -0.413349 0.140406 -2.944
## aveach 0.048671 0.021855 2.227
## private -0.181827 0.420004 -0.433
## diffach:aveach 0.009505 0.002548 3.730
##
## Correlation of Fixed Effects:
## (Intr) diffch aveach privat
## diffach -0.163
## aveach -0.981 0.176
## private -0.631 0.000 0.525
## diffach:vch 0.160 -0.984 -0.179 0.000Find values for probing interactions. Let’s check the within-school and between-school standard deviations of math achievement.
out4m0 <- lmer(ach ~ 1 + (1|schoolid), data=dat4, REML=FALSE)
summary(out4m0)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: ach ~ 1 + (1 | schoolid)
## Data: dat4
##
## AIC BIC logLik deviance df.resid
## 14905.9 14922.7 -7449.9 14899.9 1988
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.4893 -0.6918 0.0152 0.6705 3.2198
##
## Random effects:
## Groups Name Variance Std.Dev.
## schoolid (Intercept) 92.87 9.637
## Residual 94.97 9.745
## Number of obs: 1991, groups: schoolid, 50
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 54.215 1.381 39.26The result of the null model showed that the school-level standard
deviation was 9.637 and student-level standard deviation was 9.745. The
grand mean of math achievement was 54.215. Thus, aveach
values to examine the within-group effect is 55, \(55 - 10 = 45\), \(55 + 10 = 65\). Thus, diffach
values to examine the between-group effect is -10, 0, 10.
aveachval <- c(45, 55, 65)
diffachval <- c(-10, 0, 10)Use the interactions package to calculate simple slopes.
First, check the within-school effect of math achievement at each level
of schools’ math achievement averages.
library(interactions)
ss41 <- sim_slopes(model=out4m1a, pred=diffach, modx=aveach, modx.values=aveachval)ss41## JOHNSON-NEYMAN INTERVAL
##
## When aveach is OUTSIDE the interval [29.52, 49.25], the slope of diffach is
## p < .05.
##
## Note: The range of observed values of aveach is [38.42, 73.67]
##
## SIMPLE SLOPES ANALYSIS
##
## Slope of diffach when aveach = 45.00:
##
## Est. S.E. t val. p
## ------ ------ -------- ------
## 0.01 0.03 0.41 0.68
##
## Slope of diffach when aveach = 55.00:
##
## Est. S.E. t val. p
## ------ ------ -------- ------
## 0.11 0.03 4.30 0.00
##
## Slope of diffach when aveach = 65.00:
##
## Est. S.E. t val. p
## ------ ------ -------- ------
## 0.20 0.04 5.41 0.00Use the interactions package can visualize the
interaction.
interact_plot(model=out4m1a, pred=diffach, modx=aveach, modx.values=aveachval)
Next, check the simple slope of between-school effect of math achievement at each level of within-school deviation of math achievement.
ss42 <- sim_slopes(model=out4m1a, pred=aveach, modx=diffach, modx.values=diffachval)ss42## JOHNSON-NEYMAN INTERVAL
##
## When diffach is OUTSIDE the interval [-14.77, -0.55], the slope of aveach
## is p < .05.
##
## Note: The range of observed values of diffach is [-44.09, 31.29]
##
## SIMPLE SLOPES ANALYSIS
##
## Slope of aveach when diffach = -10.00:
##
## Est. S.E. t val. p
## ------- ------ -------- ------
## -0.05 0.04 -1.24 0.22
##
## Slope of aveach when diffach = 0.00:
##
## Est. S.E. t val. p
## ------ ------ -------- ------
## 0.05 0.02 2.15 0.04
##
## Slope of aveach when diffach = 10.00:
##
## Est. S.E. t val. p
## ------ ------ -------- ------
## 0.14 0.03 4.58 0.00The interaction can be plotted as well.
interact_plot(model=out4m1a, pred=aveach, modx=diffach, modx.values=diffachval)
Example 1: Customers Satisfaction within Each Table in Dining Services
Let’s examine another example of cross-level interactions. Instead of
age, the cross-level interaction of female is
used. Calculate the group means and group-mean-centered
female.
dat3$avefemale <- ave(dat3$female, dat3$tableid)
dat3$difffemale <- dat3$female - dat3$avefemaleFirst, check whether the effect of female is varied
across tables.
out3g1 <- lmer(sat ~ 1 + female + avefemale
+ (1|tableid), data=dat3, REML=FALSE)
out3g2 <- lmer(sat ~ 1 + female + avefemale
+ (1 + female|tableid), data=dat3, REML=FALSE)
anova(out3g1, out3g2)## Data: dat3
## Models:
## out3g1: sat ~ 1 + female + avefemale + (1 | tableid)
## out3g2: sat ~ 1 + female + avefemale + (1 + female | tableid)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## out3g1 5 16940 16969 -8465.1 16930
## out3g2 7 16937 16977 -8461.3 16923 7.704 2 0.02124 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Next, check whether the table proportions of females can predict the
varying effect of female.
out3g3 <- lmer(sat ~ 1 + female + avefemale + female:avefemale
+ (1 + female|tableid), data=dat3, REML=FALSE)
anova(out3g2, out3g3)## Data: dat3
## Models:
## out3g2: sat ~ 1 + female + avefemale + (1 + female | tableid)
## out3g3: sat ~ 1 + female + avefemale + female:avefemale + (1 + female | tableid)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## out3g2 7 16937 16977 -8461.3 16923
## out3g3 8 16934 16980 -8459.2 16918 4.1516 1 0.04159 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The interaction effect was significant. Check the output.
summary(out3g3)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: sat ~ 1 + female + avefemale + female:avefemale + (1 + female |
## tableid)
## Data: dat3
##
## AIC BIC logLik deviance df.resid
## 16934.4 16980.2 -8459.2 16918.4 2254
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.09389 -0.61477 -0.00542 0.61149 2.74930
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## tableid (Intercept) 55.48 7.448
## female 10.22 3.197 0.06
## Residual 72.89 8.538
## Number of obs: 2262, groups: tableid, 500
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 66.5712 0.9541 69.774
## female -3.8363 1.2455 -3.080
## avefemale -3.6460 1.9185 -1.900
## female:avefemale 4.6837 2.2916 2.044
##
## Correlation of Fixed Effects:
## (Intr) female avefml
## female -0.382
## avefemale -0.884 0.444
## female:vfml 0.417 -0.934 -0.560As explained in the lecture, the between-table effect of females is actually curvilinear. Let’s check the graph.
avefemale <- seq(0, 1, 0.01)
est <- coef(summary(out3g3))[,"Estimate"]
cust <- est[1] + ((est[2] + est[3]) * avefemale) + (est[4]*avefemale*avefemale)
plot(avefemale, cust, type="l", xlim=c(0, 1), ylim=c(60, 70),
xlab="Female Proportion in a Table", ylab="Customer Satisfaction")
The values of 20%, 50%, and 80% females in each table are used to probe the interaction.
avefemaleval <- c(0.2, 0.5, 0.8)Check the simple slopes of the within-table effect of females at each level of the table proportions of females.
ss31 <- sim_slopes(model=out3g3, pred=female, modx=avefemale, modx.values=avefemaleval)ss31## JOHNSON-NEYMAN INTERVAL
##
## When avefemale is OUTSIDE the interval [0.61, 8.24], the slope of female is
## p < .05.
##
## Note: The range of observed values of avefemale is [0.00, 1.00]
##
## SIMPLE SLOPES ANALYSIS
##
## Slope of female when avefemale = 0.20:
##
## Est. S.E. t val. p
## ------- ------ -------- ------
## -2.90 0.84 -3.46 0.00
##
## Slope of female when avefemale = 0.50:
##
## Est. S.E. t val. p
## ------- ------ -------- ------
## -1.49 0.45 -3.35 0.00
##
## Slope of female when avefemale = 0.80:
##
## Est. S.E. t val. p
## ------- ------ -------- ------
## -0.09 0.81 -0.11 0.91Check the interaction plot.
interact_plot(model=out3g3, pred=female, modx=avefemale, modx.values=avefemaleval)
Please not that the effect of the interaction of no-centered females and the table proportions of females is not only represent the moderating within-table effect but also the between-table effect too. The model is equally constrained the within- and between-table moderating effects.
Thus, this model is not easy to understand. The better model in this case is to use group-mean-centered females. First, check the random intercept model.
out3f1 <- lmer(sat ~ 1 + difffemale + avefemale + (1|tableid),
data=dat3, REML=FALSE,
control = lmerControl(optimizer ="Nelder_Mead"))
summary(out3f1)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: sat ~ 1 + difffemale + avefemale + (1 | tableid)
## Data: dat3
## Control: lmerControl(optimizer = "Nelder_Mead")
##
## AIC BIC logLik deviance df.resid
## 16940.3 16968.9 -8465.1 16930.3 2257
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.15838 -0.63360 -0.00586 0.61192 2.75216
##
## Random effects:
## Groups Name Variance Std.Dev.
## tableid (Intercept) 60.12 7.753
## Residual 75.45 8.686
## Number of obs: 2262, groups: tableid, 500
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 65.6444 0.8822 74.409
## difffemale -1.5293 0.4175 -3.663
## avefemale -2.6546 1.5320 -1.733
##
## Correlation of Fixed Effects:
## (Intr) dfffml
## difffemale 0.000
## avefemale -0.892 0.000Next, check whether the within-table effect of females is varied across tables.
out3f2 <- lmer(sat ~ 1 + difffemale + avefemale
+ (1 + difffemale|tableid), data=dat3, REML=FALSE,
control = lmerControl(optimizer ="Nelder_Mead"))
anova(out3f1, out3f2)## Data: dat3
## Models:
## out3f1: sat ~ 1 + difffemale + avefemale + (1 | tableid)
## out3f2: sat ~ 1 + difffemale + avefemale + (1 + difffemale | tableid)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## out3f1 5 16940 16969 -8465.1 16930
## out3f2 7 16938 16978 -8461.9 16924 6.4746 2 0.03927 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Next, check whether the cross-level interaction is significant.
out3f3 <- lmer(sat ~ 1 + difffemale + avefemale + difffemale:avefemale
+ (1 + difffemale|tableid), data=dat3, REML=FALSE,
control = lmerControl(optimizer ="Nelder_Mead"))
anova(out3f2, out3f3)## Data: dat3
## Models:
## out3f2: sat ~ 1 + difffemale + avefemale + (1 + difffemale | tableid)
## out3f3: sat ~ 1 + difffemale + avefemale + difffemale:avefemale + (1 + difffemale | tableid)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## out3f2 7 16938 16978 -8461.9 16924
## out3f3 8 16939 16985 -8461.3 16923 1.1424 1 0.2851The cross-level interaction was not significant. Thus, the variation of within-table effect of females between tables was not significantly explained by the proportion of females in each table. However, the cross-level interaction was significant in the previous analysis using no-centered females. It is possible that the interaction between levels was not equal. The within-table effect did not have interaction whereas the between-table effect has the interaction. Thus, add the quadratic term of tables’ female proportion in the random slope model.
out3f4 <- lmer(sat ~ 1 + difffemale + avefemale + I(avefemale^2)
+ (1 + difffemale|tableid), data=dat3, REML=FALSE,
control = lmerControl(optimizer ="Nelder_Mead"))
anova(out3f2, out3f4)## Data: dat3
## Models:
## out3f2: sat ~ 1 + difffemale + avefemale + (1 + difffemale | tableid)
## out3f4: sat ~ 1 + difffemale + avefemale + I(avefemale^2) + (1 + difffemale | tableid)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## out3f2 7 16938 16978 -8461.9 16924
## out3f4 8 16935 16981 -8459.4 16919 5.0536 1 0.02457 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The quadratic term was significant. Make sure whether the quadratic term of tables’ female proportion predict the random slopes of the within-table effect of females.
out3f5 <- lmer(sat ~ 1 + difffemale + avefemale + I(avefemale^2)
+ avefemale:difffemale + I(avefemale^2):difffemale
+ (1 + difffemale|tableid), data=dat3, REML=FALSE,
control = lmerControl(optimizer ="Nelder_Mead"))
anova(out3f4, out3f5)## Data: dat3
## Models:
## out3f4: sat ~ 1 + difffemale + avefemale + I(avefemale^2) + (1 + difffemale | tableid)
## out3f5: sat ~ 1 + difffemale + avefemale + I(avefemale^2) + avefemale:difffemale + I(avefemale^2):difffemale + (1 + difffemale | tableid)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## out3f4 8 16935 16981 -8459.4 16919
## out3f5 10 16938 16995 -8458.8 16918 1.2531 2 0.5344The interaction effect was not significant. Let’s check the output of the random slope model with quadratic terms.
summary(out3f4)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: sat ~ 1 + difffemale + avefemale + I(avefemale^2) + (1 + difffemale |
## tableid)
## Data: dat3
## Control: lmerControl(optimizer = "Nelder_Mead")
##
## AIC BIC logLik deviance df.resid
## 16934.8 16980.6 -8459.4 16918.8 2254
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.10038 -0.62581 -0.01142 0.60715 2.78700
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## tableid (Intercept) 59.99 7.745
## difffemale 10.38 3.223 0.26
## Residual 72.91 8.538
## Number of obs: 2262, groups: tableid, 500
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 67.6715 1.2386 54.634
## difffemale -1.4859 0.4449 -3.340
## avefemale -13.1943 4.8760 -2.706
## I(avefemale^2) 10.2163 4.5320 2.254
##
## Correlation of Fixed Effects:
## (Intr) dfffml avefml
## difffemale 0.002
## avefemale -0.868 0.032
## I(avefml^2) 0.706 -0.035 -0.950Let’s visualize the quadratic effect of the tables’ female proportion.
est <- coef(summary(out3f4))
avefemale <- seq(0, 1, 0.01)
cust <- est[1] + (est[3] * avefemale) + (est[4]*avefemale*avefemale)
plot(avefemale, cust, type="l", xlim=c(0, 1), ylim=c(60, 70),
xlab="Female Proportion in a Table", ylab="Customer Satisfaction")