library(readxl)
## Warning: package 'readxl' was built under R version 3.6.3
library(nlme)
library(lattice)
library(lme4)
## Warning: package 'lme4' was built under R version 3.6.3
## Loading required package: Matrix
## 
## Attaching package: 'lme4'
## The following object is masked from 'package:nlme':
## 
##     lmList
OurData<-read.csv("C:\\Users\\Owner\\GLMM\\HapData.csv")
attach(OurData)
#Entry Analysis
hist(OurData$HappinessScore,xlab = "Happiness Score", ylab = "Number of Countries",main = "Frequency of Happiness Scores")

use a lmm because normal distribution.

BeerMod <- lmer(HappinessScore~Beer_PerCapita + (1|Region), REML = FALSE)
WineMod <- lmer(HappinessScore~Wine_PerCapita +Beer_PerCapita + (1|Region), REML = FALSE)
SpiritMod <- lmer(HappinessScore~Spirit_PerCapita+Beer_PerCapita + (1|Region), REML = FALSE)
anova(BeerMod,SpiritMod,test="Chisq")[8]
##           Pr(>Chisq)
## BeerMod             
## SpiritMod     0.4519
anova(BeerMod,WineMod,test="Chisq")[8]
##         Pr(>Chisq)
## BeerMod           
## WineMod     0.8667

p values are high so beer mod is better

Finalmod<-lmer(HappinessScore~Beer_PerCapita + (1|Region))
summary(Finalmod)
## Linear mixed model fit by REML ['lmerMod']
## Formula: HappinessScore ~ Beer_PerCapita + (1 | Region)
## 
## REML criterion at convergence: 290.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3747 -0.5512  0.0171  0.6884  2.4861 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Region   (Intercept) 0.6431   0.8019  
##  Residual             0.4762   0.6901  
## Number of obs: 122, groups:  Region, 9
## 
## Fixed effects:
##                Estimate Std. Error t value
## (Intercept)    5.482686   0.301682  18.174
## Beer_PerCapita 0.002645   0.000730   3.623
## 
## Correlation of Fixed Effects:
##             (Intr)
## Beer_PerCpt -0.340
nullmod<-lm(HappinessScore~Beer_PerCapita,data=OurData)
teststat<-2*as.numeric(logLik(Finalmod)-logLik(nullmod))
pchisq(teststat,df=1,lower.tail = FALSE)
## [1] 1.44288e-13

pvalue is below sig level of .05 meaning it is better to use lmm than lm.

randslope <-lmer(HappinessScore~Beer_PerCapita+(Beer_PerCapita|Region))
## boundary (singular) fit: see ?isSingular
teststat2<-2*as.numeric(logLik(randslope)-logLik(Finalmod))
pchisq(teststat2,df=1,lower.tail = FALSE)
## [1] 0.09342283

this proves what that paragraph we have says. p val is wrong on here for soe reason should be .09

 confint(Finalmod,level=.95)
## Computing profile confidence intervals ...
##                      2.5 %      97.5 %
## .sig01         0.452946052 1.362866550
## .sigma         0.606457135 0.787852692
## (Intercept)    4.888299499 6.122649181
## Beer_PerCapita 0.001224304 0.004146477

We are 95% confident that an increase of 1 beer per year per capita leads to a multiplicative change in happiness between 0.0012243 and .004164

xyplot(HappinessScore~Wine_PerCapita|Region,ylab = "Happiness Score", xlab = "Wine Consumption")

xyplot(HappinessScore~Spirit_PerCapita|Region,ylab = "Happiness Score", xlab = "Spirit Consumption")

xyplot(HappinessScore~Beer_PerCapita|Region,ylab = "Happiness Score", xlab = "Beer Consumption")

here are the xyplots