## (C) (cc by-sa) Wouter van Atteveldt, file generated juni 23 2015Note on the data used in this howto: This data can be downloaded from http://piketty.pse.ens.fr/files/capital21c/en/xls/, but the excel format is a bit difficult to parse at it is meant to be human readable, with multiple header rows etc. For that reason, I’ve extracted csv files for some interesting tables that I’ve uploaded to http://vanatteveldt.com/uploads/rcourse/data
Basic Modeling
In this hands-on we continue with the capital variable created in the transforming data howto. You can also download this variable from the course pages:
load("data/capital.rdata")
head(capital)## Year Country Public Private Total
## 1 1970 Australia 0.61 3.30 3.91
## 2 1970 Canada 0.37 2.47 2.84
## 3 1970 France 0.41 3.10 3.51
## 4 1970 Germany 0.88 2.25 3.13
## 5 1970 Italy 0.20 2.39 2.59
## 6 1970 Japan 0.61 2.99 3.60T-tests
First, let’s split our countries into two groups, anglo-saxon countries and european countries (plus Japan): We can use the ifelse command here combined with the %in% operator
anglo = c("U.S.", "U.K.", "Canada", "Australia")
capital$Group = ifelse(capital$Country %in% anglo, "anglo", "european")
table(capital$Group)##
## anglo european
## 164 205Now, let’s see whether capital accumulation is different between these two groups.
t.test(capital$Private ~ capital$Group)##
## Welch Two Sample t-test
##
## data: capital$Private by capital$Group
## t = -4.6664, df = 289.34, p-value = 4.692e-06
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.7775339 -0.3162154
## sample estimates:
## mean in group anglo mean in group european
## 3.748232 4.295106So, according to this test capital accumulation is indeed significantly higher in European countries than in Anglo-Saxon countries. Of course, the data here are not independently distributed since the data in the same year in different countries is related (as are data in subsequent years in the same country, but let’s ignore that for the moment) We could also do a paired t-test of average accumulation per year per group by first using the cast command to aggregate the data. Note that we first remove the NA values (for Spain).
library(reshape2)
capital = na.omit(capital)
pergroup = dcast(capital, Year ~ Group, value.var="Private", fun.aggregate=mean)
head(pergroup)## Year anglo european
## 1 1970 3.0625 2.6825
## 2 1971 3.1475 2.7425
## 3 1972 3.2450 2.9000
## 4 1973 3.1800 2.9500
## 5 1974 3.1125 3.0025
## 6 1975 3.0300 3.1325Let’s plot the data to have a look at the lines:
plot(pergroup$Year, pergroup$european, type="l", xlab="Year", ylab="Capital accumulation")
lines(pergroup$Year, pergroup$anglo, lty=2)
legend("top", lty=c(1,2), legend=c("European", "Anglo-Saxon"))
So initially capital is higher in the Anglo-Saxon countries, but the European countries overtake quickly and stay higher.
Now, a paired-sample t-test again shows a significant difference between the two:
t.test(pergroup$anglo, pergroup$european, paired=T)##
## Paired t-test
##
## data: pergroup$anglo and pergroup$european
## t = -6.5332, df = 40, p-value = 8.424e-08
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.6007073 -0.3168537
## sample estimates:
## mean of the differences
## -0.4587805Anova
We can also use a one-way Anova to see whether accumulation differs per country. Let’s first do a box-plot to see how different the countries are. Plot by default gives a box plot of a formula with a nominal independeny variable
plot(capital$Private ~ capital$Country)
So, it seems that in fact a lot of countries are quite similar, with some extreme cases of high capital accumulation. (also, it seems that including Japan in the Europeaqn countries might have been a mistake). We use the aov function for this, the anova function is meant to analyze already fitted models, as will be shown below.
m = aov(capital$Private ~ capital$Country)
summary(m)## Df Sum Sq Mean Sq F value Pr(>F)
## capital$Country 8 201.3 25.158 30.78 <2e-16 ***
## Residuals 343 280.3 0.817
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1So in fact there is a significant difference. We can use pairwise.t.test to perform
posthoc = pairwise.t.test(capital$Private, capital$Country, p.adj = "bonf")
round(posthoc$p.value, 2)## U.S. Japan Germany France U.K. Italy Canada Australia
## Japan 0.00 NA NA NA NA NA NA NA
## Germany 0.00 0.00 NA NA NA NA NA NA
## France 1.00 0.00 0.06 NA NA NA NA NA
## U.K. 1.00 0.00 0.00 1.00 NA NA NA NA
## Italy 0.02 0.01 0.00 0.00 0.31 NA NA NA
## Canada 0.02 0.00 1.00 0.27 0.00 0.00 NA NA
## Australia 1.00 0.00 0.00 1.00 1.00 0.46 0 NA
## Spain 0.00 1.00 0.00 0.00 0.00 0.01 0 0Linear models
A more generic way of fitting models is using the lm command. In fact, aov is a wrapper around lm. Let’s model private capital as a function of country and public capital:
m = lm(Private ~ Country + Public, data=capital)
summary(m)##
## Call:
## lm(formula = Private ~ Country + Public, data = capital)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.4457 -0.4091 -0.1076 0.2601 2.8346
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.7241 0.1468 32.184 < 2e-16 ***
## CountryJapan 1.8183 0.1753 10.374 < 2e-16 ***
## CountryGermany -0.7776 0.1721 -4.518 8.63e-06 ***
## CountryFrance -0.3337 0.1729 -1.930 0.054420 .
## CountryU.K. 0.3910 0.1733 2.257 0.024643 *
## CountryItaly -0.7911 0.2198 -3.600 0.000366 ***
## CountryCanada -1.6928 0.1949 -8.685 < 2e-16 ***
## CountryAustralia 0.6421 0.1768 3.633 0.000323 ***
## CountrySpain 0.8331 0.2120 3.930 0.000103 ***
## Public -1.8144 0.1660 -10.933 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7793 on 342 degrees of freedom
## Multiple R-squared: 0.5687, Adjusted R-squared: 0.5573
## F-statistic: 50.1 on 9 and 342 DF, p-value: < 2.2e-16As you can see, R automatically creates dummy values for nominal values, using the first value (U.S. in this case) as reference category. An alternative is to remove the intercept and create a dummy for each country:
m = lm(Private ~ Country + Public - 1, data=capital)
summary(m)##
## Call:
## lm(formula = Private ~ Country + Public - 1, data = capital)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.4457 -0.4091 -0.1076 0.2601 2.8346
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## CountryU.S. 4.7241 0.1468 32.18 <2e-16 ***
## CountryJapan 6.5424 0.1676 39.05 <2e-16 ***
## CountryGermany 3.9465 0.1474 26.77 <2e-16 ***
## CountryFrance 4.3904 0.1383 31.74 <2e-16 ***
## CountryU.K. 5.1151 0.1587 32.23 <2e-16 ***
## CountryItaly 3.9330 0.1334 29.48 <2e-16 ***
## CountryCanada 3.0313 0.1221 24.83 <2e-16 ***
## CountryAustralia 5.3662 0.1725 31.11 <2e-16 ***
## CountrySpain 5.5572 0.1596 34.82 <2e-16 ***
## Public -1.8144 0.1660 -10.93 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7793 on 342 degrees of freedom
## Multiple R-squared: 0.9666, Adjusted R-squared: 0.9657
## F-statistic: 991.2 on 10 and 342 DF, p-value: < 2.2e-16(- 1 removes the intercept because there is an implicit +1 constant for the intercept in the regression formula)
You can also introduce interaction terms by using either the : operator (which only creates the interaction term) or the * (which creates a full model including the main effects). To keep the model somewhat parsimonious, let’s use the country group rather than the country itself
m1 = lm(Private ~ Group + Public, data=capital)
m2 = lm(Private ~ Group + Public + Group:Public, data=capital)A nice package to display multiple regression results side by side is the screenreg function from the texreg package:
library(texreg)## Version: 1.35
## Date: 2015-04-25
## Author: Philip Leifeld (University of Konstanz)
##
## Please cite the JSS article in your publications -- see citation("texreg").screenreg(list(m1, m2))##
## ============================================
## Model 1 Model 2
## --------------------------------------------
## (Intercept) 3.97 *** 3.75 ***
## (0.11) (0.13)
## Groupeuropean 0.47 *** 0.78 ***
## (0.12) (0.16)
## Public -0.49 *** -0.01
## (0.14) (0.22)
## Groupeuropean:Public -0.83 **
## (0.28)
## --------------------------------------------
## R^2 0.09 0.11
## Adj. R^2 0.08 0.10
## Num. obs. 352 352
## RMSE 1.12 1.11
## ============================================
## *** p < 0.001, ** p < 0.01, * p < 0.05So, there is a significant interaction effect which displaces the main effect of public wealth.
Comparing and diagnosing models
A relevant question can be whether a model with an interaction effect is in fact a better model than the model without the interaction. This can be investigated with an anova of the model fits of the two models:
m1 = lm(Private ~ Group + Public, data=capital)
m2 = lm(Private ~ Group + Public + Group:Public, data=capital)
anova(m1, m2)## Analysis of Variance Table
##
## Model 1: Private ~ Group + Public
## Model 2: Private ~ Group + Public + Group:Public
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 349 440.02
## 2 348 429.36 1 10.661 8.641 0.003506 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1So, the interaction term is in fact a significant improvement of the model. Apparently, in European countries private capital is accumulated faster in those times that the government goes into depth.
After doing a linear model it is a good idea to do some diagnostics. We can ask R for a set of standard plots by simply calling plot on the model fit. We use the parameter (par) mfrow here to put the four plots this produces side by side.
old.settings = par(mfrow=c(2,2))
plot(m)
par(old.settings)See http://www.statmethods.net/stats/rdiagnostics.html for a more exhausitve list of model diagnostics.