Lab 4 - Linear Regression Contd.

Naive Model

1. Outcome = posttest
2. Predictor = treat

Let us run a simple linear regression to estimate the naive model

naive.model <- lm(posttest~treat, data=csiw)
summ(naive.model)

## MODEL INFO:
## Observations: 325 (82 missing obs. deleted)
## Dependent Variable: posttest
## Type: OLS linear regression 
## 
## MODEL FIT:
## F(1,323) = 36.10, p = 0.00
## R² = 0.10
## Adj. R² = 0.10 
## 
## Standard errors: OLS
## -----------------------------------------------
##                     Est.   S.E.   t val.      p
## ----------------- ------ ------ -------- ------
## (Intercept)         0.71   0.35     2.06   0.04
## treat               2.44   0.41     6.01   0.00
## -----------------------------------------------

ANCOVA

In ANCOVA the idea is that we have a continuous variable as the outcome (posttest) and a discrete variable as the predictor of interest (treat) and other covariates to be controlled for (pretest).

In R, ANCOVA is also implemented in the lm() funcion we use for linear regression.

ancova.model <- lm(posttest ~ treat + pretest, data=csiw)
summ(ancova.model)

## MODEL INFO:
## Observations: 272 (135 missing obs. deleted)
## Dependent Variable: posttest
## Type: OLS linear regression 
## 
## MODEL FIT:
## F(2,269) = 24.30, p = 0.00
## R² = 0.15
## Adj. R² = 0.15 
## 
## Standard errors: OLS
## -----------------------------------------------
##                     Est.   S.E.   t val.      p
## ----------------- ------ ------ -------- ------
## (Intercept)         0.22   0.40     0.53   0.59
## treat               2.82   0.45     6.28   0.00
## pretest             0.43   0.14     3.03   0.00
## -----------------------------------------------

Now we can also plot the predicted values as a function of pretest and treat (group membership).

# creating a new column in csiw dataset for the predicted values from ANCOVA model
# csiw$ancova.predicted <- predict(ancova.model)

This error happened because while doing linear regression, R drops the rows where column values are NA. So while doing the analysis, R is dropping 135 observations due to missing data.

Let us remove those rows and redo the analysis.

csiw <- csiw[which(complete.cases(csiw)), ]

Now we will do the ANOVA steps again to see if we get this error

ancova.model <- lm(posttest ~ treat + pretest, data=csiw)
summ(ancova.model)

## MODEL INFO:
## Observations: 258
## Dependent Variable: posttest
## Type: OLS linear regression 
## 
## MODEL FIT:
## F(2,255) = 23.94, p = 0.00
## R² = 0.16
## Adj. R² = 0.15 
## 
## Standard errors: OLS
## -----------------------------------------------
##                     Est.   S.E.   t val.      p
## ----------------- ------ ------ -------- ------
## (Intercept)         0.19   0.42     0.44   0.66
## treat               2.95   0.47     6.34   0.00
## pretest             0.43   0.15     2.90   0.00
## -----------------------------------------------

Dropping the missing values seems to change the treatment effect, this might be an indication that the data is not missing at random. But we will continue with the analysis for now

# creating a new column in csiw dataset for the predicted values from ANCOVA model
csiw$ancova.predicted <- predict(ancova.model)

As you can see now we do not have the error message we had earlier.

Now let us implement the required plot using ggplot2 pakage.

ggplot(csiw, aes(x=pretest, y=ancova.predicted, group=factor(treat),color=factor(treat))) + geom_point()+geom_smooth(method="loess", color="darkred")

Assessing linearity

The assignment asks to plot the residuals from the ANCOVA model against the pretest covariate. We will first calculate the residuals and store it in a new column in csiw and then do the plotting

# storing residuals as new column in csiw
csiw$ancova.residuals <- residuals(ancova.model)

Now we can plot the required graph. Let us add labels to the x and y axis in this plot

ggplot(csiw, aes(x=pretest, y=ancova.residuals, color=factor(treat))) +geom_point()+geom_smooth(method="loess", color="darkred")+
  xlab("Pretest Scores") + ylab("Residuals from ANCOVA model") + labs(title="Residuals vs Pretest - ANCOVA model")

How do we add a quadratic fit? We can create a new variable (pretest.sq) for the pretest variable’s squared value and then replot the graph with the squared variable as the covariate.

csiw$pretest.sq <- csiw$pretest^2

Now let us plot the graph again

ggplot(csiw, aes(x=pretest.sq, y=ancova.residuals, color=factor(treat))) +geom_point()+geom_smooth(method="loess", color="darkred")+
  xlab("Pretest Scores Squared") + ylab("Residuals from ANCOVA model") + labs(title="Residuals vs Pretest squared - ANCOVA model")

We can see that linearity seems to be satisfied, because the residuals are close to zero on average and seem to be distributed almost equally above and below the curve. we have heterosketastity

Quadratic model

Now we have to center the variable pretest and create another variable for the squared value of the centered pretest variable.

# We will use the mutate function to create new variables
pretest.mean <- mean(csiw$pretest)
pretest.mean

## [1] 0.9186047

csiw$pretest.c <- csiw$pretest-pretest.mean
csiw$pretest.csq <- (csiw$pretest.c)^2

Now let us run the ANCOVA again with the centered variable and its square

quadratic.ancova.model <- lm(posttest ~ treat+pretest.c+pretest.csq, data=csiw)
summ(quadratic.ancova.model)

## MODEL INFO:
## Observations: 258
## Dependent Variable: posttest
## Type: OLS linear regression 
## 
## MODEL FIT:
## F(3,254) = 15.98, p = 0.00
## R² = 0.16
## Adj. R² = 0.15 
## 
## Standard errors: OLS
## -----------------------------------------------
##                     Est.   S.E.   t val.      p
## ----------------- ------ ------ -------- ------
## (Intercept)         0.53   0.41     1.29   0.20
## treat               2.94   0.47     6.28   0.00
## pretest.c           0.35   0.23     1.49   0.14
## pretest.csq         0.03   0.07     0.46   0.65
## -----------------------------------------------

Let us look at the predicted values and residual plots of this updated model

csiw$quadratic.ancova.predicted <- predict(quadratic.ancova.model)
csiw$quadratic.ancova.residuals <- residuals(quadratic.ancova.model)

Now the plots

par(mfrow=c(1,2))

ggplot(csiw, aes(x=pretest.c, y = quadratic.ancova.predicted, color=factor(treat))) + geom_point()+
  geom_smooth(method="loess", color="darkred")+  xlab("Pretest Scores Centered") + ylab("Predicted values from Quadratic ANCOVA model") + labs(title="Predicted values vs Pretest squared - Quadratic ANCOVA model")

ggplot(csiw, aes(x=pretest.c, y = quadratic.ancova.residuals, color=factor(treat))) + geom_point()+
  geom_smooth(method="loess", color="darkred")+  xlab("Pretest Scores Centered") + ylab("Residuals from Quadratic ANCOVA model") + labs(title="Residuals vs Pretest squared - Quadratic ANCOVA model")

Searching for confounders

In order to identify confounders we need three kinds of tests:

1. ANOVA - aov() - for continuous outcome on discrete predictor
   
2. Linear Regression - lm () - for continuous outcome on continuous predictor
3. Chi Square test - chisq.test() - for discrete outcome on discrete predictor

The list of confounders to check:

1. Dummies from achievement level
2. Dummy for grade level

Let us see some examples. Let us see

# Dummies from achievement level eg. achievement level = high
## corresponding dummy in the dataset = group_high

### outcome vs. predictor = posttest vs. group_high = continuous vs. discrete = ANOVA

summary(aov(posttest~group_high, data=csiw))

##              Df Sum Sq Mean Sq F value   Pr(>F)    
## group_high    1  474.4   474.4    43.2 2.75e-10 ***
## Residuals   256 2811.2    11.0                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(aov(posttest~group_average, data=csiw))

##                Df Sum Sq Mean Sq F value Pr(>F)
## group_average   1      2   2.322   0.181  0.671
## Residuals     256   3283  12.825

summary(aov(posttest~group_low, data=csiw))

##              Df Sum Sq Mean Sq F value Pr(>F)
## group_low     1     19   19.05   1.493  0.223
## Residuals   256   3267   12.76

summary(aov(posttest~group_learndis, data=csiw))

##                 Df Sum Sq Mean Sq F value   Pr(>F)    
## group_learndis   1  416.5   416.5   37.16 3.99e-09 ***
## Residuals      256 2869.2    11.2                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

### treatment vs. predictor = treat vs group_learndis = Chi Square test
csiw$treat.factor <- factor(csiw$treat)
chisq.test(csiw$treat.factor,csiw$group_high)

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  csiw$treat.factor and csiw$group_high
## X-squared = 9.5627, df = 1, p-value = 0.001986

chisq.test(csiw$treat.factor,csiw$group_learndis)

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  csiw$treat.factor and csiw$group_learndis
## X-squared = 7.7297, df = 1, p-value = 0.005432

From the results above we can see that this is a potential confounder. Now let us add this to the model and see the output

confounder.check.model <- lm(posttest ~ treat+pretest.c+pretest.csq+group_high + group_learndis, data=csiw)
summ(confounder.check.model)

## MODEL INFO:
## Observations: 258
## Dependent Variable: posttest
## Type: OLS linear regression 
## 
## MODEL FIT:
## F(5,252) = 20.67, p = 0.00
## R² = 0.29
## Adj. R² = 0.28 
## 
## Standard errors: OLS
## ---------------------------------------------------
##                         Est.   S.E.   t val.      p
## -------------------- ------- ------ -------- ------
## (Intercept)             0.82   0.43     1.92   0.06
## treat                   2.17   0.45     4.86   0.00
## pretest.c               0.07   0.22     0.32   0.75
## pretest.csq             0.08   0.06     1.35   0.18
## group_high              2.12   0.46     4.62   0.00
## group_learndis         -1.73   0.50    -3.48   0.00
## ---------------------------------------------------

Let us look at the predicted values and residual plots of this updated model

csiw$confounder.check.model <- predict(confounder.check.model)
csiw$confounder.check.residuals <- residuals(confounder.check.model)

Now the plots

par(mfrow=c(1,2))

ggplot(csiw, aes(x=pretest.c, y = confounder.check.model, color=factor(treat))) + geom_point()+
  geom_smooth(method="loess", color="darkred")+  xlab("Pretest Scores Centered") + ylab("Predicted values from Quadratic ANCOVA model") + labs(title="Predicted values vs Pretest squared - Confounder Check model")

ggplot(csiw, aes(x=pretest.c, y = confounder.check.residuals, color=factor(treat))) + geom_point()+
  geom_smooth(method="loess", color="darkred")+  xlab("Pretest Scores Centered") + ylab("Residuals from Quadratic ANCOVA model") + labs(title="Residuals vs Pretest squared - Confounder Check model")

Now let us look at the results we have gotten so far. For this we will use the export_summs() function from the jtools package

export_summs(naive.model, ancova.model, quadratic.ancova.model, confounder.check.model)

Since the remaining list of covariates to check for confounding are all discrete, we can repeat the steps from above with different covariates.

Exercise: Do the confounder check for grade dummy called grade_2.  
          grade_2 = 1, if Grade = 5 and grade_2 = 0 if Grade = 4.

Heterogeneity

interaction.model <- lm(posttest ~ treat*grade_1 +pretest.c+pretest.csq+group_high + group_learndis, data=csiw)
summ(interaction.model)

## MODEL INFO:
## Observations: 258
## Dependent Variable: posttest
## Type: OLS linear regression 
## 
## MODEL FIT:
## F(7,250) = 14.75, p = 0.00
## R² = 0.29
## Adj. R² = 0.27 
## 
## Standard errors: OLS
## ---------------------------------------------------
##                         Est.   S.E.   t val.      p
## -------------------- ------- ------ -------- ------
## (Intercept)             1.14   0.64     1.79   0.07
## treat                   1.80   0.69     2.63   0.01
## grade_1                -0.52   0.77    -0.67   0.50
## pretest.c               0.09   0.23     0.42   0.68
## pretest.csq             0.08   0.06     1.27   0.21
## group_high              2.13   0.46     4.62   0.00
## group_learndis         -1.77   0.50    -3.51   0.00
## treat:grade_1           0.63   0.89     0.70   0.48
## ---------------------------------------------------

interaction.model2 <- lm(posttest ~ treat*grade_2 +pretest.c+pretest.csq+group_high + group_learndis, data=csiw)
summ(interaction.model2)

## MODEL INFO:
## Observations: 258
## Dependent Variable: posttest
## Type: OLS linear regression 
## 
## MODEL FIT:
## F(7,250) = 14.75, p = 0.00
## R² = 0.29
## Adj. R² = 0.27 
## 
## Standard errors: OLS
## ---------------------------------------------------
##                         Est.   S.E.   t val.      p
## -------------------- ------- ------ -------- ------
## (Intercept)             0.63   0.52     1.20   0.23
## treat                   2.43   0.58     4.16   0.00
## grade_2                 0.52   0.77     0.67   0.50
## pretest.c               0.09   0.23     0.42   0.68
## pretest.csq             0.08   0.06     1.27   0.21
## group_high              2.13   0.46     4.62   0.00
## group_learndis         -1.77   0.50    -3.51   0.00
## treat:grade_2          -0.63   0.89    -0.70   0.48
## ---------------------------------------------------

Exercise: Do the interaction check for treatment and prior achievement level on 4 groups.

Checking assumptions about the random part of model

plot(interaction.model2)