--- title: "Regression with a Single Categorical Explanatory Variable" subtitle: Analysis template author: Sunny Avry date: "January 27th, 2020" output: html_document: toc: false theme: united --- This R notebook can be used as a template to compute regression with a single categorical explanatory variable. It is usable in Rstudio. The content is adapated from the mooc Multilevel modelling online course from the Centre for Mutlilevel Modelling at the University of Bristol (http://www.bristol.ac.uk/cmm/learning/online-course/). # Packages {r message=FALSE, warning=FALSE} library(xlsx) #read xlsx files  # Data importation {r} #Text file mydata <- read.table(file = "3.2.txt", sep = ",", header = TRUE) #put the file in the same directory as this notebook #Excel file #mydata <- read.xlsx("myfile.xlsx", sheetName = "mysheetname")  # Data exploration {r} dim(mydata) #number of rows and columns  {r} str(mydata) #basic description of the dataframe  {r} mydata[1:20, ] #first 20 rows of the data  # Dummy categorical variable Having viewed the data we will now examine score (ratio variable : point score calculated from awards in Standard grades. Scores range from 0 to 75, with a higher score indicating a higher attainment ) and female (categorical variable : girls 1; boys 0 ## Dependent variable (score) ### Histogram of the data {r} #continuous variable hist(mydata$score, xlim = c(0,80))  ### Summary of the data {r} summary(mydata$score)  ### Standard deviation of the data {r} sd(mydata$score)  ## Independent (explanatory) variable (female) ### Frequencies, percentages and cumulated percentages {r} mytable <- table(mydata$female) mytablecomb <- cbind(mytable, prop.table(mytable), cumsum(prop.table(mytable))) colnames(mytablecomb) <- c("Freq", "Perc", "Cum") mytablecomb  ## Relationship between explanatory and dependant variables ### Frequence, mean and sd of the dependent variable for each value of the explanatory variable {r} l <- tapply(mydata$score, factor(mydata$female), length) m <- tapply(mydata$score, factor(mydata$female), mean) s <- tapply(mydata$score, factor(mydata$female), sd) tableScore <- cbind("Freq" = l, "mean(score)" = m, "sd(score)" = s) tableScore <- rbind(tableScore, c(sum(tableScore[, 1]), mean(mydata$score), sd(mydata$score))) rownames(tableScore)[dim(tableScore)[1]] <- "Total" tableScore  ## Linear regression $$dependent_i = \beta_0 + \beta_1 * explanatory_i + e_i$$ {r} fit <- lm(score ~ female, data = mydata) summary(fit)  $$\hat{score}_i = 29.83 + 2.41 * female_i$$ ### Residual (unexplained) variance {r} model <- summary(fit) model$sigma^2  ### % Explained variance {r} model$r.squared  ## Hypothesis testing The test statistic is calculated by dividing the estimated slope by its standard error {r} model$coefficients[2,1] / model$coefficients[2,2] #significant if larger than 1.96 or -1.96, the critical value at the 5% level of significance  # Categorical variable with more than 2 categories We will now examine score (ratio variable : point score calculated from awards in Standard grades. Scores range from 0 to 75, with a higher score indicating a higher attainment) and parental social class (categorical variable : managerial and professional 1; intermediate 2; working 3; unclassified 4) ## Dependent variable (score) ### Histogram of the data {r} #continuous variable hist(mydata$score, xlim = c(0,80))  ### Summary of the data {r} summary(mydata$score)  ### Standard deviation of the data {r} sd(mydata$score)  ## Independent (explanatory) variable (parental social class) ### Frequencies, percentages and cumulated percentages {r} mytable <- table(mydata$sclass) mytablecomb <- cbind(mytable, prop.table(mytable), cumsum(prop.table(mytable))) colnames(mytablecomb) <- c("Freq", "Perc", "Cum") mytablecomb  ## Relationship between explanatory and dependant variables ### Frequence, mean and sd of the dependent variable for each value of the explanatory variable {r} l <- tapply(mydata$score, factor(mydata$sclass), length) m <- tapply(mydata$score, factor(mydata$sclass), mean) s <- tapply(mydata$score, factor(mydata$sclass), sd) tableScore <- cbind("Freq" = l, "mean(score)" = m, "sd(score)" = s) tableScore <- rbind(tableScore, c(sum(tableScore[, 1]), mean(mydata$score), sd(mydata$score))) rownames(tableScore)[dim(tableScore)[1]] <- "Total" tableScore  ### Creation of dummy variables {r} mydata$sclass1 <- mydata$sclass == 1 mydata$sclass2 <- mydata$sclass == 2 mydata$sclass3 <- mydata$sclass == 3 mydata$sclass4 <- mydata$sclass == 4 mydata[1:20, ]  ## Linear regression ### A first model with sclass3 as reference $$score_i = \beta_0 + \beta_1 * sclass1_i + \beta_2 * sclass2_i + \beta_3 * sclass4_i + e_i$$ {r} fit <- lm(score ~ sclass1 + sclass2 + sclass4, data = mydata) summary(fit)  $$\hat{score}_i = 24.1 + 15.14 * sclass1_i + 7.77 * sclass2_i - 2.75 * sclass4_i$$ ### Residual (unexplained) variance {r} model <- summary(fit) model$sigma^2  ### % Explained variance {r} model$r.squared  ## Hypothesis testing The test statistic is calculated by dividing the estimated slope by its standard error {r} #sclass3 vs.sclass1 model$coefficients[2,1] / model$coefficients[2,2] #sclass3 vs.sclass2 model$coefficients[3,1] / model$coefficients[3,2] #sclass3 vs.sclass4 model$coefficients[4,1] / model$coefficients[4,2] #significant if larger than 1.96 or -1.96, the critical value at the 5% level of significance  ### A second model with sclass2 as reference $$score_i = \beta_0 + \beta_1 * sclass1_i + \beta_2 * sclass3_i + \beta_3 * sclass4_i + e_i$$ {r} fit <- lm(score ~ sclass1 + sclass3 + sclass4, data = mydata) summary(fit)  $$\hat{score}_i = 31.87 + 7.37 * sclass1_i - 7.77 * sclass3_i - 10.52 * sclass4_i$$ ### Residual (unexplained) variance {r} model <- summary(fit) model$sigma^2  ### % Explained variance {r} model$r.squared  ## Hypothesis testing The test statistic is calculated by dividing the estimated slope by its standard error {r} #sclass2 vs.sclass1 model$coefficients[2,1] / model$coefficients[2,2] #sclass2 vs.sclass3 model$coefficients[3,1] / model$coefficients[3,2] #sclass2 vs.sclass4 model$coefficients[4,1] / model$coefficients[4,2] #significant if larger than 1.96 or -1.96, the critical value at the 5% level of significance  ## Non-linear regression We will now examine score (ratio variable : point score calculated from awards in Standard grades. Scores range from 0 to 75, with a higher score indicating a higher attainment ) and cohort90 taken as a categorical variable (cohorts: 1984, 1986, 1988, 1990, 1996 and 1998). The cohort90 variable is calculated by subtracting 1990 from each value. Thus values range from -6 (corresponding to 1984) to 8 (1998), with 1990 coded as zero. ### Frequencies {r} table(mydata$cohort90)  ### Creation of dummy variables {r} mydata$cohort90yr84 <- mydata$cohort90 == -6 mydata$cohort90yr86 <- mydata$cohort90 == -4 mydata$cohort90yr88 <- mydata$cohort90 == -2 mydata$cohort90yr90 <- mydata$cohort90 == 0 mydata$cohort90yr96 <- mydata$cohort90 == 6 mydata$cohort90yr98 <- mydata\$cohort90 == 8 mydata[1:20, ]  ### A model with cohort90yr88 as a reference $$score_i = \beta_0 + \beta_1 * cohort90yr84_i + \beta_2 * cohort90yr86_i + \\ \beta_3 * cohort90yr90_i + \beta_4 * cohort90yr96_i + \beta_5 * cohort90yr98_i + e_i$$ {r} fit <- lm(score ~ cohort90yr84 + cohort90yr86 + cohort90yr90 + cohort90yr96 + cohort90yr98, data = mydata) summary(fit)  $$\hat{score}_i = 28.52 + -4.87 * cohort90yr84_i - 3.75 * cohort90yr86_i + \\ 0.58 * cohort90yr86_i + 10.91 * cohort90yr96_i + 12.81 * cohort90yr98_i$$ The coefficients increase with year, as assumed in a linear model. We can treat cohort90 as continuous in all subsequent analyses. # Sources - Multilevel modelling online course from the Centre for Mutlilevel Modelling of the University of Bristol (http://www.bristol.ac.uk/cmm/learning/online-course/).