title: Mapping a function to each participant of a study
body: <p>I have been trying to apply a function to this data set that consists of data from single studies AB designed. It consists of 6 variables: a. Tier (number of AB cycle, in this particular data set is always 1 because all studies are AB designed); b. ID (participant ID code), c. Scale (type of scale on which participant is evaluated), d. Time (occasion of the rating), e. Phase (A = baseline, B = intervention), f. Score (participant's score).</p>
  
  <pre><code>library(purrr)
  library(dplyr)
  
  Tier &lt;- rep(c(1), 36)
  ID &lt;- 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 &lt;- rep(c(1, 2), each = 18)
  Time &lt;- 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 &lt;- 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 &lt;- 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 &lt;- data.frame(Scale, ID, Time, Phase, Score, Tier)
  </code></pre>
  
  <p>These are two functions, developed by Manolov and Rochat (2015). They work perfectly until it comes to calculate more than one participant at one time. So I tried to change some details of the second one (MPD_meta) to do so. </p>
  
  <pre><code>MPD &lt;- function(a, b) { #Manolov &amp; Rochat, 2015
  
  # Obtain phase length
  n_a &lt;- length(a[!is.na(a)])
  n_b &lt;- length(b[!is.na(b)])
  meanA &lt;- mean(a[!is.na(a)]) 
  
  # Estimate baseline trend
  base_diff &lt;- rep(0,(n_a - 1))
  
  for (i in 1:(n_a - 1)) base_diff[i] &lt;- a[!is.na(a)][i + 1] - a[!is.na(a)][i]
  
  trendA &lt;- mean(base_diff)
  
  # Predict baseline data
  baseline_pred &lt;- rep(0, n_a)
  
  midX &lt;- median(c(1:n_a))
  midY &lt;- median(a[!is.na(a)])
  
  is.wholenumber &lt;- function(x, tol = .Machine$double.eps ^ 0.5)  abs(x - round(x)) &lt; tol
  
  if (is.wholenumber(midX)) {
  baseline_pred[midX] &lt;- midY
  for (i in (midX - 1): 1) baseline_pred[i] &lt;- baseline_pred[i + 1] - trendA
  for (i in (midX + 1): n_a) baseline_pred[i] &lt;- baseline_pred[i - 1] + trendA
  }
  
  #baseline_pred
  
  if (!is.wholenumber(midX)) {
  baseline_pred[midX - 0.5] &lt;- midY - (0.5 * trendA)
  baseline_pred[midX + 0.5] &lt;- midY + (0.5 * trendA)
  for (i in (midX - 1.5) : 1) baseline_pred[i] &lt;- baseline_pred[i + 1] - trendA
  for (i in (midX + 1.5) : n_a) baseline_pred[i] &lt;- baseline_pred[i - 1] + trendA
  }
  
  
  # Project baseline trend to the treatment phase
  treatment_pred &lt;- rep(0, n_b)
  treatment_pred[1] &lt;- baseline_pred[n_a] + trendA
  
  for (i in 2:n_b) 
      treatment_pred[i] &lt;- treatment_pred[i - 1] + trendA
  
  diff_absolute &lt;- rep(0, n_b)
  diff_relative &lt;- rep(0, n_b)
  
  for (i in 1:n_b) {
  diff_absolute[i] &lt;- b[!is.na(b)][i] - treatment_pred[i]
  
  if (treatment_pred[i]!=0) diff_relative[i] &lt;- ((b[!is.na(b)][i] - treatment_pred[i]) / abs(treatment_pred[i])) * 100
      }
  
  mpd_value &lt;- mean(diff_absolute)
  mpd_relative &lt;- mean(diff_relative)
  mpd_relsd &lt;- mpd_value/sd(a[!is.na(a)])
  
  # Compute the residual (MSE): difference between actual and predicted baseline data
  residual &lt;- 0
  
  for (i in 1:n_a)
  
    residual &lt;- residual + (a[!is.na(a)][i] - baseline_pred[i]) * (a[!is.na(a)][i] - baseline_pred[i])
    residual_mse &lt;- 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 &lt;- function(DBR, n) {
  
  MBD &lt;- db %&gt;%
    dplyr::filter(Scale == DBR, ID == n)
  
  # Create the objects necessary for the calculations and the graphs
  # Dimensions of the data set
  tiers &lt;- max(MBD$Tier)
  max.length &lt;- max(MBD$Time)
  max.score &lt;- max(MBD$Score)
  min.score &lt;- min(MBD$Score)
  
  # Vectors for keeping the main results
  # Series lengths per tier
  obs &lt;- rep(0,tiers)
  # Raw MPD effect sizes per tier
  results &lt;- rep(0,tiers)
  # Weights per tier
  tier.weights &lt;- rep(0,tiers)
  # Standardized MPD results per tier
  stdresults &lt;- rep(0,tiers)
  
  # Obtain the values calling the MPD function
  for (i in 1:tiers) {
  
    tempo &lt;- subset(MBD, Tier == i)
    # Number of observations for the tier
    obs[i] &lt;- max(tempo$Time)
  
    # Data corresponding to phase A
    Atemp &lt;- subset(tempo, Phase == "A")
    Aphase &lt;- Atemp$Score
    # Data corresponding to phase B
    Btemp &lt;- subset(tempo, Phase == "B")
    Bphase &lt;- Btemp$Score
  
    # Obtain the MPD
    results[i] &lt;- MPD(Aphase, Bphase)[[1]]
    # Obtain the weight for the MPD
    tier.weights[i] &lt;-  obs[i]
    # Obtain the standardized-version of the MPD
    stdresults[i] &lt;- MPD(Aphase, Bphase)[[2]]
  }  
  
  # Obtain the weighted average: single effect size per study
  one_ES_num &lt;- 0
  
  for (i in 1:tiers)
    one_ES_num &lt;- one_ES_num + results[i]*tier.weights[i]
  one_ES_den &lt;- sum(tier.weights)
  one_ES &lt;- 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 &lt;- sum(obs)
  variation &lt;- 0
  if (tiers == 1) weight2 = 1
  if (tiers != 1) {
    for (i in 1:tiers)
      variation &lt;- variation + (stdresults[i]-one_stdES)*(stdresults[i]-one_stdES)
    # Weight two: variability of the outcomes
    coefvar &lt;- (sqrt(variation/tiers))/abs(one_stdES)
    weight2 &lt;- 1/coefvar
  }
  # Weight
  weight_std &lt;- weight + weight2
  
  # Obtain the standardized-version of the MPD values for each tier
  
  # Obtain the standardized-version of the weighted average
  one_stdES_num &lt;- 0
  one_stdES_den &lt;- 0
  for (i in 1:tiers)
  
    if (abs(stdresults[i])!=Inf) 
    { 
       one_stdES_num &lt;-  one_stdES_num + stdresults[i]*tier.weights[i]
       one_stdES_den &lt;- one_stdES_den + tier.weights[i] 
    }
  
  one_stdES &lt;- one_stdES_num / one_stdES_den
  
  data.frame(Id = n, ES = one_ES, stdES = one_stdES, weight = weight_std)
  
  }
  </code></pre>
  
  <p>As you can see, when I apply the function to one single participant, it works with no problem. <code>MPD_meta</code> requires two arguments, i.e. the number of the scale and the participant ID. When applied, it gives the result I want.</p>
  
  <pre><code>MPD_meta(1, "C1")
  #  Id        ES stdES weight
  #1 C1 -13.79167 0.325      7 
  </code></pre>
  
  <p>However, when I tried to map this function to each participant for one specific scale, I failed to obtain the result I wanted. It gives me the same line repeated n times for each participant. Note the results are still correct. Obviously, I could say <code>%&gt;% unique()</code> and it would select only one but it takes lots of time when I deal with a larger data set.</p>
  
  <pre><code>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
  </code></pre>