DOI

Load the library

library(conflicted)
library(rstatix)
library(gridExtra)
library(dplyr)
library(ggpubr)
library(rstatix)
library(broom)
library(readr)
library(tidyverse)

Load the dataset

df58=read_csv("data58.csv")
df58
## # A tibble: 58 × 3
##    VLEs  pretest posttest
##    <chr>   <dbl>    <dbl>
##  1 VLE1     43.3     48.3
##  2 VLE1     48.3     90.6
##  3 VLE1     51.7     55  
##  4 VLE1     50       43.3
##  5 VLE1     38.3     51.7
##  6 VLE1     40       50.6
##  7 VLE1     63.9     65  
##  8 VLE1     56.7     76.7
##  9 VLE1     35       41.7
## 10 VLE1     40       70  
## # ℹ 48 more rows

Summary statistics

Summary statistics for dependent variable posttest

df58 %>% group_by(VLEs) %>%  get_summary_stats(posttest, type="common")
## # A tibble: 3 × 11
##   VLEs  variable     n   min   max median   iqr  mean    sd    se    ci
##   <chr> <chr>    <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 VLE1  posttest    20  40    90.6   55.3  17.8  57.6  13.3  2.98  6.23
## 2 VLE2  posttest    19  23.3  98.3   54.4  10    56.5  20.8  4.76 10.0 
## 3 VLE3  posttest    19  57.2  98.3   77.8  24.4  77.5  13.5  3.11  6.53

Summary statistics for covariate pretest

df58 %>% group_by(VLEs) %>%  get_summary_stats(pretest, type="common")
## # A tibble: 3 × 11
##   VLEs  variable     n   min   max median   iqr  mean    sd    se    ci
##   <chr> <chr>    <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 VLE1  pretest     20  27.8  63.9   42.2  14.3  44.5  9.82  2.19  4.59
## 2 VLE2  pretest     19  27.8  95     48.3  23.9  54.2 22.4   5.13 10.8 
## 3 VLE3  pretest     19  36.1  91.7   80    27.8  70.5 17.8   4.07  8.56

Data visualization

p1 <- ggplot(df58, aes(pretest, posttest, colour = VLEs)) + geom_point(size = 3) + theme(legend.position="top")
p2 <- ggplot(df58, aes(x = VLEs, y = posttest, col = VLEs)) + geom_boxplot(outlier.shape = NA) + geom_jitter(width = 0.2) + theme(legend.position="top")
p3 <- ggplot(df58, aes(x = VLEs, y = pretest, fill = VLEs)) + geom_boxplot(outlier.shape = NA) + geom_jitter(width = 0.2) + theme(legend.position="top")
grid.arrange(p1, p2, p3, ncol=3)

Model fit

Fit the model, the covariate goes first

model <- lm(posttest ~ pretest + VLEs, data = df58)

Inspect the model diagnostic metrics and remove details

model.metrics <- augment(model) %>%
  select(-.hat, -.sigma, -.fitted)
head(model.metrics, 3)
## # A tibble: 3 × 6
##   posttest pretest VLEs  .resid .cooksd .std.resid
##      <dbl>   <dbl> <chr>  <dbl>   <dbl>      <dbl>
## 1     48.3    43.3 VLE1   -8.50 0.00811     -0.784
## 2     90.6    48.3 VLE1   30.3  0.105        2.80 
## 3     55      51.7 VLE1   -7.55 0.00683     -0.698

Assumption tests

Outliers above 3.00

model.metrics %>% 
  dplyr::filter(abs(.std.resid) > 3) %>%
  as.data.frame()
## [1] posttest   pretest    VLEs       .resid     .cooksd    .std.resid
## <0 rows> (or 0-length row.names)

Outliers above 2.50

model.metrics %>% 
  dplyr::filter(abs(.std.resid) > 2.5) %>%
  as.data.frame()
##   posttest pretest VLEs  .resid   .cooksd .std.resid
## 1    90.56   48.33 VLE1 30.3035 0.1049588   2.797879

Linearity

ggscatter(
  df58, x = "pretest", y = "posttest",
  color = "VLEs", add = "reg.line"
)+
  stat_regline_equation(
    aes(label =  paste(..eq.label.., ..rr.label.., sep = "~~~~"), color = VLEs)
  )

Homogeneity of regression slopes

df58 %>% anova_test(posttest ~ VLEs*pretest)
## ANOVA Table (type II tests)
## 
##         Effect DFn DFd      F        p p<.05   ges
## 1         VLEs   2  52  3.881 2.70e-02     * 0.130
## 2      pretest   1  52 60.491 2.86e-10     * 0.538
## 3 VLEs:pretest   2  52  0.028 9.73e-01       0.001

Normality of residuals

shapiro_test(model.metrics$.resid)
## # A tibble: 1 × 3
##   variable             statistic p.value
##   <chr>                    <dbl>   <dbl>
## 1 model.metrics$.resid     0.990   0.901

Homogeneity of variances

model.metrics %>% levene_test(.resid ~ VLEs)
## # A tibble: 1 × 4
##     df1   df2 statistic      p
##   <int> <int>     <dbl>  <dbl>
## 1     2    55      2.86 0.0656

ANCOVA

ANCOVA table

anova_test(data = df58, formula = posttest ~ pretest + VLEs, type = 3, detailed = TRUE)
## ANOVA Table (type III tests)
## 
##        Effect      SSn      SSd DFn DFd      F        p p<.05   ges
## 1 (Intercept) 3027.477 6674.372   1  54 24.494 7.67e-06     * 0.312
## 2     pretest 7755.854 6674.372   1  54 62.750 1.32e-10     * 0.537
## 3        VLEs  995.137 6674.372   2  54  4.026 2.30e-02     * 0.130

Estimated Marginal Means (EMM)

adj_means <- emmeans_test(data = df58, formula = posttest ~ VLEs, covariate = pretest)
get_emmeans(adj_means)
## # A tibble: 3 × 8
##   pretest VLEs  emmean    se    df conf.low conf.high method      
##     <dbl> <fct>  <dbl> <dbl> <dbl>    <dbl>     <dbl> <chr>       
## 1    56.2 VLE1    65.6  2.68    54     60.3      71.0 Emmeans test
## 2    56.2 VLE2    57.8  2.56    54     52.7      63.0 Emmeans test
## 3    56.2 VLE3    67.7  2.83    54     62.0      73.3 Emmeans test

Post-hoc test: bonferroni

emmeans_test(data = df58, formula = posttest ~ VLEs, covariate = pretest, p.adjust.method = "bonferroni")
## # A tibble: 3 × 9
##   term         .y.      group1 group2    df statistic      p  p.adj p.adj.signif
## * <chr>        <chr>    <chr>  <chr>  <dbl>     <dbl>  <dbl>  <dbl> <chr>       
## 1 pretest*VLEs posttest VLE1   VLE2      54     2.13  0.0377 0.113  ns          
## 2 pretest*VLEs posttest VLE1   VLE3      54    -0.481 0.632  1      ns          
## 3 pretest*VLEs posttest VLE2   VLE3      54    -2.54  0.0141 0.0423 *

Referensi

  1. Ancova in R
  2. Ancova using R and Python

Sitasi dokumen ini

Santyadiputra, G. S., Purnomo, & Juniantari, M. (2023, November 29). Ancova satu jalur menggunakan R. RPubs. https://doi.org/10.5281/zenodo.10219200

DOI