P15 - Model Moderator Berganda

Library

library(foreign)
library(psych)
library(lm.beta)
library(ggplot2)
library(dplyr)

Dataset

d1 <- read.spss(file ="D:/UNY/MySta/SEM 5/AMS/dataset_ams/multiple_moderators.sav", to.data.frame=TRUE)

## re-encoding from CP1252

str(d1)

## 'data.frame':    1774 obs. of  4 variables:
##  $ age     : Factor w/ 9 levels "9","10","11",..: 4 2 4 3 2 2 3 2 2 2 ...
##  $ rum     : num  NA 1 2 1.25 2 2 2 1 2 3.75 ...
##  $ stress  : num  1 1 3.75 2.25 1.75 1.5 1.75 2 1.75 3.25 ...
##  $ selfharm: num  1 1 1 1 1 1 1.5 1 1 1 ...
##  - attr(*, "variable.labels")= Named chr [1:4] "Age" "Rumination" "Stress" "Self-harm"
##   ..- attr(*, "names")= chr [1:4] "age" "rum" "stress" "selfharm"
##  - attr(*, "codepage")= int 1252

d1$agec <- ifelse(d1$age == "13"| d1$age == "14"| d1$age == "15"| d1$age == "16"| d1$age == "17", 1, 0)
head(d1)

##   age  rum stress selfharm agec
## 1  12   NA   1.00        1    0
## 2  10 1.00   1.00        1    0
## 3  12 2.00   3.75        1    0
## 4  11 1.25   2.25        1    0
## 5  10 2.00   1.75        1    0
## 6  10 2.00   1.50        1    0

Statistik Deskriptif

des <- describe(d1)
print(des, digits = 4)

##          vars    n   mean     sd median trimmed    mad min max range   skew
## age*        1 1774 4.1246 1.7328   4.00  4.0394 2.9652   1   8     7 0.1662
## rum         2 1757 2.3775 0.8364   2.25  2.3316 0.7413   1   5     4 0.5037
## stress      3 1732 2.2823 0.9186   2.00  2.1890 0.7413   1   5     4 0.8542
## selfharm    4 1750 1.1849 0.5444   1.00  1.0407 0.0000   1   5     4 4.0808
## agec        5 1774 0.3963 0.4893   0.00  0.3704 0.0000   0   1     1 0.4237
##          kurtosis     se
## age*      -1.2608 0.0411
## rum       -0.0927 0.0200
## stress     0.3153 0.0221
## selfharm  19.4117 0.0130
## agec      -1.8215 0.0116

Model 1

model1 <- lm(selfharm ~ rum + stress + agec, data = d1)
summary(lm.beta(model1))

## 
## Call:
## lm(formula = selfharm ~ rum + stress + agec, data = d1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.8393 -0.2231 -0.0855  0.0489  3.6273 
## 
## Coefficients:
##             Estimate Standardized Std. Error t value Pr(>|t|)    
## (Intercept)  0.54026           NA    0.04124  13.100  < 2e-16 ***
## rum          0.12078      0.18861    0.01590   7.598 4.96e-14 ***
## stress       0.13903      0.23771    0.01454   9.564  < 2e-16 ***
## agec         0.09365      0.08548    0.02471   3.789 0.000156 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4994 on 1697 degrees of freedom
##   (73 observations deleted due to missingness)
## Multiple R-squared:  0.139,  Adjusted R-squared:  0.1375 
## F-statistic: 91.35 on 3 and 1697 DF,  p-value: < 2.2e-16

with(d1,cor.test(rum,selfharm))

## 
##  Pearson's product-moment correlation
## 
## data:  rum and selfharm
## t = 12.559, df = 1735, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.2449768 0.3312195
## sample estimates:
##       cor 
## 0.2886837

Model 2

model2 <- lm(selfharm ~ rum + stress + agec + rum*agec + rum*stress + agec*stress, data = d1)
summary(lm.beta(model2))

## 
## Call:
## lm(formula = selfharm ~ rum + stress + agec + rum * agec + rum * 
##     stress + agec * stress, data = d1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.4103 -0.1799 -0.0749 -0.0076  3.7317 
## 
## Coefficients:
##             Estimate Standardized Std. Error t value Pr(>|t|)    
## (Intercept)  1.08623           NA    0.09686  11.215  < 2e-16 ***
## rum         -0.07666     -0.11971    0.03896  -1.968 0.049257 *  
## stress      -0.02218     -0.03793    0.04082  -0.543 0.586893    
## agec        -0.52185     -0.47632    0.08065  -6.471 1.27e-10 ***
## rum:agec     0.18737      0.44528    0.03244   5.775 9.12e-09 ***
## rum:stress   0.05055      0.35941    0.01422   3.556 0.000387 ***
## stress:agec  0.07352      0.17479    0.02988   2.461 0.013965 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4881 on 1694 degrees of freedom
##   (73 observations deleted due to missingness)
## Multiple R-squared:  0.179,  Adjusted R-squared:  0.1761 
## F-statistic: 61.55 on 6 and 1694 DF,  p-value: < 2.2e-16

Model 3

model3 <- lm(selfharm ~ rum + stress + agec + rum*agec + rum*stress + agec*stress + agec*stress*rum, data = d1)
summary(lm.beta(model3))

## 
## Call:
## lm(formula = selfharm ~ rum + stress + agec + rum * agec + rum * 
##     stress + agec * stress + agec * stress * rum, data = d1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.8189 -0.1763 -0.0797  0.0046  3.7165 
## 
## Coefficients:
##                 Estimate Standardized Std. Error t value Pr(>|t|)    
## (Intercept)      0.72583           NA    0.11408   6.362 2.55e-10 ***
## rum              0.07036      0.10987    0.04608   1.527  0.12695    
## stress           0.13848      0.23676    0.04891   2.831  0.00469 ** 
## agec             0.50637      0.46219    0.19344   2.618  0.00893 ** 
## rum:agec        -0.23774     -0.56498    0.07961  -2.986  0.00287 ** 
## rum:stress      -0.01147     -0.08152    0.01764  -0.650  0.51577    
## stress:agec     -0.36531     -0.86850    0.08080  -4.521 6.58e-06 ***
## rum:stress:agec  0.17088      1.20519    0.02928   5.836 6.39e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4834 on 1693 degrees of freedom
##   (73 observations deleted due to missingness)
## Multiple R-squared:  0.1952, Adjusted R-squared:  0.1919 
## F-statistic: 58.65 on 7 and 1693 DF,  p-value: < 2.2e-16

Output statistik regresi untuk grup remaja lebih muda

# Young Model 1
young <- subset(d1, agec == "0")
young.mod1 <- lm(selfharm ~ stress + rum, data = young)
summary(lm.beta(young.mod1))

## 
## Call:
## lm(formula = selfharm ~ stress + rum, data = young)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.5522 -0.1667 -0.0856 -0.0148  3.7137 
## 
## Coefficients:
##             Estimate Standardized Std. Error t value Pr(>|t|)    
## (Intercept)  0.79246           NA    0.04533  17.483  < 2e-16 ***
## stress       0.10878      0.22582    0.01579   6.890 9.79e-12 ***
## rum          0.04318      0.08077    0.01752   2.464   0.0139 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4376 on 1011 degrees of freedom
##   (57 observations deleted due to missingness)
## Multiple R-squared:  0.0714, Adjusted R-squared:  0.06956 
## F-statistic: 38.87 on 2 and 1011 DF,  p-value: < 2.2e-16

# Young Model 2
young.mod2 <- lm(selfharm ~ stress + rum + stress*rum, data = young)
summary(lm.beta(young.mod2))

## 
## Call:
## lm(formula = selfharm ~ stress + rum + stress * rum, data = young)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.4834 -0.1724 -0.0894 -0.0110  3.7165 
## 
## Coefficients:
##             Estimate Standardized Std. Error t value Pr(>|t|)    
## (Intercept)  0.72583           NA    0.10330   7.026  3.9e-12 ***
## stress       0.13848      0.28748    0.04429   3.127  0.00182 ** 
## rum          0.07036      0.13161    0.04172   1.686  0.09204 .  
## stress:rum  -0.01147     -0.09617    0.01597  -0.718  0.47301    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4377 on 1010 degrees of freedom
##   (57 observations deleted due to missingness)
## Multiple R-squared:  0.07187,    Adjusted R-squared:  0.06912 
## F-statistic: 26.07 on 3 and 1010 DF,  p-value: 3.018e-16

Output statistik regresi untuk grup remaja lebih tua

# Old Model 1
old <- subset(d1, agec == "1")
old.mod1 <- lm(selfharm ~ stress + rum, data = old)
summary(lm.beta(old.mod1))

## 
## Call:
## lm(formula = selfharm ~ stress + rum, data = old)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.2206 -0.2748 -0.0801  0.1300  3.3500 
## 
## Coefficients:
##             Estimate Standardized Std. Error t value Pr(>|t|)    
## (Intercept)  0.25430           NA    0.07164   3.550 0.000412 ***
## stress       0.18049      0.25018    0.02759   6.543 1.18e-10 ***
## rum          0.23983      0.31134    0.02946   8.142 1.84e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.558 on 684 degrees of freedom
##   (16 observations deleted due to missingness)
## Multiple R-squared:  0.2351, Adjusted R-squared:  0.2328 
## F-statistic: 105.1 on 2 and 684 DF,  p-value: < 2.2e-16

# Old Model 2
old.mod2 <- lm(selfharm ~ stress + rum + stress*rum, data = old)
summary(lm.beta(old.mod2))

## 
## Call:
## lm(formula = selfharm ~ stress + rum + stress * rum, data = old)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.8189 -0.1950 -0.0534  0.0184  3.3753 
## 
## Coefficients:
##             Estimate Standardized Std. Error t value Pr(>|t|)    
## (Intercept)  1.23219           NA    0.17579   7.009 5.76e-12 ***
## stress      -0.22682     -0.31440    0.07238  -3.134   0.0018 ** 
## rum         -0.16738     -0.21728    0.07306  -2.291   0.0223 *  
## stress:rum   0.15941      0.96264    0.02630   6.062 2.23e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.544 on 683 degrees of freedom
##   (16 observations deleted due to missingness)
## Multiple R-squared:  0.2741, Adjusted R-squared:  0.2709 
## F-statistic: 85.97 on 3 and 683 DF,  p-value: < 2.2e-16

Grup lebih muda

Hitung rata-rata dari pengaruh interaksi stress dan perenungan pada menyakiti diri sendiri bagi grup lebih muda

#  PLOT MODERASI UNTUK GRUP LEBIH MUDA


# Stress value low/med/high (±1 SD)
stress.mean <- mean(young$stress, na.rm = TRUE)
stress.sd   <- sd(young$stress, na.rm = TRUE)

stress.levels <- data.frame(
  stress = c(stress.mean - stress.sd,
             stress.mean,
             stress.mean + stress.sd),
  StressCat = c("low", "med", "high")
)

# level rumination low/med/high
rum.mean <- mean(young$rum, na.rm = TRUE)
rum.sd   <- sd(young$rum, na.rm = TRUE)

rum.levels <- data.frame(
  rum = c(rum.mean - rum.sd,
          rum.mean,
          rum.mean + rum.sd),
  RumCat = c("low", "med", "high")
)

# Kombinasi level stress × rumination
pred.grid <- merge(stress.levels, rum.levels)

# Prediksi self-harm berdasarkan model
pred.grid$pred <- predict(young.mod2, newdata = pred.grid)

# kategori terurut
pred.grid$StressCat <- factor(pred.grid$StressCat, levels = c("low","med","high"))
pred.grid$RumCat    <- factor(pred.grid$RumCat, levels = c("low","med","high"))


# PLOT
ggplot(pred.grid, aes(x = StressCat, y = pred, 
                      group = RumCat, shape = RumCat)) +
  geom_line() +
  geom_point(size = 3) +
  scale_y_continuous(limits = c(0, 1.4), breaks = seq(0, 1.4, by = 0.2)) +
  labs(
    title = "Moderation by Rumination",
    x = "Stress",
    y = "Self-harm",
    shape = "Rumination"
  ) +
  theme_minimal(base_size = 14) +
  theme(
    plot.title = element_text(hjust = 0.5), 
    legend.position = c(0.85, 0.25),      
    legend.background = element_rect(
      fill = "white", color = "black", size = 0.3 
    )
  )

## Warning: The `size` argument of `element_rect()` is deprecated as of ggplot2 3.4.0.
## ℹ Please use the `linewidth` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: A numeric `legend.position` argument in `theme()` was deprecated in ggplot2
## 3.5.0.
## ℹ Please use the `legend.position.inside` argument of `theme()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Interpretasi: Pengaruh moderasi tidak signifikan sehingga terlihat garis-garis adalah parallel.

Grup Lebih Tua

Hitung rata-rata dari pengaruh interaksi stress dan perenungan pada menyakiti diri sendiri bagi grup lebih tua

# Stress value low/med/high (±1 SD)
stress.mean.old <- mean(old$stress, na.rm = TRUE)
stress.sd.old   <- sd(old$stress, na.rm = TRUE)

stress.levels.old <- data.frame(
  stress = c(stress.mean.old - stress.sd.old,
             stress.mean.old,
             stress.mean.old + stress.sd.old),
  StressCat = c("low", "med", "high")
)

# level rumination low/med/high
rum.mean.old <- mean(old$rum, na.rm = TRUE)
rum.sd.old   <- sd(old$rum, na.rm = TRUE)

rum.levels.old <- data.frame(
  rum = c(rum.mean.old - rum.sd.old,
          rum.mean.old,
          rum.mean.old + rum.sd.old),
  RumCat = c("low", "med", "high")
)

# Kombinasi grid prediksi
pred.grid.old <- merge(stress.levels.old, rum.levels.old)
pred.grid.old$pred <- predict(old.mod2, newdata = pred.grid.old)

# kategori terurut
pred.grid.old$StressCat <- factor(pred.grid.old$StressCat, 
                                  levels = c("low","med","high"))
pred.grid.old$RumCat <- factor(pred.grid.old$RumCat, 
                               levels = c("low","med","high"))

# PLOT
ggplot(pred.grid.old, 
       aes(x = StressCat, y = pred, 
           group = RumCat, shape = RumCat)) +
  
  geom_line() +
  geom_point(size = 3) +
  scale_y_continuous(limits = c(0.8, 1.7),
                     breaks = seq(0.8, 1.7, by = 0.1)) +
  
  labs(
    title = "Moderation by Rumination",
    x = "Stress",
    y = "Self harm",
    shape = "Rumination"
  ) +
  
  theme_minimal(base_size = 14) +
  theme(
    plot.title = element_text(hjust = 0.5), 
    legend.position = c(0.25, 0.75),
    legend.background = element_rect(
      fill = "white", color = "black", size = 0.4
    )
  )

Interpretasi:

Gambar ini menunjukkan dinamika di mana remaja yang lebih tua yang merenung pada tingkat tinggi memperlihatkan hubungan yang lebih kuat antara stres dan menyakiti diri sendiri daripada mereka yang kurang merenung.
Analisis slope sederhana akan menunjukkan bahwa slope untuk grup perenungan sedang dan tinggi secara statistik berbeda dari nol, tetapi grup perenungan rendah tampaknya menunjukkan slope yang tidak signifikan.