Сопоставление функции каждому участнику исследования

Я пытался применить функцию к этому набору данных, который состоит из данных из отдельных исследований, разработанных AB. Он состоит из 6 переменных: a. Уровень (количество циклов AB в этом конкретном наборе данных всегда равно 1, поскольку все исследования разработаны AB); б. ID (идентификационный код участника), c. Шкала (тип шкалы, по которой оценивается участник), d. Время (случай рейтинга), т.е. Фаза (A = исходный уровень, B = вмешательство), f. Score (счет участника).

library(purrr)
library(dplyr)

Tier <- rep(c(1), 36)
ID <- c("C1", "C1", "C1", "C1", "C1", "C1", "C1", "C1", "C1", "C1", "C1", "C2", "C2", "C2", "C2", "C2", "C2", "C2", "C2", "C1", "C1", "C1", "C1", "C1", "C1", "C1", "C1", "C2", "C2", "C2", "C2", "C2", "C2", "C2", "C2", "C2")
Scale <- rep(c(1, 2), each = 18)
Time <- c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 1, 2, 3, 4, 5, 6, 7, 8, 1, 2, 3, 4, 5, 6, 7, 8, 1, 2, 3, 4, 5, 6, 7, 8, 9)
Phase <- c("A", "A", "A", "A", "A", "B", "B", "B", "B", "B", "B", "A", "A", "A", "A", "B", "B", "B", "B", "A", "A", "A", "A", "B", "B", "B", "B", "A", "A", "A", "A", "B", "B", "B", "B", "B")
Score <- c(4, 8, 10, 12, 15, 7, 7, 9, 14, 15, 16, 4, 3, 2, 2, 7, 7, 9, 14, 2, 3, 6, 6, 7, 5, 9, 11, 5, 6, 3, 4, 8, 7, 9, 3, 3)
db <- data.frame(Scale, ID, Time, Phase, Score, Tier)

Это две функции, разработанные Manolov and Rochat (2015). Они отлично работают до тех пор, пока не дойдет до расчета более одного участника за один раз. Поэтому я попытался изменить некоторые детали второго (MPD_meta), чтобы сделать это.

MPD <- function(a, b) { #Manolov & Rochat, 2015

# Obtain phase length
n_a <- length(a[!is.na(a)])
n_b <- length(b[!is.na(b)])
meanA <- mean(a[!is.na(a)]) 

# Estimate baseline trend
base_diff <- rep(0,(n_a - 1))

for (i in 1:(n_a - 1)) base_diff[i] <- a[!is.na(a)][i + 1] - a[!is.na(a)][i]

trendA <- mean(base_diff)

# Predict baseline data
baseline_pred <- rep(0, n_a)

midX <- median(c(1:n_a))
midY <- median(a[!is.na(a)])

is.wholenumber <- function(x, tol = .Machine$double.eps ^ 0.5)  abs(x - round(x)) < tol

if (is.wholenumber(midX)) {
baseline_pred[midX] <- midY
for (i in (midX - 1): 1) baseline_pred[i] <- baseline_pred[i + 1] - trendA
for (i in (midX + 1): n_a) baseline_pred[i] <- baseline_pred[i - 1] + trendA
}

#baseline_pred

if (!is.wholenumber(midX)) {
baseline_pred[midX - 0.5] <- midY - (0.5 * trendA)
baseline_pred[midX + 0.5] <- midY + (0.5 * trendA)
for (i in (midX - 1.5) : 1) baseline_pred[i] <- baseline_pred[i + 1] - trendA
for (i in (midX + 1.5) : n_a) baseline_pred[i] <- baseline_pred[i - 1] + trendA
}


# Project baseline trend to the treatment phase
treatment_pred <- rep(0, n_b)
treatment_pred[1] <- baseline_pred[n_a] + trendA

for (i in 2:n_b) 
    treatment_pred[i] <- treatment_pred[i - 1] + trendA

diff_absolute <- rep(0, n_b)
diff_relative <- rep(0, n_b)

for (i in 1:n_b) {
diff_absolute[i] <- b[!is.na(b)][i] - treatment_pred[i]

if (treatment_pred[i]!=0) diff_relative[i] <- ((b[!is.na(b)][i] - treatment_pred[i]) / abs(treatment_pred[i])) * 100
    }

mpd_value <- mean(diff_absolute)
mpd_relative <- mean(diff_relative)
mpd_relsd <- mpd_value/sd(a[!is.na(a)])

# Compute the residual (MSE): difference between actual and predicted baseline data
residual <- 0

for (i in 1:n_a)

  residual <- residual + (a[!is.na(a)][i] - baseline_pred[i]) * (a[!is.na(a)][i] - baseline_pred[i])
  residual_mse <- residual/n_a

list(MPD_raw = mpd_value, Residual_MSE = residual_mse, MPD_sd = mpd_relsd) #, MPD_percent = mpd_relative, Baseline_pred = baseline_pred

}

MPD_meta <- function(DBR, n) {

MBD <- db %>%
  dplyr::filter(Scale == DBR, ID == n)

# Create the objects necessary for the calculations and the graphs
# Dimensions of the data set
tiers <- max(MBD$Tier)
max.length <- max(MBD$Time)
max.score <- max(MBD$Score)
min.score <- min(MBD$Score)

# Vectors for keeping the main results
# Series lengths per tier
obs <- rep(0,tiers)
# Raw MPD effect sizes per tier
results <- rep(0,tiers)
# Weights per tier
tier.weights <- rep(0,tiers)
# Standardized MPD results per tier
stdresults <- rep(0,tiers)

# Obtain the values calling the MPD function
for (i in 1:tiers) {

  tempo <- subset(MBD, Tier == i)
  # Number of observations for the tier
  obs[i] <- max(tempo$Time)

  # Data corresponding to phase A
  Atemp <- subset(tempo, Phase == "A")
  Aphase <- Atemp$Score
  # Data corresponding to phase B
  Btemp <- subset(tempo, Phase == "B")
  Bphase <- Btemp$Score

  # Obtain the MPD
  results[i] <- MPD(Aphase, Bphase)[[1]]
  # Obtain the weight for the MPD
  tier.weights[i] <-  obs[i]
  # Obtain the standardized-version of the MPD
  stdresults[i] <- MPD(Aphase, Bphase)[[2]]
}  

# Obtain the weighted average: single effect size per study
one_ES_num <- 0

for (i in 1:tiers)
  one_ES_num <- one_ES_num + results[i]*tier.weights[i]
one_ES_den <- sum(tier.weights)
one_ES <- one_ES_num / one_ES_den

# Obtain the weight of the single effect size - it's a weight for the meta-analysis
# Weight 1: series length
weight <- sum(obs)
variation <- 0
if (tiers == 1) weight2 = 1
if (tiers != 1) {
  for (i in 1:tiers)
    variation <- variation + (stdresults[i]-one_stdES)*(stdresults[i]-one_stdES)
  # Weight two: variability of the outcomes
  coefvar <- (sqrt(variation/tiers))/abs(one_stdES)
  weight2 <- 1/coefvar
}
# Weight
weight_std <- weight + weight2

# Obtain the standardized-version of the MPD values for each tier

# Obtain the standardized-version of the weighted average
one_stdES_num <- 0
one_stdES_den <- 0
for (i in 1:tiers)

  if (abs(stdresults[i])!=Inf) 
  { 
     one_stdES_num <-  one_stdES_num + stdresults[i]*tier.weights[i]
     one_stdES_den <- one_stdES_den + tier.weights[i] 
  }

one_stdES <- one_stdES_num / one_stdES_den

data.frame(Id = n, ES = one_ES, stdES = one_stdES, weight = weight_std)

}

Как видите, когда я применяю функцию к одному участнику, она работает без проблем. MPD_meta требует два аргумента, т. е. номер весов и идентификатор участника. При применении дает нужный мне результат.

MPD_meta(1, "C1")
#  Id        ES stdES weight
#1 C1 -13.79167 0.325      7 

Однако, когда я попытался сопоставить эту функцию каждому участнику для одной конкретной шкалы, я не получил желаемого результата. Это дает мне одну и ту же строку, повторяющуюся n раз для каждого участника. Обратите внимание, что результаты по-прежнему верны. Очевидно, я мог бы сказать %>% unique(), и он выбрал бы только один, но это занимает много времени, когда я работаю с большим набором данных.

map2_df(.x = db$Scale, .y = db$ID, ~ MPD_meta(.x, .y))
#   Id         ES     stdES weight
#1  C1 -13.791667 0.3250000      7
#2  C1 -13.791667 0.3250000      7
#3  C1 -13.791667 0.3250000      7
#4  C1 -13.791667 0.3250000      7
#5  C1 -13.791667 0.3250000      7
#6  C1 -13.791667 0.3250000      7
#7  C1 -13.791667 0.3250000      7