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

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

library(ggplot2)
ggplot(oind, aes(x = Hours, y = Income, color=Gender)) +
    geom_point() +
    geom_function(fun=function(x) mod1$coefficients["(Intercept)"] +
                                  mod1$coefficients["Hours"]*x, size=1, color="red") +
    geom_function(fun=function(x) mod1$coefficients["(Intercept)"] +
                                  mod1$coefficients["Genderfemale"] + 
                                  mod1$coefficients["Hours"]*x, size=1, color="cyan")
    ggsave("predscatR2.png", units="cm", width=20, height=11)

nlsw88 <- read.dta("https://teaching.sociology.ul.ie/so4046/nlsw88.dta")
summary(mod2 <- 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) mod2$coefficients["(Intercept)"] +
                                  mod2$coefficients["ttl_exp"]*x + 
                                  mod2$coefficients["grade"]*5, size=0.4, aes(col="Grade =  5")) +
    geom_function(fun=function(x) mod2$coefficients["(Intercept)"] +
                                  mod2$coefficients["ttl_exp"]*x + 
                                  mod2$coefficients["grade"]*10, size=0.4, aes(col="Grade = 10")) +
    geom_function(fun=function(x) mod2$coefficients["(Intercept)"] +
                                  mod2$coefficients["ttl_exp"]*x + 
                                  mod2$coefficients["grade"]*15, size=0.4, aes(col="Grade = 15")) +
    geom_function(fun=function(x) mod2$coefficients["(Intercept)"] +
                                  mod2$coefficients["ttl_exp"]*x + 
                                  mod2$coefficients["grade"]*20, size=0.4, aes(col="Grade = 20")) +
    labs(col="Predictions")
    ggsave("continterR1.png", units="cm", width=20, height=11)

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

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

modbase <- lm(data=oind, Income ~ Hours)
modplus <- lm(data=oind, Income ~ Hours + JobScore + Gender)
anova(modbase, modplus)

summary(lm(formula = wage ~ ttl_exp + grade + collgrad, data = nlsw88))

lm(formula = wage ~ ttl_exp + occupation, data = nlsw88)