Intro to Multilevel Model: Lecture 6 Example Code
Sunthud Pornprasertmanit
23/01/2024
Cross-level Interactions
Example 1: Interviewees’ Rating Scores within Interviewers
Read the data set from Lecture 4.
dat2 <- read.table("lecture4ex2.csv", sep=",", header=TRUE)First, check whether the interviewers’ sex predict the random intercept.
library(lme4)
out2m1 <- lmer(score ~ 1 + eesex + I(iq - 100) + (1 + eesex|erid), data=dat2, REML=FALSE)
out2m1a <- lmer(score ~ 1 + eesex + I(iq - 100) + ersex + (1 + eesex|erid), data=dat2, REML=FALSE)
anova(out2m1, out2m1a)## Data: dat2
## Models:
## out2m1: score ~ 1 + eesex + I(iq - 100) + (1 + eesex | erid)
## out2m1a: score ~ 1 + eesex + I(iq - 100) + ersex + (1 + eesex | erid)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## out2m1 7 61669 61719 -30828 61655
## out2m1a 8 61669 61727 -30826 61653 2.0971 1 0.1476Even though the interviewers’ sex effect on interview scores was not
significant, the interviewers’ sex will be used to predict random slopes
of interviewees’ sex. For running likelihood ratio test, in the full
model, out2m1b, interviewers’ sex predicts both random
intercept and random slope while, in the restricted model,
out2m1a, interviewers’ sex predicts only random
intercept.
out2m1b <- lmer(score ~ 1 + eesex + I(iq - 100) + ersex + eesex:ersex
+ (1 + eesex|erid), data=dat2, REML=FALSE)
anova(out2m1a, out2m1b)## Data: dat2
## Models:
## out2m1a: score ~ 1 + eesex + I(iq - 100) + ersex + (1 + eesex | erid)
## out2m1b: score ~ 1 + eesex + I(iq - 100) + ersex + eesex:ersex + (1 + eesex | erid)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## out2m1a 8 61669 61727 -30826 61653
## out2m1b 9 61641 61706 -30812 61623 29.87 1 4.619e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Run the summary function, the fixed effect of
interaction was significant.
summary(out2m1b)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: score ~ 1 + eesex + I(iq - 100) + ersex + eesex:ersex + (1 +
## eesex | erid)
## Data: dat2
##
## AIC BIC logLik deviance df.resid
## 61641.0 61705.9 -30811.5 61623.0 9991
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.7413 -0.6367 -0.0048 0.6410 3.4320
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## erid (Intercept) 20.1152 4.4850
## eesex 0.7979 0.8933 0.12
## Residual 21.7645 4.6652
## Number of obs: 10000, groups: erid, 1000
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 75.538700 0.221217 341.469
## eesex -0.413234 0.137868 -2.997
## I(iq - 100) 0.029418 0.003299 8.918
## ersex 0.009887 0.312847 0.032
## eesex:ersex 1.073648 0.194981 5.506
##
## Correlation of Fixed Effects:
## (Intr) eesex I(-100 ersex
## eesex -0.253
## I(iq - 100) -0.003 0.001
## ersex -0.707 0.179 0.002
## eesex:ersex 0.179 -0.707 -0.008 -0.253The effects package is used to calcuate the expected
means of interview ratings in each combination of interviewer’s and
interviewee’s sex.
library(effects)
eff2m1b <- effect(term="eesex:ersex", mod=out2m1b,
xlevels=list(ersex=c(0, 1), eesex=c(0, 1)))
summary(eff2m1b)##
## eesex*ersex effect
## ersex
## eesex 0 1
## 0 75.54466 75.55455
## 1 75.13143 76.21496
##
## Lower 95 Percent Confidence Limits
## ersex
## eesex 0 1
## 0 75.11103 75.12092
## 1 74.68215 75.76568
##
## Upper 95 Percent Confidence Limits
## ersex
## eesex 0 1
## 0 75.97829 75.98818
## 1 75.58070 76.66424The expected means are used to make a clustered bar graph using the
ggplot2 package.
library(ggplot2)
mydat <- as.data.frame(eff2m1b)
mydat$eesex <- factor(mydat$eesex, labels=c("Male", "Female"))
mydat$ersex <- factor(mydat$ersex, labels=c("Male", "Female"))
ggplot(mydat, aes(factor(eesex), fit, fill = ersex)) +
geom_bar(stat="identity", position = "dodge") +
labs(x = "Interviewee", y = "Interview Ratings", fill = "Interviewer")
Use the interactions package to calculate simple slopes.
First, check the effect of interviewees’ sex at each sex of
interviewers.
library(interactions)
ss21 <- sim_slopes(model=out2m1b, pred=eesex, modx=ersex)ss21## JOHNSON-NEYMAN INTERVAL
##
## When ersex is OUTSIDE the interval [0.17, 0.56], the slope of eesex is p <
## .05.
##
## Note: The range of observed values of ersex is [0.00, 1.00]
##
## SIMPLE SLOPES ANALYSIS
##
## Slope of eesex when ersex = 0.00 (0):
##
## Est. S.E. t val. p
## ------- ------ -------- ------
## -0.41 0.14 -2.99 0.00
##
## Slope of eesex when ersex = 1.00 (1):
##
## Est. S.E. t val. p
## ------ ------ -------- ------
## 0.66 0.14 4.78 0.00The interactions package can visualize the interaction
as well. However, it provides the line plot which might not be
appropriate when either the predictor or the moderator is
categorical.
interact_plot(model=out2m1b, pred=eesex, modx=ersex)## Using data dat2 from global environment. This could cause incorrect results
## if dat2 has been altered since the model was fit. You can manually provide
## the data to the "data =" argument.
Next, check the simple slope effect of interviewers’ sex at each sex of interviewees.
ss22 <- sim_slopes(model=out2m1b, pred=ersex, modx=eesex)ss22## JOHNSON-NEYMAN INTERVAL
##
## When eesex is OUTSIDE the interval [-0.68, 0.55], the slope of ersex is p <
## .05.
##
## Note: The range of observed values of eesex is [0.00, 1.00]
##
## SIMPLE SLOPES ANALYSIS
##
## Slope of ersex when eesex = 0.00 (0):
##
## Est. S.E. t val. p
## ------ ------ -------- ------
## 0.01 0.31 0.03 0.97
##
## Slope of ersex when eesex = 1.00 (1):
##
## Est. S.E. t val. p
## ------ ------ -------- ------
## 1.08 0.32 3.34 0.00The interaction can be plotted as well.
interact_plot(model=out2m1b, pred=ersex, modx=eesex)## Using data dat2 from global environment. This could cause incorrect results
## if dat2 has been altered since the model was fit. You can manually provide
## the data to the "data =" argument.
Example 2: Goldfish Food Consumption within Fish Tanks
Read the data set from Lecture 4.
dat1 <-read.table("lecture4ex1.csv", sep=",", header=TRUE)First, check whether the tank volume and the number of goldfishes predict the random intercept.
out1m1 <- lmer(consume ~ 1 + I(length - 15) + goldcolor + (1 + I(length - 15)|groupid),
data=dat1, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
out1m1intcept <- lmer(consume ~ 1 + I(length - 15) + goldcolor + I(volume - 1) + I(nfish - 8)
+ (1 + I(length - 15)|groupid), data=dat1, REML=FALSE,
control = lmerControl(optimizer ="Nelder_Mead"))
anova(out1m1, out1m1intcept)## Data: dat1
## Models:
## out1m1: consume ~ 1 + I(length - 15) + goldcolor + (1 + I(length - 15) | groupid)
## out1m1intcept: consume ~ 1 + I(length - 15) + goldcolor + I(volume - 1) + I(nfish - 8) + (1 + I(length - 15) | groupid)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## out1m1 7 4109.4 4142.3 -2047.7 4095.4
## out1m1intcept 9 4072.2 4114.4 -2027.1 4054.2 41.214 2 1.123e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Use the tank volume and the number of goldfishes to predict random
slopes of the fish lengths. For running likelihood ratio test, in the
full model, out1m1slope, the tank-level predictors predict
both random intercept and random slope while, in the restricted model,
out1m1intcept, the tank-level predictors predict only
random intercept.
out1m1slope <- lmer(consume ~ 1 + I(length - 15) + goldcolor + I(volume - 1) + I(nfish - 8)
+ I(length - 15):I(volume - 1) + I(length - 15):I(nfish - 8)
+ (1 + I(length - 15)|groupid),
data=dat1, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
anova(out1m1intcept, out1m1slope)## Data: dat1
## Models:
## out1m1intcept: consume ~ 1 + I(length - 15) + goldcolor + I(volume - 1) + I(nfish - 8) + (1 + I(length - 15) | groupid)
## out1m1slope: consume ~ 1 + I(length - 15) + goldcolor + I(volume - 1) + I(nfish - 8) + I(length - 15):I(volume - 1) + I(length - 15):I(nfish - 8) + (1 + I(length - 15) | groupid)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## out1m1intcept 9 4072.2 4114.4 -2027.1 4054.2
## out1m1slope 11 4070.8 4122.5 -2024.4 4048.8 5.3731 2 0.06811 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Even though the likelihood ratio test was not significant, let’s check each interaction term.
summary(out1m1slope)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: consume ~ 1 + I(length - 15) + goldcolor + I(volume - 1) + I(nfish -
## 8) + I(length - 15):I(volume - 1) + I(length - 15):I(nfish -
## 8) + (1 + I(length - 15) | groupid)
## Data: dat1
## Control: lmerControl(optimizer = "Nelder_Mead")
##
## AIC BIC logLik deviance df.resid
## 4070.8 4122.5 -2024.4 4048.8 800
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2657 -0.6586 0.0523 0.6363 2.8170
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## groupid (Intercept) 12.34246 3.5132
## I(length - 15) 0.01195 0.1093 0.67
## Residual 5.98796 2.4470
## Number of obs: 811, groups: groupid, 100
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 43.786852 0.717966 60.987
## I(length - 15) 0.928289 0.043326 21.426
## goldcolor -0.490354 0.183360 -2.674
## I(volume - 1) 15.117734 2.051729 7.368
## I(nfish - 8) -1.117353 0.211192 -5.291
## I(length - 15):I(volume - 1) 0.224905 0.123540 1.821
## I(length - 15):I(nfish - 8) -0.004754 0.013094 -0.363
##
## Correlation of Fixed Effects:
## (Intr) I(-15) gldclr I(v-1) I(n-8) I(-1-1
## I(lngth-15) 0.352
## goldcolor -0.135 0.008
## I(volume-1) 0.851 0.313 -0.008
## I(nfish-8) -0.669 -0.249 0.004 -0.759
## I(-15):I(-1 0.321 0.843 -0.016 0.396 -0.291
## I(-15):I(-8 -0.238 -0.750 -0.017 -0.276 0.360 -0.762The the tank volume alone significantly predicted the random slope but the number of fishes could not significantly predict the random slope. Then, the tank volume effect was retained and the effect of the number of fishes was dropped.
out1m1trimmed <- lmer(consume ~ 1 + I(length - 15) + goldcolor + I(volume - 1) + I(nfish - 8)
+ I(volume - 1):I(length - 15) + (1 + I(length - 15)|groupid),
data=dat1, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
anova(out1m1intcept, out1m1trimmed)## Data: dat1
## Models:
## out1m1intcept: consume ~ 1 + I(length - 15) + goldcolor + I(volume - 1) + I(nfish - 8) + (1 + I(length - 15) | groupid)
## out1m1trimmed: consume ~ 1 + I(length - 15) + goldcolor + I(volume - 1) + I(nfish - 8) + I(volume - 1):I(length - 15) + (1 + I(length - 15) | groupid)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## out1m1intcept 9 4072.2 4114.4 -2027.1 4054.2
## out1m1trimmed 10 4068.9 4115.9 -2024.5 4048.9 5.2431 1 0.02203 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Check the output of the cross-level interaction model.
summary(out1m1trimmed)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: consume ~ 1 + I(length - 15) + goldcolor + I(volume - 1) + I(nfish -
## 8) + I(volume - 1):I(length - 15) + (1 + I(length - 15) | groupid)
## Data: dat1
## Control: lmerControl(optimizer = "Nelder_Mead")
##
## AIC BIC logLik deviance df.resid
## 4068.9 4115.9 -2024.5 4048.9 801
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2580 -0.6570 0.0492 0.6325 2.8096
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## groupid (Intercept) 12.36109 3.5158
## I(length - 15) 0.01178 0.1085 0.67
## Residual 5.98870 2.4472
## Number of obs: 811, groups: groupid, 100
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 43.72636 0.69806 62.640
## I(length - 15) 0.91655 0.02861 32.038
## goldcolor -0.49156 0.18333 -2.681
## I(volume - 1) 14.91728 1.97460 7.555
## I(nfish - 8) -1.09049 0.19743 -5.523
## I(length - 15):I(volume - 1) 0.19096 0.07989 2.390
##
## Correlation of Fixed Effects:
## (Intr) I(-15) gldclr I(v-1) I(n-8)
## I(lngth-15) 0.266
## goldcolor -0.143 -0.008
## I(volume-1) 0.841 0.165 -0.013
## I(nfish-8) -0.644 0.033 0.011 -0.736
## I(-15):I(-1 0.220 0.634 -0.046 0.297 -0.028Next, find the values of fish lengthes to probe the interaction. The values of 9, 14, and 19 are used because they are close to M - SD, M, and M + SD.
mean(dat1$length) - sd(dat1$length)## [1] 9.056527mean(dat1$length)## [1] 14.06846mean(dat1$length) + sd(dat1$length)## [1] 19.08039lengthval <- c(9, 14, 19)Using the same logic, the tank volume values of 0.5, 0.8, and 1.1 are used.
mean(dat1$volume) - sd(dat1$volume)## [1] 0.4918327mean(dat1$volume)## [1] 0.7783384mean(dat1$volume) + sd(dat1$volume)## [1] 1.064844volumeval <- c(0.5, 0.8, 1.1)Because the interactions package does not allow
predictors or moderators in I(), users need to transform
the values in the data frame or simply use the original scale. In this
case, the original scales of length and volume
are used. Note that nfish is still centered to make the
expected means meaningful.
o1m1 <- lmer(consume ~ 1 + goldcolor + I(nfish - 8)
+ length*volume + (1 + length|groupid),
data=dat1, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))Use the sim_slopes package to calculate simple slopes.
First, check the effect of fish tank volumes at each level of fish
lengths.
ss11 <- sim_slopes(model=o1m1, pred=volume, modx=length, modx.values = lengthval)ss11## JOHNSON-NEYMAN INTERVAL
##
## When length is OUTSIDE the interval [-391.49, -28.00], the slope of volume
## is p < .05.
##
## Note: The range of observed values of length is [3.17, 24.33]
##
## SIMPLE SLOPES ANALYSIS
##
## Slope of volume when length = 9.00:
##
## Est. S.E. t val. p
## ------- ------ -------- ------
## 13.77 1.93 7.13 0.00
##
## Slope of volume when length = 14.00:
##
## Est. S.E. t val. p
## ------- ------ -------- ------
## 14.73 2.00 7.38 0.00
##
## Slope of volume when length = 19.00:
##
## Est. S.E. t val. p
## ------- ------ -------- ------
## 15.68 2.14 7.34 0.00Check the interaction plot.
interact_plot(model=o1m1, pred=volume, modx=length, modx.values = lengthval)## Using data dat1 from global environment. This could cause incorrect results
## if dat1 has been altered since the model was fit. You can manually provide
## the data to the "data =" argument.
Next, check the simple slope effect of fish lengths at each level of fish tank volume.
ss12 <- sim_slopes(model=o1m1, pred=length, modx=volume, modx.values = volumeval)ss12## JOHNSON-NEYMAN INTERVAL
##
## When volume is OUTSIDE the interval [-24.81, -1.73], the slope of length is
## p < .05.
##
## Note: The range of observed values of volume is [0.04, 1.60]
##
## SIMPLE SLOPES ANALYSIS
##
## Slope of length when volume = 0.50:
##
## Est. S.E. t val. p
## ------ ------ -------- ------
## 0.82 0.03 26.01 0.00
##
## Slope of length when volume = 0.80:
##
## Est. S.E. t val. p
## ------ ------ -------- ------
## 0.88 0.02 38.90 0.00
##
## Slope of length when volume = 1.10:
##
## Est. S.E. t val. p
## ------ ------ -------- ------
## 0.94 0.03 26.87 0.00The interaction can be plotted as well.
interact_plot(model=out2m1b, pred=ersex, modx=eesex)## Using data dat2 from global environment. This could cause incorrect results
## if dat2 has been altered since the model was fit. You can manually provide
## the data to the "data =" argument.
Other Interaction Types
Example 3: Customers Satisfaction within Each Table in Dining Services
Read the data set from Lecture 4.
dat3 <-read.table("lecture4ex3.csv", sep=",", header=TRUE)Run the model with the interaction between age and sex.
out3m1 <- lmer(sat ~ 1 + I(age - 40)*female + I(numperson-4) + outdoor
+ (1|tableid), data=dat3, REML=FALSE)
summary(out3m1)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: sat ~ 1 + I(age - 40) * female + I(numperson - 4) + outdoor +
## (1 | tableid)
## Data: dat3
##
## AIC BIC logLik deviance df.resid
## 16541.1 16586.9 -8262.5 16525.1 2254
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.00219 -0.59995 -0.00421 0.61277 2.87979
##
## Random effects:
## Groups Name Variance Std.Dev.
## tableid (Intercept) 65.93 8.119
## Residual 59.60 7.720
## Number of obs: 2262, groups: tableid, 500
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 65.02384 0.62221 104.504
## I(age - 40) 0.24028 0.01816 13.234
## female -1.98218 0.37652 -5.264
## I(numperson - 4) -0.26536 0.19641 -1.351
## outdoor 3.01743 0.80863 3.732
## I(age - 40):female 0.03998 0.02504 1.597
##
## Correlation of Fixed Effects:
## (Intr) I(g-40) female I(n-4) outdor
## I(age - 40) 0.143
## female -0.322 -0.250
## I(nmprsn-4) -0.217 -0.002 0.001
## outdoor -0.647 0.008 0.000 -0.008
## I(g-40):fml -0.110 -0.730 0.278 -0.004 0.006The interaction was significant. Next, the interactions
package is used. Remove I() and rerun the analysis.
out3m1a <- lmer(sat ~ 1 + age*female + I(numperson-4) + outdoor
+ (1|tableid), data=dat3, REML=FALSE)Use the sim_slopes to calculate simple slopes. First,
check the effect of age at each sex.
ss31 <- sim_slopes(model=out3m1a, pred=age, modx=female, modx.values = c(0, 1))ss31## JOHNSON-NEYMAN INTERVAL
##
## When female is INSIDE the interval [-2.38, 29.20], the slope of age is p <
## .05.
##
## Note: The range of observed values of female is [0.00, 1.00]
##
## SIMPLE SLOPES ANALYSIS
##
## Slope of age when female = 0.00:
##
## Est. S.E. t val. p
## ------ ------ -------- ------
## 0.24 0.02 13.22 0.00
##
## Slope of age when female = 1.00:
##
## Est. S.E. t val. p
## ------ ------ -------- ------
## 0.28 0.02 16.36 0.00Check the interaction plot.
interact_plot(model=out3m1a, pred=age, modx=female, modx.values = c(0, 1))## Using data dat3 from global environment. This could cause incorrect results
## if dat3 has been altered since the model was fit. You can manually provide
## the data to the "data =" argument.
Next, check the sex differences at each level of age. Before probing
the interaction, find the appropriate values of age to probe. Let’s try
the quantile function to get 25th,
50th, and 75th percentiles. The values of 25, 35,
and 45 are picked because they are close to the quantile
result.
quantile(dat3$age, c(0.25, 0.50, 0.75))## 25% 50% 75%
## 25.00 36.00 46.75ageval <- c(25, 35, 45)Next, check the simple effect of sex difference at each level of age.
ss32 <- sim_slopes(model=out3m1a, pred=female, modx=age, modx.values = ageval)ss32## JOHNSON-NEYMAN INTERVAL
##
## When age is INSIDE the interval [-197.59, 57.54], the slope of female is p
## < .05.
##
## Note: The range of observed values of age is [1.00, 75.00]
##
## SIMPLE SLOPES ANALYSIS
##
## Slope of female when age = 25.00:
##
## Est. S.E. t val. p
## ------- ------ -------- ------
## -2.58 0.45 -5.71 0.00
##
## Slope of female when age = 35.00:
##
## Est. S.E. t val. p
## ------- ------ -------- ------
## -2.18 0.36 -6.02 0.00
##
## Slope of female when age = 45.00:
##
## Est. S.E. t val. p
## ------- ------ -------- ------
## -1.78 0.43 -4.15 0.00The interaction can be plotted as well.
interact_plot(model=out3m1a, pred=female, modx=age, modx.values = ageval)## Using data dat3 from global environment. This could cause incorrect results
## if dat3 has been altered since the model was fit. You can manually provide
## the data to the "data =" argument.
Let’s make the clustered bar graphs. First, the expected means of males and females at the ages of 25, 35, and 45 are calculated. First, the simple intercepts and slopes are needed.
ss32$ints## Value of age Est. S.E. 2.5% 97.5% t val. p
## 1 25 60.15991 1.031141 58.14564 62.17418 58.34305 1.022540e-214
## 2 35 62.76776 1.021448 60.77251 64.76301 61.44982 4.724007e-219
## 3 45 65.37562 1.026823 63.36982 67.38142 63.66786 4.814870e-228ss32$slopes## Value of age Est. S.E. 2.5% 97.5% t val. p
## 1 25 -2.581921 0.4525202 -3.467833 -1.696008 -5.705649 1.338296e-08
## 2 35 -2.182094 0.3627457 -2.892226 -1.471962 -6.015493 2.136203e-09
## 3 45 -1.782268 0.4289849 -2.622084 -0.942451 -4.154616 3.399653e-05The simple intercepts are the expected male satisfaction at each age. The sums of simple intercepts and simple slopes are the expected female satisfaction at each age. Next, a data frame containing expected means is created.
mydat <- data.frame(sex=c(rep("Male", 3), rep("Female", 3)),
age=rep(ageval, 2),
sat=c(ss32$ints[,"Est."], ss32$ints[,"Est."] + ss32$slopes[,"Est."]))
mydat## sex age sat
## 1 Male 25 60.15991
## 2 Male 35 62.76776
## 3 Male 45 65.37562
## 4 Female 25 57.57799
## 5 Female 35 60.58567
## 6 Female 45 63.59335Similar to the first example, the clustered bar graph is created by
the ggplot2 package.
library(ggplot2)
ggplot(mydat, aes(factor(age), sat, fill = sex)) +
geom_bar(stat="identity", position = "dodge") +
labs(x = "Age", y = "Customer Satisfaction", fill = "Sex")