Please, don’t change the structure of the document. Submit your knitted .html file no later that 11:59 p.m. May, 15 as Home task 2 project in LMS.

We have a dataset containing information about happiness, age, and gender on an individual level and gender-equality on country level. (hometask2.csv) The dataset was artificially generated for this exercise.

country - country #
age - scaled age
gender - 0 - female, 1 - male
g.eq - scaled gender equality index (country-level variable)
happy = level of happiness, scaled

We are interested to know how the level of individual happiness depends on age, gender, and gender equality.

Task 1

Make a plot showing the difference in happiness distribution across countries (1 point)

data<-read.csv("C:/Users/ASUS/Desktop/multilevel regression 2020/hometask2.csv")
install.packages("lme4")
## Installing package into 'C:/Users/ASUS/Documents/R/win-library/3.6'
## (as 'lib' is unspecified)
## package 'lme4' successfully unpacked and MD5 sums checked
## Warning: cannot remove prior installation of package 'lme4'
## Warning in file.copy(savedcopy, lib, recursive = TRUE): problem copying C:
## \Users\ASUS\Documents\R\win-library\3.6\00LOCK\lme4\libs\x64\lme4.dll to C:
## \Users\ASUS\Documents\R\win-library\3.6\lme4\libs\x64\lme4.dll: Permission
## denied
## Warning: restored 'lme4'
## 
## The downloaded binary packages are in
##  C:\Users\ASUS\AppData\Local\Temp\RtmpsB43qt\downloaded_packages
install.packages("sjPlot")
## Installing package into 'C:/Users/ASUS/Documents/R/win-library/3.6'
## (as 'lib' is unspecified)
## package 'sjPlot' successfully unpacked and MD5 sums checked
## 
## The downloaded binary packages are in
##  C:\Users\ASUS\AppData\Local\Temp\RtmpsB43qt\downloaded_packages
install.packages("ICC")
## Installing package into 'C:/Users/ASUS/Documents/R/win-library/3.6'
## (as 'lib' is unspecified)
## package 'ICC' successfully unpacked and MD5 sums checked
## 
## The downloaded binary packages are in
##  C:\Users\ASUS\AppData\Local\Temp\RtmpsB43qt\downloaded_packages
install.packages("ggplot2")
## Installing package into 'C:/Users/ASUS/Documents/R/win-library/3.6'
## (as 'lib' is unspecified)
## package 'ggplot2' successfully unpacked and MD5 sums checked
## 
## The downloaded binary packages are in
##  C:\Users\ASUS\AppData\Local\Temp\RtmpsB43qt\downloaded_packages
install.packages("glmmTMB")
## Installing package into 'C:/Users/ASUS/Documents/R/win-library/3.6'
## (as 'lib' is unspecified)
## package 'glmmTMB' successfully unpacked and MD5 sums checked
## Warning: cannot remove prior installation of package 'glmmTMB'
## Warning in file.copy(savedcopy, lib, recursive = TRUE): problem copying C:
## \Users\ASUS\Documents\R\win-library\3.6\00LOCK\glmmTMB\libs\x64\glmmTMB.dll
## to C:\Users\ASUS\Documents\R\win-library\3.6\glmmTMB\libs\x64\glmmTMB.dll:
## Permission denied
## Warning: restored 'glmmTMB'
## 
## The downloaded binary packages are in
##  C:\Users\ASUS\AppData\Local\Temp\RtmpsB43qt\downloaded_packages
install.packages("mitml")
## Installing package into 'C:/Users/ASUS/Documents/R/win-library/3.6'
## (as 'lib' is unspecified)
## package 'mitml' successfully unpacked and MD5 sums checked
## 
## The downloaded binary packages are in
##  C:\Users\ASUS\AppData\Local\Temp\RtmpsB43qt\downloaded_packages
library(glmmTMB)
## Warning: package 'glmmTMB' was built under R version 3.6.3
## Warning in checkMatrixPackageVersion(): Package version inconsistency detected.
## TMB was built with Matrix version 1.2.18
## Current Matrix version is 1.2.17
## Please re-install 'TMB' from source using install.packages('TMB', type = 'source') or ask CRAN for a binary version of 'TMB' matching CRAN's 'Matrix' package
library(lme4)
## Warning: package 'lme4' was built under R version 3.6.3
## Loading required package: Matrix
library(mitml)
## Warning: package 'mitml' was built under R version 3.6.3
## *** This is beta software. Please report any bugs!
## *** See the NEWS file for recent changes.
library(sjPlot)
## Warning: package 'sjPlot' was built under R version 3.6.3
## Install package "strengejacke" from GitHub (`devtools::install_github("strengejacke/strengejacke")`) to load all sj-packages at once!
library(ICC)
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 3.6.3

Task 2

Estimate an empty model with happy as DV and country as a grouping variable. Describe the result. What is the happiest country? (1 point)

moodel <- lmer(happy ~ (1 | country), data = data)
plot_model(moodel, type = "re", sort.est = "sort.all", grid = F)

# the graph shows that the 15 country is the most happiest in comparison with the rest of samples. ## Task 3 Calculate ICC for the empty model. Decide, whether multilevel modeling is necessary. (1 point)

ICCbare(country,happy,data)
## Warning in ICCbare(country, happy, data): 'x' has been coerced to a factor
## [1] 0.5487436
###[1] 0.5487436 the ICC is > than 0.1 thus multilevel regression modeling is necessary to fit the model.

Task 4

Add gender and age as fixed effects. Decide whether it is a good idea to include those predictors. Interpret fixed effects coefficients. (1 point)

model1<-lmer(happy~gender+(1|country),data=data)
model2<-lmer(happy~age+ (1|country),data = data)
anova(moodel,model1)
## refitting model(s) with ML (instead of REML)
## Data: data
## Models:
## moodel: happy ~ (1 | country)
## model1: happy ~ gender + (1 | country)
##        npar   AIC   BIC logLik deviance  Chisq Df Pr(>Chisq)    
## moodel    3 53077 53099 -26536    53071                         
## model1    4 52652 52680 -26322    52644 427.83  1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(moodel,model1)
## refitting model(s) with ML (instead of REML)
## Data: data
## Models:
## moodel: happy ~ (1 | country)
## model1: happy ~ gender + (1 | country)
##        npar   AIC   BIC logLik deviance  Chisq Df Pr(>Chisq)    
## moodel    3 53077 53099 -26536    53071                         
## model1    4 52652 52680 -26322    52644 427.83  1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# ANOVA test shows the P-value of moodel and comparing with model1- model2 is statistically significant,and the AIC and the BIC for model1 and model2 are small it means that the quality of each variable is good and   both can be include in the model  
model3<-lmer(happy~ gender + age + (1|country),data=data)
anova(model3,moodel)
## refitting model(s) with ML (instead of REML)
## Data: data
## Models:
## moodel: happy ~ (1 | country)
## model3: happy ~ gender + age + (1 | country)
##        npar   AIC   BIC   logLik deviance Chisq Df Pr(>Chisq)    
## moodel    3 53077 53099 -26535.7    53071                        
## model3    5 16873 16910  -8431.7    16863 36208  2  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(model3,model1)
## refitting model(s) with ML (instead of REML)
## Data: data
## Models:
## model1: happy ~ gender + (1 | country)
## model3: happy ~ gender + age + (1 | country)
##        npar   AIC   BIC   logLik deviance Chisq Df Pr(>Chisq)    
## model1    4 52652 52680 -26321.8    52644                        
## model3    5 16873 16910  -8431.7    16863 35780  1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(model3,model2)
## refitting model(s) with ML (instead of REML)
## Data: data
## Models:
## model2: happy ~ age + (1 | country)
## model3: happy ~ gender + age + (1 | country)
##        npar   AIC   BIC   logLik deviance  Chisq Df Pr(>Chisq)    
## model2    4 25813 25842 -12902.5    25805                         
## model3    5 16874 16910  -8431.7    16864 8941.6  1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot_model(model3, type = "re", grid = T)

summary(model3)
## Linear mixed model fit by REML ['lmerMod']
## Formula: happy ~ gender + age + (1 | country)
##    Data: data
## 
## REML criterion at convergence: 16878.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.6648 -0.6380  0.0063  0.6371  5.6805 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  country  (Intercept) 1.7609   1.3270  
##  Residual             0.3053   0.5525  
## Number of obs: 10000, groups:  country, 50
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 0.113409   0.187758   0.604
## gender      3.112471   0.025860 120.358
## age         5.002661   0.008465 590.983
## 
## Correlation of Fixed Effects:
##        (Intr) gender
## gender -0.007       
## age     0.010 -0.009
# the general pattern for the sample shows that men have 3.046025 times more happiness than women

Task 5

Randomize the effect of gender, decide whether it is a good idea to randomize it. What is the country with the strongest relationship between gender and happiness? (1 point)

model4<-lmer(happy ~ gender + age + ( 1 + gender | country), data = data)
summary(model4)
## Linear mixed model fit by REML ['lmerMod']
## Formula: happy ~ gender + age + (1 + gender | country)
##    Data: data
## 
## REML criterion at convergence: 15080.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.1621 -0.6738  0.0066  0.6732  3.3843 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr
##  country  (Intercept) 1.6672   1.2912       
##           gender      1.2083   1.0992   0.59
##  Residual             0.2507   0.5007       
## Number of obs: 10000, groups:  country, 50
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 0.113612   0.182684   0.622
## gender      3.045882   0.157376  19.354
## age         5.004128   0.007693 650.475
## 
## Correlation of Fixed Effects:
##        (Intr) gender
## gender  0.585       
## age     0.010 -0.001
ranef(model4)
## $country
##    (Intercept)      gender
## 1  -1.68658601 -0.07449271
## 2   0.24967652  0.44601154
## 3   1.03500165  0.53109340
## 4  -1.06296966 -1.31083016
## 5   1.73777350  0.84135572
## 6   2.31998777  1.97340155
## 7  -2.66337820 -1.23285059
## 8  -2.10180112 -1.17995499
## 9  -0.43674386 -0.01665233
## 10  0.21490861 -0.20137466
## 11  3.68600387  1.61096290
## 12 -0.05530219 -0.95392938
## 13  0.64624802 -0.34445153
## 14  0.86935708  2.61333036
## 15  1.76849986  1.00341605
## 16 -1.25939373 -1.34918842
## 17 -1.67360352  0.66953774
## 18 -3.24971760 -2.03091308
## 19  0.81638441  0.81197612
## 20  1.16753603  0.56527005
## 21 -0.91719750 -1.83526370
## 22  0.46468258 -0.62348734
## 23 -1.14019157 -1.48634259
## 24  1.97755623  1.03943098
## 25  1.02478555  0.77980511
## 26  1.07767045  1.29789962
## 27 -0.45940534  0.23761986
## 28 -0.03316577  0.28034836
## 29  0.66329818 -1.54174965
## 30  0.90464021  1.85781281
## 31 -1.35929689  0.06464082
## 32 -0.80230619 -0.63416163
## 33  0.55121549 -0.88126003
## 34 -1.15489413  0.17918365
## 35 -0.10433178 -1.27251552
## 36 -0.31688958 -0.01524034
## 37 -0.87723667 -1.20330606
## 38  1.06229022  1.06719060
## 39  0.59931275 -0.17800674
## 40 -0.90939754  0.86470738
## 41 -0.23107571  0.71092788
## 42 -0.20109613 -1.01352238
## 43  0.41470055  0.06641976
## 44  1.44852210  0.58661363
## 45 -0.87112000 -0.58946807
## 46  0.06868160  0.14787867
## 47 -0.26758685  0.93593707
## 48 -0.82127426  0.14884527
## 49  0.35938601 -2.02747707
## 50 -0.47215746  0.66482209
## 
## with conditional variances for "country"
# the random effect shows the correlation between gender and happiness like the country number 7,8 and 18

##Task 6 Visualize and briefly describe the distribution of random effects (1 point)

plot_model(model3,type = "re", grid = T)

##Task 7 Add gender equality index as a group-level predictor. Interpret the results. (1 point)

model5<-lmer(happy~ gender + age + g.eq+(1|country),data = data)
summary(model5)
## Linear mixed model fit by REML ['lmerMod']
## Formula: happy ~ gender + age + g.eq + (1 | country)
##    Data: data
## 
## REML criterion at convergence: 16855.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.6609 -0.6387  0.0056  0.6379  5.6767 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  country  (Intercept) 1.0898   1.0439  
##  Residual             0.3053   0.5525  
## Number of obs: 10000, groups:  country, 50
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 0.051061   0.148179   0.345
## gender      3.112440   0.025860 120.357
## age         5.002881   0.008462 591.185
## g.eq        0.930071   0.166698   5.579
## 
## Correlation of Fixed Effects:
##        (Intr) gender age   
## gender -0.008              
## age     0.013 -0.009       
## g.eq   -0.075 -0.001  0.008

##Task 8 We are interested to know how the relationship between gender and happiness depends on the level of gender equality. Add the cross-level interaction effect to learn that. Visualize the interaction, interpret the results. (1 point)

model6<-lmer(happy ~age + g.eq*gender + ( 1 + gender | country), data = data)
summary(model6,detail=T, digits=2)
## Warning in summary.merMod(model6, detail = T, digits = 2): additional arguments
## ignored
## Linear mixed model fit by REML ['lmerMod']
## Formula: happy ~ age + g.eq * gender + (1 + gender | country)
##    Data: data
## 
## REML criterion at convergence: 14912.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.1616 -0.6733  0.0067  0.6700  3.4006 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  country  (Intercept) 1.081520 1.03996       
##           gender      0.004012 0.06334  -0.12
##  Residual             0.250873 0.50087       
## Number of obs: 10000, groups:  country, 50
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 0.055420   0.147593   0.375
## age         5.005126   0.007676 652.079
## g.eq        0.870945   0.166046   5.245
## gender      2.973543   0.025398 117.079
## g.eq:gender 1.220313   0.028257  43.186
## 
## Correlation of Fixed Effects:
##             (Intr) age    g.eq   gender
## age          0.011                     
## g.eq        -0.075  0.007              
## gender      -0.049 -0.009  0.004       
## g.eq:gender  0.004  0.006 -0.050 -0.120
plot_model(model6,type = "pred",terms = c("g.eq","gender"))

##Task 9 Write down a system of equations for the last model (with i’s and j’s specified correctly). You can use latex equations or insert an image. (1 point, you can earn an extra point here!)

##Task 10 Create a table combining outputs of all estimated models (1 point)

tab_model(moodel,model1, model2, model3, model4, model5, model6,  show.ci = F)
  happy happy happy happy happy happy happy
Predictors Estimates p Estimates p Estimates p Estimates p Estimates p Estimates p Estimates p
(Intercept) -0.89 0.094 -1.04 0.049 0.27 0.161 0.11 0.546 0.11 0.534 0.05 0.730 0.06 0.707
gender 3.25 <0.001 3.11 <0.001 3.05 <0.001 3.11 <0.001 2.97 <0.001
age 5.01 <0.001 5.00 <0.001 5.00 <0.001 5.00 <0.001 5.01 <0.001
g.eq 0.93 <0.001 0.87 <0.001
g.eq * gender 1.22 <0.001
Random Effects
σ2 11.49 11.01 0.75 0.31 0.25 0.31 0.25
τ00 13.98 country 13.93 country 1.79 country 1.76 country 1.67 country 1.09 country 1.08 country
τ11         1.21 country.gender   0.00 country.gender
ρ01         0.59 country   -0.12 country
ICC 0.55 0.56 0.70 0.85 0.88 0.78 0.81
N 50 country 50 country 50 country 50 country 50 country 50 country 50 country
Observations 10000 10000 10000 10000 10000 10000 10000
Marginal R2 / Conditional R2 0.000 / 0.549 0.019 / 0.567 0.904 / 0.972 0.922 / 0.988 0.922 / 0.991 0.945 / 0.988 0.947 / 0.990