ADCS

まず、df3 が巨大なデータフレームなので、ADCSとその他必要なデータだけの小さなデータフレームを作成します。

最初に、行数をそろえるために、インデックス列をコピーします。最初に使う列をすべて指定してもよいのですが、ここではインデックス列のみとしました。なお、最初にすべてコピーする場合は、

df3ADCS <- df3[c("Patientennummer", "Gruppe", "Messzeitpunkt_Nummer", "ADCS_ges")]

とします。列数が一つの時は文字列“A”ですが、複数の場合は文字列リストc(“A”, “B”, “C”) になります。

Time (Messzeitpunkt_Nummer) は、因子型にもとれますが、geom_smooth() を使うためには x軸も数値型である必要があるので、数値型とします。

library(ggplot2)

df3ADCS <- df3["Patientennummer"]
df3ADCS$Group <- factor(df3$Gruppe, level=c(1,2))
df3ADCS$Time <- as.numeric(df3$Messzeitpunkt_Nummer)
df3ADCS$ADCS <- as.numeric(df3$ADCS_ges)

平均値と標準誤差を手動で計算する方法です。

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

df3Plot <- df3ADCS %>%
              group_by(Group, Time) %>%
              summarise(ADCS_means = mean(ADCS, na.rm=T),
                        ADCS_SEs = sd(ADCS, na.rm=T)/sqrt(length(ADCS)))

## `summarise()` has grouped output by 'Group'. You can override using the `.groups` argument.

ggplot(data = df3Plot, aes(x = Time, y = ADCS_means, group = Group, colour = Group)) +
  geom_point(position = position_dodge(width=0.2)) +
  geom_line(position = position_dodge(width=0.2)) + 
  geom_errorbar(aes(ymin=ADCS_means-ADCS_SEs, ymax=ADCS_means+ADCS_SEs), width=.1, position = position_dodge(width=0.2))+
  theme_bw()

平均値と標準誤差を自動的にプロットする方法です。

library(ggpubr)
ggline(df3ADCS, x = "Time", y = "ADCS", add = "mean_se",
          color = "Group")+
  stat_compare_means(aes(group = Group), label = "p.signif")

## Warning: Removed 4 rows containing non-finite values (stat_summary).

## Warning: Removed 4 rows containing non-finite values (stat_compare_means).

ggplot(data = df3ADCS, aes(x = Time, y = ADCS, colour = Group)) +
  geom_smooth(method="lm")

## `geom_smooth()` using formula 'y ~ x'

## Warning: Removed 4 rows containing non-finite values (stat_smooth).

two way repeated measures ANOVA

データのない行を削除し、 外れ値を確認します。

library(rstatix)

## 
## Attaching package: 'rstatix'

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

df3ADCSnoNA <- na.omit(df3ADCS)
identify_outliers(group_by(df3ADCSnoNA, "Group", "Time"), "ADCS")

df3ADCSnoNA <- na.omit(df3ADCS)
shapiro_test(group_by(df3ADCSnoNA, "Group", "Time"), "ADCS")

残念ながら、これは以下のエラーが出ます。

Error: Problem with mutate() column data. i data = map(.data$data, .f, ...). x Can’t subset columns that don’t exist. x Column "ADCS" doesn’t exist. Run rlang::last_error() to see where the error occurred.

library(rstatix)
aovDf3ADCS <- anova_test(
  data = df3ADCSnoNA, 
  dv = "ADCS", 
  wid = "Patientennummer",
  within = c("Group","Time"))
get_anova_table(aovDf3ADCS)

残念ながら、これは以下のエラーが出ます。

Error in lm.fit(x, y, offset = offset, singular.ok = singular.ok, …) : 0 (non-NA) cases

repeated measures ANOVA は、計測されていないデータがあると実行することができません。 例えば、Patientennummber = 11 の T3 におけるデータがありません（他にも何か所か欠測があります）。 このため、Patientennumber = 11 全体を削除するか、代入 (imputation) を必要があります。

介入後のフォローアップを伴う解析を予定している場合は、特に欠測が発生しやすいので注意が必要です。

削除するか、代入するかは、研究計画段階で決める必要があります。

Mixed model

一方、欠測値があっても実行できるのが　Mixed Model の強みです。

library(lme4)
library(sjPlot)
lmeModeldf3ADCS = lmer(ADCS ~ Group*Time + (1|Patientennummer), data=df3ADCS)
anova(lmeModeldf3ADCS)

tab_model(lmeModeldf3ADCS)

library(emmeans)
emmeans(lmeModeldf3ADCS, list(pairwise ~ Group), adjust = "tukey")

## $`emmeans of Group`
##  Group emmean   SE   df lower.CL upper.CL
##  1       62.0 2.74 28.1     56.4     67.6
##  2       56.3 2.74 27.9     50.7     61.9
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## 
## $`pairwise differences of Group`
##  1     estimate   SE df t.ratio p.value
##  1 - 2     5.69 3.87 28   1.469  0.1529
## 
## Degrees-of-freedom method: kenward-roger

ここでは、 p = 0.1529 より、交互作用があったとは言えないと解釈します。

NPI

library(ggplot2)


df3NPI <- df3["Patientennummer"]

df3NPI$Group <- factor(df3$Gruppe, level=c(1,2))
df3NPI$Time <- as.numeric(df3$Messzeitpunkt_Nummer)
df3NPI$NPI <- as.numeric(df3$NPI_ges)
  
#df3NPI <- subset(df3, !is.na(df3$NPI))

ggplot(data = df3NPI, aes(x = Time, y = NPI, colour = Group)) +
  geom_point() +
  geom_smooth(method = "lm")

library(lme4)
library(sjPlot)
lmeModeldf3NPI = lmer(NPI ~ Group*Time + (1|Patientennummer), data=df3NPI)
anova(lmeModeldf3NPI)

tab_model(lmeModeldf3NPI)

library(emmeans)
emmeans(lmeModeldf3NPI, list(pairwise ~ Group), adjust = "tukey")

## $`emmeans of Group`
##  Group emmean   SE   df lower.CL upper.CL
##  1       8.36 1.92 28.2     4.43     12.3
##  2      16.61 1.91 27.8    12.69     20.5
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## 
## $`pairwise differences of Group`
##  1     estimate   SE df t.ratio p.value
##  1 - 2    -8.24 2.71 28  -3.044  0.0050
## 
## Degrees-of-freedom method: kenward-roger

p = 0.0050 なので、交互作用があったと解釈します。

サンプルサイズ

library(WebPower)

## Loading required package: MASS

## 
## Attaching package: 'MASS'

## The following object is masked from 'package:rstatix':
## 
##     select

## The following object is masked from 'package:dplyr':
## 
##     select

## Loading required package: lavaan

## This is lavaan 0.6-9
## lavaan is FREE software! Please report any bugs.

## Loading required package: parallel

## Loading required package: PearsonDS

wp.rmanova(ng = 2,
           nm = 3,
           power = 0.80,
           alpha = 0.05,
           f = 0.32,
           type = "2")

## Repeated-measures ANOVA analysis
## 
##            n    f ng nm nscor alpha power
##     95.61069 0.32  2  3     1  0.05   0.8
## 
## NOTE: Power analysis for interaction-effect test
## URL: http://psychstat.org/rmanova