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$aveage
Calculate 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.738
Check 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.000
Let’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.203
Check 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.000
Before 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$avesage
Let’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.738
out3m2 <- 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.000
out3m3 <- 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.203
out3m4 <- 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.000
The 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.808
The 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.014
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.01
Calculate the group means and group-mean-centered math achievement scores.
dat4$aveach <- ave(dat4$ach, dat4$schoolid)
dat4$diffach <- dat4$ach - dat4$aveach
Predict 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.000
Find 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.26
The 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.00
Use 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.00
The interaction can be plotted as well.
interact_plot(model=out4m1a, pred=aveach, modx=diffach, modx.values=diffachval)
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$avefemale
First, 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 ' ' 1
Next, 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 ' ' 1
The 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.560
As 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.91
Check 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.000
Next, 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 ' ' 1
Next, 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.2851
The 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 ' ' 1
The 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.5344
The 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.950
Let’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")