2017meathappy
Alice Melvin
16/10/2020
New plots using the new meat consumption data. Has 140 data points.
#this is to set a working directory
setwd("C:\\Users\\melvo\\Documents\\CCMF\\Quantitive methods\\Coursework\\Meat and Happiness")
#this checks the working directory
getwd()## [1] "C:/Users/melvo/Documents/CCMF/Quantitive methods/Coursework/Meat and Happiness"#selecting your dataset
meathap <- read.csv("meathappy_plot2017.csv")Here comes a graph!
library(ggplot2)
ggplot(meathap, aes(x = meat, y = happiness)) +
geom_smooth(method = "lm", se = FALSE) +
geom_point()+theme_minimal()+
xlab("Meat Consumption kg/cap")+ylab("Happiness Score in Cantril Ladder") + ggtitle("Relationship between meat consumption and happiness 2017")## `geom_smooth()` using formula 'y ~ x'
reg1=lm(happiness~meat,meathap)
summary(reg1)##
## Call:
## lm(formula = happiness ~ meat, data = meathap)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.3124 -0.6721 0.1232 0.7356 2.6435
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.731108 0.168088 28.147 < 2e-16 ***
## meat 0.015068 0.002953 5.103 1.09e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.046 on 138 degrees of freedom
## Multiple R-squared: 0.1587, Adjusted R-squared: 0.1526
## F-statistic: 26.04 on 1 and 138 DF, p-value: 1.086e-06 coeftest(reg1, vcov=vcovHC)##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.7311081 0.1580845 29.9277 < 2.2e-16 ***
## meat 0.0150685 0.0029554 5.0987 1.107e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Note: according to the internet vcov is a specification of the covariance matrix of the estimated coefficients. Here by specifying vcovHC a function for extracting the covariance matrix from x is supplied from the package sandwich- The HC is the heteroskedasticity-consistent covariance matrix estimation of the covariance matrix of the coefficient estimates in a linear regression model. So this is what we use for a fitted model in the class lm as we have here. btw heteroscedasticity refers to the circumstance is which the variability of a variable is unequal across the range of values of a second variable that predicts it.
Also, the fact that the relationship between meat and happiness seems to disappear when controlling for gdp per capita suggests an endogeneity problem in the original model where there was a relationship between the x variable and the error term in the model. (I think this is right?) Our confounding variable is therefore gdp per capita. We can also look into other potential confounding variables? Draw arrow diagram to show this? Another suggestion is the urbanisation of the country, more urbanised meaning more meat consumption and potentially having an effect on perceived happiness too? Also meat consumption is suggested to increase in western countries, so could run a west vs east comparison. More women in the workforce could also mean more meat consumption “because meat is a “convenient food” insofar as it is readily available as a mainstay of fast-food menus" and would the fact that their working make them happier or unhappier could be either really? Then we essentially come to the conclusion that once you take into account anything confounding which is affecting the happiness as well as the meat consumption, then the initial findings are void and actually happiness is not affected by the meat you eat.
library(ggplot2)
ggplot(meathap, aes(x = gdppc, y = happiness)) +
geom_smooth(method = "lm", se = FALSE) +
geom_point()+theme_minimal()+
xlab("GDP/cap USD")+ylab("Happiness Score") + ggtitle("Relationship between GDP per capita and happiness 2017")## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 1 rows containing non-finite values (stat_smooth).## Warning: Removed 1 rows containing missing values (geom_point).
reg2=lm(happiness~gdppc,meathap)
summary(reg2)##
## Call:
## lm(formula = happiness ~ gdppc, data = meathap)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.655 -0.775 0.171 0.711 2.422
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.136e+00 1.111e-01 46.224 < 2e-16 ***
## gdppc 2.241e-05 4.524e-06 4.953 2.12e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.054 on 137 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.1519, Adjusted R-squared: 0.1457
## F-statistic: 24.53 on 1 and 137 DF, p-value: 2.117e-06 coeftest(reg2, vcov=vcovHC)##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.1365e+00 1.0939e-01 46.954 < 2.2e-16 ***
## gdppc 2.2409e-05 5.8678e-06 3.819 0.0002023 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1By the way, the warning about the “missing rows” is because there was no GDP/cap for Venezuela.
reg3=lm(happiness~meat+gdppc,meathap)
summary(reg3)##
## Call:
## lm(formula = happiness ~ meat + gdppc, data = meathap)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.4144 -0.6655 0.1242 0.7549 2.6207
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.819e+00 1.713e-01 28.137 <2e-16 ***
## meat 9.422e-03 3.910e-03 2.410 0.0173 *
## gdppc 1.290e-05 5.946e-06 2.169 0.0318 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.036 on 136 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.1866, Adjusted R-squared: 0.1746
## F-statistic: 15.6 on 2 and 136 DF, p-value: 7.954e-07 coeftest(reg3, vcov=vcovHC)##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.8186e+00 1.6150e-01 29.8368 < 2e-16 ***
## meat 9.4217e-03 3.6801e-03 2.5602 0.01155 *
## gdppc 1.2899e-05 6.9822e-06 1.8474 0.06686 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Now let’s look at each of the following 6 regions: Europe, Africa, Asia & Pacific, South America, North America, Middle East
eumeathap <- read.csv("Eu_meathappy_plot2017.csv")
library(ggplot2)
ggplot(eumeathap, aes(x = meat, y = happiness)) +
geom_smooth(method = "lm", se = FALSE) +
geom_point()+theme_minimal()+
xlab("Meat Consumption kg/cap")+ylab("Happiness Score in Cantril Ladder") + ggtitle("Relationship between meat consumption and happiness Europe 2017")## `geom_smooth()` using formula 'y ~ x'
library(ggplot2)
ggplot(eumeathap, aes(x = gdppc, y = happiness)) +
geom_smooth(method = "lm", se = FALSE) +
geom_point()+theme_minimal()+
xlab("GDP/cap USD")+ylab("Happiness Score") + ggtitle("Relationship between GDP per capita and happiness Europe 2017")## `geom_smooth()` using formula 'y ~ x'
reg4=lm(happiness~meat,eumeathap)
summary(reg4)##
## Call:
## lm(formula = happiness ~ meat, data = eumeathap)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.96615 -1.03228 0.09871 1.06116 2.02390
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.628080 0.740911 7.596 2.07e-09 ***
## meat 0.001641 0.010429 0.157 0.876
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.216 on 42 degrees of freedom
## Multiple R-squared: 0.0005892, Adjusted R-squared: -0.02321
## F-statistic: 0.02476 on 1 and 42 DF, p-value: 0.8757 coeftest(reg4, vcov=vcovHC)##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.628080 0.668496 8.419 1.467e-10 ***
## meat 0.001641 0.010005 0.164 0.8705
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 reg5=lm(happiness~meat+gdppc,eumeathap)
summary(reg5)##
## Call:
## lm(formula = happiness ~ meat + gdppc, data = eumeathap)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.3738 -0.9365 0.3635 0.8088 1.9966
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.957e+00 7.178e-01 8.299 2.59e-10 ***
## meat -1.075e-02 1.123e-02 -0.957 0.3439
## gdppc 1.875e-05 7.982e-06 2.349 0.0237 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.155 on 41 degrees of freedom
## Multiple R-squared: 0.1191, Adjusted R-squared: 0.07616
## F-statistic: 2.772 on 2 and 41 DF, p-value: 0.07425 coeftest(reg5, vcov=vcovHC)##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.9565e+00 6.3816e-01 9.3340 1.066e-11 ***
## meat -1.0748e-02 1.0513e-02 -1.0224 0.31260
## gdppc 1.8748e-05 8.6048e-06 2.1788 0.03515 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1afmeathap <- read.csv("Af_meathappy_plot2017.csv")
library(ggplot2)
ggplot(afmeathap, aes(x = meat, y = happiness)) +
geom_smooth(method = "lm", se = FALSE) +
geom_point()+theme_minimal()+
xlab("Meat Consumption kg/cap")+ylab("Happiness Score in Cantril Ladder") + ggtitle("Relationship between meat consumption and happiness Africa 2017")## `geom_smooth()` using formula 'y ~ x'
library(ggplot2)
ggplot(afmeathap, aes(x = gdppc, y = happiness)) +
geom_smooth(method = "lm", se = FALSE) +
geom_point()+theme_minimal()+
xlab("GDP/cap USD")+ylab("Happiness Score") + ggtitle("Relationship between GDP per capita and happiness Africa 2017")## `geom_smooth()` using formula 'y ~ x'
reg6=lm(happiness~meat,afmeathap)
summary(reg6)##
## Call:
## lm(formula = happiness ~ meat, data = afmeathap)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.51299 -0.53529 -0.09317 0.63643 1.64525
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.554581 0.258165 17.642 <2e-16 ***
## meat 0.007355 0.010335 0.712 0.481
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8666 on 36 degrees of freedom
## Multiple R-squared: 0.01387, Adjusted R-squared: -0.01352
## F-statistic: 0.5064 on 1 and 36 DF, p-value: 0.4813 coeftest(reg6, vcov=vcovHC)##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.5545807 0.2947824 15.4507 <2e-16 ***
## meat 0.0073545 0.0143187 0.5136 0.6106
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 reg7=lm(happiness~meat+gdppc,afmeathap)
summary(reg7)##
## Call:
## lm(formula = happiness ~ meat + gdppc, data = afmeathap)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.51029 -0.60861 -0.01451 0.60179 1.67165
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.536e+00 2.628e-01 17.260 <2e-16 ***
## meat 1.315e-02 1.467e-02 0.897 0.376
## gdppc -4.854e-05 8.625e-05 -0.563 0.577
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.875 on 35 degrees of freedom
## Multiple R-squared: 0.02272, Adjusted R-squared: -0.03313
## F-statistic: 0.4068 on 2 and 35 DF, p-value: 0.6689 coeftest(reg7, vcov=vcovHC)##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.5358e+00 3.0662e-01 14.7929 <2e-16 ***
## meat 1.3155e-02 2.0539e-02 0.6405 0.5260
## gdppc -4.8542e-05 1.0965e-04 -0.4427 0.6607
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1apacmeathap <- read.csv("Apac_meathappy_plot2017.csv")
library(ggplot2)
ggplot(apacmeathap, aes(x = meat, y = happiness)) +
geom_smooth(method = "lm", se = FALSE) +
geom_point()+theme_minimal()+
xlab("Meat Consumption kg/cap")+ylab("Happiness Score in Cantril Ladder") + ggtitle("Relationship between meat consumption and happiness Asia & Pacific 2017")## `geom_smooth()` using formula 'y ~ x'
library(ggplot2)
ggplot(apacmeathap, aes(x = gdppc, y = happiness)) +
geom_smooth(method = "lm", se = FALSE) +
geom_point()+theme_minimal()+
xlab("GDP/cap USD")+ylab("Happiness Score") + ggtitle("Relationship between GDP per capita and happiness Asia & Pacific 2017")## `geom_smooth()` using formula 'y ~ x'
reg8=lm(happiness~meat,apacmeathap)
summary(reg8)##
## Call:
## lm(formula = happiness ~ meat, data = apacmeathap)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.12547 -0.74134 -0.06837 0.54105 2.14199
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.371659 0.357754 15.015 7.53e-11 ***
## meat 0.004804 0.005591 0.859 0.403
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9482 on 16 degrees of freedom
## Multiple R-squared: 0.04411, Adjusted R-squared: -0.01563
## F-statistic: 0.7383 on 1 and 16 DF, p-value: 0.4029 coeftest(reg8, vcov=vcovHC)##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.3716588 0.4085543 13.1480 5.417e-10 ***
## meat 0.0048041 0.0065737 0.7308 0.4755
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 reg9=lm(happiness~meat+gdppc,apacmeathap)
summary(reg9)##
## Call:
## lm(formula = happiness ~ meat + gdppc, data = apacmeathap)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.0296 -0.7240 -0.1398 0.4913 2.1368
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.322e+00 3.701e-01 14.379 3.52e-10 ***
## meat 9.816e-03 9.058e-03 1.084 0.296
## gdppc -1.415e-05 1.992e-05 -0.710 0.488
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9632 on 15 degrees of freedom
## Multiple R-squared: 0.07522, Adjusted R-squared: -0.04808
## F-statistic: 0.61 on 2 and 15 DF, p-value: 0.5563 coeftest(reg9, vcov=vcovHC)##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.3219e+00 4.0452e-01 13.1559 1.218e-09 ***
## meat 9.8161e-03 7.9446e-03 1.2356 0.2356
## gdppc -1.4152e-05 2.1526e-05 -0.6574 0.5209
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1memeathap <- read.csv("ME_meathappy_plot2017.csv")
library(ggplot2)
ggplot(memeathap, aes(x = meat, y = happiness)) +
geom_smooth(method = "lm", se = FALSE) +
geom_point()+theme_minimal()+
xlab("Meat Consumption kg/cap")+ylab("Happiness Score in Cantril Ladder") + ggtitle("Relationship between meat consumption and happiness Middle East 2017")## `geom_smooth()` using formula 'y ~ x'
library(ggplot2)
ggplot(memeathap, aes(x = gdppc, y = happiness)) +
geom_smooth(method = "lm", se = FALSE) +
geom_point()+theme_minimal()+
xlab("GDP/cap USD")+ylab("Happiness Score") + ggtitle("Relationship between GDP per capita and happiness Middle East 2017")## `geom_smooth()` using formula 'y ~ x'
reg10=lm(happiness~meat,memeathap)
summary(reg10)##
## Call:
## lm(formula = happiness ~ meat, data = memeathap)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.8239 -0.7293 -0.1477 1.0015 1.5841
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.17478 0.52725 7.918 9.77e-07 ***
## meat 0.03619 0.01129 3.206 0.00589 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.08 on 15 degrees of freedom
## Multiple R-squared: 0.4066, Adjusted R-squared: 0.3671
## F-statistic: 10.28 on 1 and 15 DF, p-value: 0.005889 coeftest(reg10, vcov=vcovHC)##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.174777 0.639044 6.5329 9.475e-06 ***
## meat 0.036186 0.011230 3.2223 0.005698 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 reg11=lm(happiness~meat+gdppc,memeathap)
summary(reg11)##
## Call:
## lm(formula = happiness ~ meat + gdppc, data = memeathap)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.8911 -0.7119 -0.2116 1.0270 1.6443
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.299e+00 6.027e-01 7.133 5.07e-06 ***
## meat 2.850e-02 2.008e-02 1.420 0.178
## gdppc 1.659e-05 3.540e-05 0.469 0.647
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.109 on 14 degrees of freedom
## Multiple R-squared: 0.4158, Adjusted R-squared: 0.3323
## F-statistic: 4.982 on 2 and 14 DF, p-value: 0.02323 coeftest(reg11, vcov=vcovHC)##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.2987e+00 7.3372e-01 5.8588 4.155e-05 ***
## meat 2.8503e-02 2.0002e-02 1.4251 0.1761
## gdppc 1.6590e-05 2.8616e-05 0.5798 0.5713
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1nameathap <- read.csv("NA_meathappy_plot2017.csv")
library(ggplot2)
ggplot(nameathap, aes(x = meat, y = happiness)) +
geom_smooth(method = "lm", se = FALSE) +
geom_point()+theme_minimal()+
xlab("Meat Consumption kg/cap")+ylab("Happiness Score in Cantril Ladder") + ggtitle("Relationship between meat consumption and happiness North America 2017")## `geom_smooth()` using formula 'y ~ x'
library(ggplot2)
ggplot(nameathap, aes(x = gdppc, y = happiness)) +
geom_smooth(method = "lm", se = FALSE) +
geom_point()+theme_minimal()+
xlab("GDP/cap USD")+ylab("Happiness Score") + ggtitle("Relationship between GDP per capita and happiness North America 2017")## `geom_smooth()` using formula 'y ~ x'
reg12=lm(happiness~meat,nameathap)
summary(reg12)##
## Call:
## lm(formula = happiness ~ meat, data = nameathap)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.00832 -0.17925 0.05381 0.31814 1.02320
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.973233 0.514993 9.657 2.19e-06 ***
## meat 0.017168 0.008288 2.071 0.0651 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8012 on 10 degrees of freedom
## Multiple R-squared: 0.3002, Adjusted R-squared: 0.2303
## F-statistic: 4.291 on 1 and 10 DF, p-value: 0.06513 coeftest(reg12, vcov=vcovHC)##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.9732327 0.2996220 16.5984 1.317e-08 ***
## meat 0.0171682 0.0055253 3.1072 0.01111 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 reg13=lm(happiness~meat+gdppc,nameathap)
summary(reg13)##
## Call:
## lm(formula = happiness ~ meat + gdppc, data = nameathap)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.44668 -0.36182 0.02285 0.42873 0.81426
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.887e+00 6.836e-01 8.611 1.22e-05 ***
## meat -1.296e-02 1.816e-02 -0.714 0.494
## gdppc 5.162e-05 2.836e-05 1.820 0.102
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.722 on 9 degrees of freedom
## Multiple R-squared: 0.4885, Adjusted R-squared: 0.3749
## F-statistic: 4.298 on 2 and 9 DF, p-value: 0.04894 coeftest(reg13, vcov=vcovHC)##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.8869e+00 7.9511e-01 7.4039 4.086e-05 ***
## meat -1.2957e-02 2.7481e-02 -0.4715 0.6485
## gdppc 5.1618e-05 4.3506e-05 1.1865 0.2658
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1sameathap <- read.csv("SA_meathappy_plot2017.csv")
library(ggplot2)
ggplot(sameathap, aes(x = meat, y = happiness)) +
geom_smooth(method = "lm", se = FALSE) +
geom_point()+theme_minimal()+
xlab("Meat Consumption kg/cap")+ylab("Happiness Score in Cantril Ladder") + ggtitle("Relationship between meat consumption and happiness South America 2017")## `geom_smooth()` using formula 'y ~ x'
library(ggplot2)
ggplot(sameathap, aes(x = gdppc, y = happiness)) +
geom_smooth(method = "lm", se = FALSE) +
geom_point()+theme_minimal()+
xlab("GDP/cap USD")+ylab("Happiness Score") + ggtitle("Relationship between GDP per capita and happiness South America 2017")## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 1 rows containing non-finite values (stat_smooth).## Warning: Removed 1 rows containing missing values (geom_point).
reg14=lm(happiness~meat,sameathap)
summary(reg14)##
## Call:
## lm(formula = happiness ~ meat, data = sameathap)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.9476 -0.4039 -0.1015 0.2032 1.7681
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.320495 0.638820 8.329 1.6e-05 ***
## meat 0.009144 0.009411 0.972 0.357
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7584 on 9 degrees of freedom
## Multiple R-squared: 0.09493, Adjusted R-squared: -0.005632
## F-statistic: 0.944 on 1 and 9 DF, p-value: 0.3566 coeftest(reg14, vcov=vcovHC)##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.3204951 0.5211547 10.2091 3.011e-06 ***
## meat 0.0091436 0.0056990 1.6044 0.1431
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 reg15=lm(happiness~meat+gdppc,sameathap)
summary(reg15)##
## Call:
## lm(formula = happiness ~ meat + gdppc, data = sameathap)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.92606 -0.24818 -0.05533 -0.00795 1.71970
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.362e+00 7.441e-01 7.206 0.000177 ***
## meat 1.977e-03 1.147e-02 0.172 0.867971
## gdppc 5.299e-05 6.228e-05 0.851 0.423006
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7807 on 7 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.1403, Adjusted R-squared: -0.1054
## F-statistic: 0.571 on 2 and 7 DF, p-value: 0.5892 coeftest(reg15, vcov=vcovHC)##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.3616e+00 6.3250e-01 8.4768 6.287e-05 ***
## meat 1.9774e-03 5.5886e-03 0.3538 0.73389
## gdppc 5.2990e-05 2.1194e-05 2.5002 0.04098 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1North and South america havn’t got many data points so im going to merge them into total america (N+s)
ameathap <- read.csv("A_meathappy_plot2017.csv")
library(ggplot2)
ggplot(ameathap, aes(x = meat, y = happiness)) +
geom_smooth(method = "lm", se = FALSE) +
geom_point()+theme_minimal()+
xlab("Meat Consumption kg/cap")+ylab("Happiness Score in Cantril Ladder") + ggtitle("Relationship between meat consumption and happiness North & South America 2017")## `geom_smooth()` using formula 'y ~ x'
library(ggplot2)
ggplot(ameathap, aes(x = gdppc, y = happiness)) +
geom_smooth(method = "lm", se = FALSE) +
geom_point()+theme_minimal()+
xlab("GDP/cap USD")+ylab("Happiness Score") + ggtitle("Relationship between GDP per capita and happiness North & South America 2017")## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 1 rows containing non-finite values (stat_smooth).## Warning: Removed 1 rows containing missing values (geom_point).
reg16=lm(happiness~meat,sameathap)
summary(reg16)##
## Call:
## lm(formula = happiness ~ meat, data = sameathap)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.9476 -0.4039 -0.1015 0.2032 1.7681
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.320495 0.638820 8.329 1.6e-05 ***
## meat 0.009144 0.009411 0.972 0.357
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7584 on 9 degrees of freedom
## Multiple R-squared: 0.09493, Adjusted R-squared: -0.005632
## F-statistic: 0.944 on 1 and 9 DF, p-value: 0.3566 coeftest(reg16, vcov=vcovHC)##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.3204951 0.5211547 10.2091 3.011e-06 ***
## meat 0.0091436 0.0056990 1.6044 0.1431
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 reg17=lm(happiness~meat+gdppc,sameathap)
summary(reg17)##
## Call:
## lm(formula = happiness ~ meat + gdppc, data = sameathap)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.92606 -0.24818 -0.05533 -0.00795 1.71970
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.362e+00 7.441e-01 7.206 0.000177 ***
## meat 1.977e-03 1.147e-02 0.172 0.867971
## gdppc 5.299e-05 6.228e-05 0.851 0.423006
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7807 on 7 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.1403, Adjusted R-squared: -0.1054
## F-statistic: 0.571 on 2 and 7 DF, p-value: 0.5892 coeftest(reg17, vcov=vcovHC)##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.3616e+00 6.3250e-01 8.4768 6.287e-05 ***
## meat 1.9774e-03 5.5886e-03 0.3538 0.73389
## gdppc 5.2990e-05 2.1194e-05 2.5002 0.04098 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Now we will use binary variables for each region to run the regression with all of them taken into account. So for example for eu the value is 1 for in europe and 0 for not in europe, and the same for all 6 of eu africa apac me northa southa
reg18=lm(happiness~meat+gdppc+eu+af+apac+me+na+sa,meathap)
summary(reg18)##
## Call:
## lm(formula = happiness ~ meat + gdppc + eu + af + apac + me +
## na + sa, data = meathap)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.6949 -0.6219 0.0650 0.6837 2.3080
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.596e+00 4.165e-01 13.434 <2e-16 ***
## meat 3.710e-03 4.479e-03 0.828 0.4091
## gdppc 1.552e-05 6.084e-06 2.551 0.0119 *
## eu -5.439e-01 3.700e-01 -1.470 0.1440
## af -9.974e-01 4.019e-01 -2.482 0.0143 *
## apac -3.892e-01 4.091e-01 -0.951 0.3431
## me -2.794e-01 4.205e-01 -0.664 0.5076
## na -1.033e-01 4.387e-01 -0.235 0.8143
## sa NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.012 on 131 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.2529, Adjusted R-squared: 0.213
## F-statistic: 6.334 on 7 and 131 DF, p-value: 2.045e-06 coeftest(reg18, vcov=vcovHC)##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.5955e+00 3.6636e-01 15.2734 < 2.2e-16 ***
## meat 3.7098e-03 4.4328e-03 0.8369 0.404179
## gdppc 1.5518e-05 7.1752e-06 2.1628 0.032374 *
## eu -5.4393e-01 3.0569e-01 -1.7793 0.077506 .
## af -9.9743e-01 3.4246e-01 -2.9125 0.004216 **
## apac -3.8920e-01 3.4893e-01 -1.1154 0.266722
## me -2.7943e-01 3.9033e-01 -0.7159 0.475334
## na -1.0326e-01 3.3754e-01 -0.3059 0.760161
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Now let’s look at urbanisation
library(ggplot2)
ggplot(meathap, aes(x = urbanpercent, y = happiness)) +
geom_smooth(method = "lm", se = FALSE) +
geom_point()+theme_minimal()+
xlab("Percentage of urbanisation of total population")+ylab("Happiness Score in Cantril Ladder") + ggtitle("Relationship between urbanisation and happiness 2017")## `geom_smooth()` using formula 'y ~ x'
library(ggplot2)
ggplot(meathap, aes(x = urbanpercent, y = meat)) +
geom_smooth(method = "lm", se = FALSE) +
geom_point()+theme_minimal()+
xlab("Percentage of urbanisation of total population")+ylab("Meat Consumption kg/cap") + ggtitle("Relationship between urbanisation and meat consumption 2017")## `geom_smooth()` using formula 'y ~ x'
reg19=lm(happiness~urbanpercent,meathap)
summary(reg19)##
## Call:
## lm(formula = happiness ~ urbanpercent, data = meathap)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.14630 -0.82143 -0.01016 0.79074 2.87957
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.3433 0.2699 16.094 < 2e-16 ***
## urbanpercent 1.8406 0.4189 4.393 2.21e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.069 on 138 degrees of freedom
## Multiple R-squared: 0.1227, Adjusted R-squared: 0.1163
## F-statistic: 19.3 on 1 and 138 DF, p-value: 2.208e-05 coeftest(reg19, vcov=vcovHC)##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.34328 0.29171 14.8892 < 2.2e-16 ***
## urbanpercent 1.84056 0.45641 4.0327 9.082e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 reg20=lm(meat~urbanpercent,meathap)
summary(reg20)##
## Call:
## lm(formula = meat ~ urbanpercent, data = meathap)
##
## Residuals:
## Min 1Q Median 3Q Max
## -51.509 -14.392 -0.736 14.893 54.841
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -10.858 5.420 -2.003 0.0471 *
## urbanpercent 97.634 8.414 11.604 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 21.46 on 138 degrees of freedom
## Multiple R-squared: 0.4939, Adjusted R-squared: 0.4902
## F-statistic: 134.7 on 1 and 138 DF, p-value: < 2.2e-16 coeftest(reg20, vcov=vcovHC)##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -10.8581 4.4307 -2.4507 0.01551 *
## urbanpercent 97.6342 7.9427 12.2923 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#using urbanisation as a control
reg21=lm(happiness~meat + urbanpercent,meathap)
summary(reg21)##
## Call:
## lm(formula = happiness ~ meat + urbanpercent, data = meathap)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.26242 -0.61748 0.08975 0.76643 2.81666
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.466814 0.267493 16.699 < 2e-16 ***
## meat 0.011377 0.004142 2.747 0.00682 **
## urbanpercent 0.729779 0.575391 1.268 0.20684
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.044 on 137 degrees of freedom
## Multiple R-squared: 0.1685, Adjusted R-squared: 0.1564
## F-statistic: 13.88 on 2 and 137 DF, p-value: 3.24e-06 coeftest(reg21, vcov=vcovHC)##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.4668139 0.2870274 15.5623 < 2e-16 ***
## meat 0.0113770 0.0036764 3.0946 0.00239 **
## urbanpercent 0.7297786 0.5654249 1.2907 0.19899
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#both gdppercap and urbanisation as controls
reg22=lm(happiness~meat + gdppc + urbanpercent,meathap)
summary(reg22)##
## Call:
## lm(formula = happiness ~ meat + gdppc + urbanpercent, data = meathap)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.2969 -0.6480 0.1260 0.7732 2.7390
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.632e+00 2.831e-01 16.362 <2e-16 ***
## meat 7.459e-03 4.577e-03 1.630 0.1055
## gdppc 1.163e-05 6.148e-06 1.891 0.0608 .
## urbanpercent 4.961e-01 5.994e-01 0.828 0.4094
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.037 on 135 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.1907, Adjusted R-squared: 0.1727
## F-statistic: 10.6 on 3 and 135 DF, p-value: 2.628e-06 coeftest(reg22, vcov=vcovHC)##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.6321e+00 3.1781e-01 14.5751 < 2e-16 ***
## meat 7.4587e-03 3.9241e-03 1.9008 0.05946 .
## gdppc 1.1626e-05 7.4001e-06 1.5711 0.11851
## urbanpercent 4.9606e-01 6.3019e-01 0.7872 0.43257
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1