Lineární modelování: vzor úkolu

Zadání

V tomto domácím úkolu využijete předdefinovaná data a vytvoříte z nich jednoduchý lineární model. Opět jde o dataframe sestavený z různých podmínek,

Balíčky

Do knihovny si připravíme následující balíčky. Cokoliv z toho nemáte, to si předtím samozřejmě ještě nainstalujte.

library(plyr)
library(dplyr)

Attaching package: 'dplyr'

The following objects are masked from 'package:plyr':

    arrange, count, desc, failwith, id, mutate, rename, summarise,
    summarize

The following objects are masked from 'package:stats':

    filter, lag

The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union

library(ggplot2)
library(lmerTest)

Warning: package 'lmerTest' was built under R version 4.3.3

Loading required package: lme4

Loading required package: Matrix

Attaching package: 'lmerTest'

The following object is masked from 'package:lme4':

    lmer

The following object is masked from 'package:stats':

    step

Tvorba dat

To už známe z minula – překopírujte tento kód tak, jak je. Vytvoří to 24 vektorů (podmínek pomyslného kvantitativního experimentu) složených z různých čísel tak, že každá podmínka je v normální distribuci a má jiný průměr než ostatní. Tentokrát však vytváříme nejen sadu vektorů od aa do ax, ale také od Aa do Ax.

set.seed(123)
x <- rnorm(100, mean = 3, sd = 1)
for (i in 1:24) {
  assign(paste0("a", letters[i]), rnorm(100, mean = 0.5 + 0.1*i, sd = 1))
}
for (i in 1:24) {
  assign(paste0("A", letters[i]), rnorm(100, mean = 3 + 0.1*i, sd = 1))
}

Tímto se vytvoří tabulka, jejíž součástí jsou podmínky aa-ax, Aa-Ax i podmínka x:

df <- data.frame()

for (i in 1:24) {
  df <- rbind(df, data.frame(condition = paste0("a", letters[i]), value = get(paste0("a", letters[i]))))}

for (i in 1:24) {
  df <- rbind(df, data.frame(condition = paste0("A", letters[i]), value = get(paste0("A", letters[i]))))}

df <- rbind(df, data.frame(condition = "x", value = x))

Na grafu vidíme, že podmínky na A mají pravidelný rozestup od podmínek na a.

boxplot(value ~ condition, data = df, col = "lightblue", main = "Boxplot of the conditions", xlab = "condition", ylab = "value")

Váš úkol

Teď budete každý postupovat podle svého zadání. Stejně jako v úkolu s t-testem, i tentokrát má každý student podle rozpisu přidělenou podmínku. Jelikož máme dva sety podmínek, budete mít jednu od “a” a jednu od “A” (například aa a Aa). Následující kódování tomu přizpůsobte, podmínka “x” však je referenční a zůstane stejná.

1. nadefinovat nový dataframe vyfiltrováním “x” a vašich dalších dvou podmínek z původního dataframu:

df2 <- df %>%
  filter(condition %in% c("x", "aa", "Aa"))
head(df2)

  condition      value
1        aa -0.1104066
2        aa  0.8568837
3        aa  0.3533081
4        aa  0.2524574
5        aa -0.3516186
6        aa  0.5549723

2. znázornit ho histogramem:

ggplot(df2, aes(x = value, fill = condition)) +
  geom_histogram(bins = 20, alpha = 0.5, position = "identity") +
  labs(title = "Histogram of x, aa and Aa", x = "value", y = "count") +
  theme_minimal()

3. vygenerovat z něj lineární model:

as.factor(df2$condition) -> df2$condition
df2$condition <- relevel(df2$condition, ref = "x")

model <- lmer(value ~ condition + (1|condition), data = df2, REML = FALSE)

boundary (singular) fit: see help('isSingular')

summary(model)

Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
  method [lmerModLmerTest]
Formula: value ~ condition + (1 | condition)
   Data: df2

     AIC      BIC   logLik deviance df.resid 
   828.9    847.4   -409.4    818.9      295 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.6178 -0.7174  0.0044  0.6981  3.5350 

Random effects:
 Groups    Name        Variance Std.Dev.
 condition (Intercept) 0.0000   0.0000  
 Residual              0.8973   0.9473  
Number of obs: 300, groups:  condition, 3

Fixed effects:
             Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)   3.09041    0.09473 300.00000  32.624   <2e-16 ***
conditionaa  -2.59795    0.13396 300.00000 -19.393   <2e-16 ***
conditionAa   0.11568    0.13396 300.00000   0.864    0.389    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) condtn
conditionaa -0.707       
conditionAa -0.707  0.500
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')

Příkazem “relevel(df2$condition, ref =”x”)” jsem si stanovila, že “x” je baseline/referenční podmínka modelu, tedy ta, se kterou model porovná ty druhé dvě. V sekci “fixed effects” ji najdu jako “(Intercept)”.

O těch dvou ostatních lze říct na základě p-value (Pr), že
1.: rozdíl mezi “x” a “aa” je statisticky významný, protože u “aa” máme p-hodnotu velmi nízkou (na mínus šestnáctou).

2.: rozdíl mezi “x” a “Aa” není statisticky významný, p-hodnota je totiž 0.389 (38,9 % pravděpodobnost platnosti nulové hypotézy), což je více než konvenčně stanovených 0.05 (5%).

Kdyby toto byly výsledky experimentu přijatelnosti v L-rexu, mohli bychom čekat, že úpravou dat “x” podle podmínky “aa” získáme nižší hodnocení od účastníků (sloupec “Estimate” ukazuje číslo v mínusu, takže jde o snížení čísel oproti baseline). Podmínka “Aa” by vzhledem ke statistické nesignifikantnosti buď neměla, nebo měla nedostatečně prokázaný efekt na číselné hodnoty hodnocení.