2015-05-25.r

d = read.csv("data/income_topdecile.csv")

head(d)
d = na.omit(d)
d


t.test(d$U.S., d$U.K., paired=T)
cor.test(d$U.S., d$U.K., paired=T)
chisq.test(d$U.S., d$U.K.)

library(reshape2)
x = melt(d, id.vars="Year")
head(x)

x = x[x$variable %in% c("U.S.", "U.K."), ]

x = x[x$variable == "U.S." & xvariable ==  "U.K.", ]
x = x[(x$variable == "U.S." | x$variable ==  "U.K.") & x$Year > 1945, ]
x = x[x$variable == "U.S." | (x$variable ==  "U.K." & x$Year > 1945), ]

table(x$variable)

t.test(x$value ~ x$variable)
t.test(value ~ variable, data=x)


x = melt(d, id.vars="Year")
head(x)

# dummy plus year
lm(value ~ variable + Year, data=x)

# dummy plus year without intercept
lm(value ~ variable + Year - 1, data=x)

# full interaction model: interaction each dummy:year except base category
lm(value ~ variable*Year, data=x)

# only interaction, no main effects: interaction each dummy:year includingt US
lm(value ~ variable:Year, data=x)

# dummy plus interaction, no main effects: interaction each dummy:year includingt US (no main dummy for US)
lm(value ~ variable + variable:Year, data=x)

# dummy plus interaction, no main effects: interaction each dummy:year includingt US (no intercept, so also main dummy US)
lm(value ~ variable + variable:Year - 1, data=x)


# dummy only for germany?
x$germany = x$variable == "Germany"
lm(value ~ Year + germany, data=x)


# calculation in model spec
lm(value ~ Year**2, data=x)
lm(value ~ log(Year), data=x)


aov(value ~ variable, data=x)


# how to get model parameters?
m = lm(value ~ Year, data=x)

# summary gives b, r etc
summary(m)

# but not beta. You want beta?
library(quantpsych)
lm.beta(m)


fitted(m)
resid(m)

plot(m)
par(mfrow=c(2,2))
plot(m)

dev.off()
plot(x=x$Year, y=x$value)

x$resid = resid(m)