R中的存储问题。替代用于创建矩阵数组和多个绘图的嵌套循环

Posted 2023-03-16

【中文标题】R中的存储问题。替代用于创建矩阵数组和多个绘图的嵌套循环

【英文标题】：storage problem in R. alternative to nested loop for creating array of matrices and then multiple plots

【发布时间】：2019-12-31 22:29:22

【问题描述】：

通过以下信息，我可以轻松创建矩阵数组

b0=data.frame(b0_1=c(11.41,11.36),b0_2=c(8.767,6.950))
b1=data.frame(b1_1=c(0.8539,0.9565),b1_2=c(-0.03179,0.06752))
b2=data.frame(b2_1=c(-0.013020 ,-0.016540),b2_2=c(-0.0002822,-0.0026720))
T.val=data.frame(T1=c(1,1),T2=c(1,2),T3=c(2,1))
dt_data=cbind(b0,b1,b2,T.val)
fu.time=seq(0,50,by=0.8)
pat=ncol(T.val) #number of T's
nit=2 #no of rows

pt.array1=array(NA, dim=c(nit,length(fu.time),pat)) 

for ( it.er in 1:nit)
  for ( ti in 1:length(fu.time))
    for (pt in 1:pat)
      pt.array1[it.er,ti,pt]=b0[it.er,T.val[it.er,pt]]+b1[it.er,T.val[it.er,pt]]*fu.time[ti]+b2[it.er,T.val[it.er,pt]]*fu.time[ti]^2
    
  


pt.array_mean=apply(pt.array1, c(3,2), mean)
pt.array_LCL=apply(pt.array1, c(3,2), quantile, prob=0.25)
pt.array_UCL=apply(pt.array1, c(3,2), quantile, prob=0.975)

现在有了这些额外的数据，我可以创建如下三个图

    mydata
       pt.ID      time IPSS
1      1  0.000000   10
2      1  1.117808    8
3      1  4.504110    5
4      1  6.410959   14
5      1 13.808220   10
6      1 19.890410    4
7      1 28.865750   15
8      1 35.112330    7
9      2  0.000000    6
10     2  1.117808    7
11     2  4.109589    8
12     2 10.093151    7
13     2 16.273973   11
14     2 18.345205   18
15     2 21.567120   14
16     2 25.808220   12
17     2 56.087670    5
18     3  0.000000    8
19     3  1.413699    3
20     3  4.405479    3
21     3 10.389041    8


pdf("plots.pdf")
par(mfrow=c(3,2))
for( pt.no in 1:pat)
  plot(IPSS[ID==pt.no]~time[ID==pt.no],xlim=c(0,57),ylim=c(0,35),type="l",col="black",
      xlab="f/u time", ylab= "",main = paste("patient", pt.no),data=mydata)
  points(IPSS[ID==pt.no]~time[ID==pt.no],data=mydata)
  lines(pt.array_mean[pt.no,]~fu.time, col="blue")
  lines(pt.array_LCL[pt.no,]~fu.time, col="green")
  lines(pt.array_UCL[pt.no,]~fu.time, col="green")

dev.off()

当每个矩阵中的行数比 10000 大得多时，就会出现问题。为 b0、b1 和 b2 中的大量行创建 pt.array1 需要花费太多计算时间。 有没有其他方法可以使用任何内置函数快速完成？ 我可以避免为pt.array1 分配存储空间，因为我不再使用它了吗？对于myplot，我只需要pt.array_mean、pt.array_UCL 和pt.array_LCL。 任何帮助表示赞赏。

【问题讨论】：

@Rui Barradas 你能检查一下吗

你想在这里实现什么？如果你能指导我一点，它可能比嵌套循环快得多

for 循环耗时过长？我可以使用任何其他函数来创建pt.array1 矩阵吗？我又发了一个类似的帖子。 ***.com/questions/57618368/… 。你可以看看

我很乐意帮助您解决您对该帖子的要求，但我对数组没有经验，抱歉

【参考方案1】：

您还可以采用其他几种方法。

首先，您的模型主要是b0 + b1*fu + b2*fu^2。因此，您可以制作系数并在事后应用fu：

ind <- expand.grid(nits = seq_len(nit), pats = seq_len(pat))
mat_ind <- cbind(ind[, 'nits'], T.val[as.matrix(ind)])

b_mat <- matrix(c(b0[mat_ind], b1[mat_ind], b2[mat_ind]), ncol = 3)

b_mat
       [,1]     [,2]       [,3]
[1,] 11.410  0.85390 -0.0130200
[2,] 11.360  0.95650 -0.0165400
[3,] 11.410  0.85390 -0.0130200
[4,]  6.950  0.06752 -0.0026720
[5,]  8.767 -0.03179 -0.0002822
[6,] 11.360  0.95650 -0.0165400

现在，如果我们将模型应用于每一行，我们将获得您的所有原始结果。唯一的问题是我们与您的原始输出不匹配 - 您数组的每个列切片都相当于我的矩阵输出的一个行切片。

pt_array <- apply(b_mat, 1, function(x) x[1] + x[2] * fu.time + x[3] * fu.time^2)

pt_array[1,]
[1] 11.410 11.360 11.410  6.950  8.767 11.360

pt.array1[, 1, ]
      [,1]  [,2]   [,3]
[1,] 11.41 11.41  8.767
[2,] 11.36  6.95 11.360

没关系，因为我们可以在获得汇总统计信息时修复它的形状 - 我们只需要将每行的 colSums 和 colQuantiles 转换为 2 x 3 矩阵：

library(matrixStats)

pt_summary = array(t(apply(pt_array,
                         1,
                         function(row) 
                           M <- matrix(row, ncol = pat)
                           c(colMeans2(M),colQuantiles(M, probs = c(0.25, 0.975))
                           )
                           
                         )),
                   dim = c(length(fu.time), pat, 3),
                   dimnames = list(NULL, paste0('pat', seq_len(pat)), c('mean', 'LCL', 'UCL'))
)

pt_summary[1, ,] #slice at time = 1

        mean      LCL      UCL
pat1 11.3850 11.37250 11.40875
pat2  9.1800  8.06500 11.29850
pat3 10.0635  9.41525 11.29518

# rm(pt.array1)

然后为了进行最终的绘图，我对其进行了简化 - data 参数可以是 subset(mydata, pt.ID == pt.no)。此外，由于汇总统计信息现在采用数组格式，matlines 允许一次完成所有操作：

par(mfrow=c(3,2))

for( pt.no in 1:pat)
  plot(IPSS~pt.ID, data=subset(mydata, pt.ID == pt.no),
       xlim=c(0,57), ylim=c(0,35),
       type="l",col="black", xlab="f/u time", ylab= "",
       main = paste("patient", pt.no)
       )

  points(IPSS~time, data=subset(mydata, pt.ID == pt.no))

  matlines(y = pt_summary[,pt.no ,], x = fu.time, col=c("blue", 'green', 'green'))

【讨论】：

与sd=c(0.5879,0.8946) 我想再添加两个函数b0 + b1*fu + b2*fu^2+2*sqrt(sd) 和b0 + b1*fu + b2*fu^2-2*sqrt(sd) 然后只取mean 值而不是LCL 和UCL。如何添加这两个功能？感谢您的帮助@Cole

您可能应该提出另一个问题。但如果这只是一个标量，看起来你可以做pt_array_sd1_pos<- pt_array +2*sqrt(sd[1])。会是劳动密集型的，但工作