R Markdown

This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.

When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:

Step. 1 Read Data

library(haven)
ANCOVA.raw <- read_sav("PISA_tawian2022_trimmed_lab.sav")
View(ANCOVA.raw)

ANCOVA<-ANCOVA.raw[c(13,20,22)]
head(ANCOVA)  #the first six rows
str(ANCOVA)   #check the scale

Step. 2 Variable Change

ANCOVA$Workpay_new<-as.factor(ANCOVA$Workpay_new)
ANCOVA$BELONG<-as.numeric(ANCOVA$BELONG)

Step. 3 Checking the homogeneity of regression coefficient (slope) assumption

TEST1    <- aov(PV1READ ~ BELONG*Workpay_new ,data=ANCOVA)
summary(TEST1) #looking at the effect of interaction term BELONG*Workpay_new

Step. 4 Graphs for regression lines per groups

library(ggplot2)
ggplot( ANCOVA, aes(BELONG, PV1READ)) +
  facet_grid(. ~ Workpay_new) +
  geom_point() +
  stat_smooth(method="lm")

ggplot( ANCOVA, aes(x = BELONG, y = PV1READ, color = Workpay_new))+
  geom_point() +
  stat_smooth(method="lm")

Step. 5 Testing the homogeneity of variance assumption

bartlett.test(PV1READ ~ Workpay_new, data=ANCOVA)

library(car)
leveneTest(PV1READ ~ Workpay_new, data=ANCOVA)

Step. 6 Remove the interaction term and run the ANCOVA model (since homogeneity of regression coefficient is met)

ANCOVA.OP <- aov(PV1READ ~ BELONG + Workpay_new, data=ANCOVA)
summary(ANCOVA.OP)

Step. 7 Alternatives ways of conducting ANCOVA, need to specify type I SS

library(rstatix)
anova_test(data = ANCOVA, formula =PV1READ ~ BELONG + Workpay_new, type = 1, detailed = TRUE)

Step. 8 Original outcome mean value for each group

aggregate(PV1READ ~ Workpay_new,mean ,data=ANCOVA)

Step. 9 Adjusted outcome mean value for each group

library(effects)
adj.mean<-effect("Workpay_new", ANCOVA.OP)
data.frame(adj.mean)

Step. 10 Alternative:Adjusted outcome mean value for each group

adj_means.1 <- emmeans_test(data = ANCOVA, formula = PV1READ ~ Workpay_new, covariate = BELONG)
get_emmeans(adj_means.1)

mean(ANCOVA$BELONG)

Step. 11 Plot the adjusted mean

plot(adj.mean)

Step. 12 Post hoc comparison using Tukey method since the main effect for group is sig

library(multcomp)
posthoc <- glht(ANCOVA.OP, linfct = mcp(Workpay_new = "Tukey"))
summary(posthoc)

Step. 13 Post hoc comparison: using bonferroni method, adjusted p-value

emmeans_test(data = ANCOVA, formula = PV1READ ~ Workpay_new, covariate = BELONG, p.adjust.method = "bonferroni")

Step. 14 Effect size

library(lsr)
etaSquared(ANCOVA.OP)

Step. 15 Comparison with ANOVA model

ANOVA<- aov(PV1READ ~ Workpay_new, data=ANCOVA)
summary(ANOVA) 
summary(ANCOVA.OP)