library(foreign)
oind <- read.dta("https://teaching.sociology.ul.ie/so5032/oind.dta")

summary(mod1 <- lm(data=oind, Income ~ Hours + Gender))

summary(mod1m <- lm(data=subset(oind, Gender=="male"), Income ~ Hours))

summary(mod1f <- lm(data=subset(oind, Gender=="female"), Income ~ Hours))

library(ggplot2)
ggplot(subset(oind,Income<=5000), aes(x = Hours, y = Income, color=Gender)) +
    geom_point() +
    geom_function(fun=function(x) mod1f$coefficients["(Intercept)"] +
                                  mod1f$coefficients["Hours"]*x, size=1, color="darkgreen") +
    geom_function(fun=function(x) mod1m$coefficients["(Intercept)"] +
                                  mod1m$coefficients["Hours"]*x, size=1, color="darkred")
    ggsave("predscatR3a.png", units="cm", width=20, height=11)

oind$intvar <- oind$Hours * (oind$Gender=="female")

summary(lm(data=oind, Income ~ Hours + Gender + intvar))

summary(lm(data=oind, Income ~ Hours*Gender))

lm(data=oind, Income ~ Hours*Gender)

nlsw88 <- read.dta("https://teaching.sociology.ul.ie/so4046/nlsw88.dta")
summary(modnls1 <- lm(data=nlsw88, wage ~ ttl_exp + grade))

library(ggplot2)
ggplot(nlsw88, aes(x = ttl_exp, y = wage)) +
    geom_point(color="red", alpha=0.1) +
    geom_function(fun=function(x) modnls1$coefficients["(Intercept)"] +
                                  modnls1$coefficients["ttl_exp"]*x + 
                                  modnls1$coefficients["grade"]*5, size=0.4, aes(col="Grade =  5")) +
    geom_function(fun=function(x) modnls1$coefficients["(Intercept)"] +
                                  modnls1$coefficients["ttl_exp"]*x + 
                                  modnls1$coefficients["grade"]*10, size=0.4, aes(col="Grade = 10")) +
    geom_function(fun=function(x) modnls1$coefficients["(Intercept)"] +
                                  modnls1$coefficients["ttl_exp"]*x + 
                                  modnls1$coefficients["grade"]*15, size=0.4, aes(col="Grade = 15")) +
    geom_function(fun=function(x) modnls1$coefficients["(Intercept)"] +
                                  modnls1$coefficients["ttl_exp"]*x + 
                                  modnls1$coefficients["grade"]*20, size=0.4, aes(col="Grade = 20")) +
    labs(col="Predictions")
    ggsave("predscat5flatR.png", units="cm", width=20, height=11)

summary(modnls2 <- lm(data=nlsw88, wage ~ ttl_exp*grade))

library(ggplot2)
ggplot(nlsw88, aes(x = ttl_exp, y = wage)) +
    geom_point(color="red", alpha=0.1) +
    geom_function(fun=function(x) modnls2$coefficients["(Intercept)"] +
                                  modnls2$coefficients["ttl_exp"]*x + 
                                  modnls2$coefficients["grade"]*5 + 
                                  modnls2$coefficients["ttl_exp:grade"]*5*x, size=0.4, aes(col="Grade =  5")) +
    geom_function(fun=function(x) modnls2$coefficients["(Intercept)"] +
                                  modnls2$coefficients["ttl_exp"]*x + 
                                  modnls2$coefficients["grade"]*10 + 
                                  modnls2$coefficients["ttl_exp:grade"]*10*x, size=0.4, aes(col="Grade = 10")) +
    geom_function(fun=function(x) modnls2$coefficients["(Intercept)"] +
                                  modnls2$coefficients["ttl_exp"]*x + 
                                  modnls2$coefficients["grade"]*15 + 
                                  modnls2$coefficients["ttl_exp:grade"]*15*x, size=0.4, aes(col="Grade = 15")) +
    geom_function(fun=function(x) modnls2$coefficients["(Intercept)"] +
                                  modnls2$coefficients["ttl_exp"]*x + 
                                  modnls2$coefficients["grade"]*20 + 
                                  modnls2$coefficients["ttl_exp:grade"]*20*x, size=0.4, aes(col="Grade = 20")) +
    labs(col="Predictions")
    ggsave("predscat5intR.png", units="cm", width=20, height=11)

summary(lm(data=oind, Income ~ Hours*Gender))