R: sum across matrix elements in a vector that are of different dimensions

user_11767393 0
Viewed 60

I have a list containing many matrices that correspond to repeated draws from a count distribution (columns represent the original count data, rows represent the simulated count data, and cell values are counts of observations from one simulation).

I'd like to be able to sum the matrices within the list to form one "super" matrix that contains counts of the observations across the simulations. I know matrices aren't supposed to be added unless they're the same dimension, so if someone has an idea how to use something like grid.expand or so that all the matrices are the same size, I'm open to it.

I should also note I'm not tied to matrices - I could also pull out unequal sized vectors of counts, but I still don't know how to (in a vectorized solution) sum over many unequal sized vectors of counts without recycling. This has probably been answered somewhere but at first glance I couldn't find anything that made sense in a vectorized solution for many simulations.

I made a toy reprex below - mtcars$carb is probably not the ideal variable for this but I wanted a Poisson-distributed variable that is small so just stick with me.

Thanks for any solutions anyone has!

library(tidyverse)

draw <- mtcars %>% 
  select(carb) %>% 
  mutate(carb2 = NA)

# initialize storage vector to hold results from 5 simulations
res <- vector(mode = "list", length = 5)

# run simulation 5 times
for (i in 1:5) { 
  
  # draw new counts once from carb ~ Poisson in carb2 with lambda = mean(carb)
  for (j in 1:nrow(draw)) {
    draw$carb2[j] <- rpois(1, mean(draw$carb))
  }
  
  # get "crosstabs" matrix of obs counts 
  # cols are original values from carb, rows are sim values from carb2; 
  # values are counts
  sim <- draw %>% 
    group_by(carb, carb2) %>% 
    tally() %>% 
    spread(carb, n) %>% 
    replace(is.na(.), 0) %>%   # replace NA with 0
    as.matrix()
  
  res[[i]] <- sim
  
}

Output below:

> res

[[1]]
     carb2 1 2 3 4 6 8
[1,]     0 1 1 0 2 0 0
[2,]     1 1 0 0 3 0 0
[3,]     2 1 4 1 2 0 1
[4,]     3 3 1 0 0 1 0
[5,]     4 0 3 0 2 0 0
[6,]     5 1 1 2 1 0 0

[[2]]
     carb2 1 2 3 4 6 8
[1,]     0 0 2 0 1 0 0
[2,]     1 1 0 0 1 1 0
[3,]     2 3 1 1 1 0 0
[4,]     3 2 2 0 3 0 1
[5,]     4 0 3 1 3 0 0
[6,]     5 0 1 1 1 0 0
[7,]     6 0 1 0 0 0 0
[8,]     7 1 0 0 0 0 0

[[3]]
     carb2 1 2 3 4 6 8
[1,]     1 1 5 0 2 0 0
[2,]     2 2 2 1 3 0 0
[3,]     3 3 2 1 1 0 1
[4,]     4 1 1 0 2 1 0
[5,]     5 0 0 0 1 0 0
[6,]     6 0 0 1 0 0 0
[7,]     8 0 0 0 1 0 0

[[4]]
     carb2 1 2 3 4 6 8
[1,]     0 1 0 0 1 0 0
[2,]     1 0 1 0 1 0 0
[3,]     2 2 2 1 1 0 0
[4,]     3 0 2 0 5 1 1
[5,]     4 1 2 2 0 0 0
[6,]     5 3 3 0 2 0 0

[[5]]
     carb2 1 2 3 4 6 8
[1,]     0 1 0 0 0 0 0
[2,]     1 1 3 1 1 0 0
[3,]     2 2 3 1 2 1 0
[4,]     3 1 0 0 3 0 1
[5,]     4 1 3 0 3 0 0
[6,]     5 1 0 1 0 0 0
[7,]     6 0 1 0 1 0 0

I would like to sum the matrices from the five simulations above into the below format:

carb2   1   2   3   4   6   8
    0   2   1   0   2   0   0
    1   5   11  1   10  1   0
    2   10  12  5   9   1   1
    3   9   7   1   12  2   4
    4   3   12  3   10  1   0
    5   5   5   4   5   0   0
    6   0   2   1   1   0   0
    7   1   0   0   0   0   0
    8   0   0   0   1   0   0
data-wrangling r matrix
1 Answers

Here is a way.

The main idea is to change res's data structure from list to 3d array. The first two dimensions are carb2 and carb, respectively, and the third dimension is the number of simulations. There's no way to guess the upper limit for rpois beforehand so in the code below I choose the 0.999 quantile, hoping that it will be enough. Then assign each iteration's simulation result sim to the array slices with the rows matching carb2. Finally, add all slices.

suppressPackageStartupMessages(
  library(tidyverse)
)
set.seed(2022)

draw <- mtcars %>% 
  select(carb) %>% 
  mutate(carb2 = NA)

cols <- c("carb2", sort(unique(mtcars$carb)))
rows <- 0:qpois(0.999, mean(draw$carb))
R <- 5L
res <- array(NA_integer_, dim = c(length(rows), length(cols), R),
             dimnames = list(rows, cols, 1:R))

for(i in 1:R) {
  draw$carb2 <- rpois(nrow(draw), mean(draw$carb))
  sim <- draw %>% 
    group_by(carb, carb2) %>% 
    tally() %>% 
    spread(carb, n) %>% 
    replace(is.na(.), 0) %>%   # replace NA with 0
    as.matrix()
  j <- match(sim[, "carb2"], rows)
  res[j, , i] <- sim
}

cbind(
  carb2 = seq_len(nrow(res)) - 1L,
  apply(res[, -1, ], 1:2, sum, na.rm = TRUE)
)
#>   carb2  1  2 3  4 6 8
#> 0     0  0  3 0  3 1 0
#> 1     1  6 11 5  6 0 0
#> 2     2  7  9 1  9 1 1
#> 3     3 11 12 4 13 1 2
#> 4     4  5  6 2 10 2 2
#> 5     5  4  6 2  7 0 0
#> 6     6  2  1 0  1 0 0
#> 7     7  0  1 1  1 0 0
#> 8     8  0  0 0  0 0 0
#> 9     9  0  1 0  0 0 0

Created on 2022-07-14 by the reprex package (v2.0.1)

user_82454060
Related