r中相同数据的多个函数

Question

有了以下信息，

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))
z=data.frame(z1=c(0.25,0.47),z2=c(0.48,0.57),z3=c(0.25,0.64))
T.val=data.frame(T1=c(1,1),T2=c(1,2),T3=c(2,1))
dt_data=cbind(b0,b1,b2,T.val,z)
fu.time=seq(0,2,by=0.8)
pat=ncol(T.val) #number of T's
nit=2 #no of rows
sd.val=c(0.48,0.65)

我可以计算三个不同的函数。第一个函数是b0 + b1*fu + b2*fu^2+z.，计算为

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+z[it.er,pt]
    }
  }
}

现在找到平均值和分位数

pt.array1.mean=apply(pt.array1,c(3,2), mean)
pt.array1.LCI=apply(pt.array1,c(3,2), quantile, prob=0.25)
pt.array1.UCI=apply(pt.array1,c(3,2), quantile, prob=0.975)

第二个函数是b0 + b1*fu + b2*fu^2+z+2*sqrt(sd.val)，计算如下 `

pt.array_UPI=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.array_UPI[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+z[it.er,pt]+2*sqrt(sd.val[it.er])
    }
  }
}

pt.array_UPI.mean=apply(pt.array_UPI, c(3,2), mean)

第三个功能是 b0 + b1*fu + b2*fu^2+z-2*sqrt(sd.val) 并计算为

pt.array_LPI=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.array_LPI[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+z[it.er,pt]+2*sqrt(sd.val[it.er])
    }
  }
}

pt.array_LPI.mean=apply(pt.array_LPI, c(3,2), mean)

`

所有代码都运行良好。我的问题是， 我可以在一个循环中或使用任何其他函数计算所有这些函数吗？ 任何帮助表示赞赏。

Answer 1

继续上一个答案：

storage problem in R. alternative to nested loop for creating array of matrices and then multiple plots

我们可以将您的循环简化为：

library(matrixStats)

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], z[mat_ind], sd.val[ind$nits]), ncol = 5)

colnames(b_mat) <- c('b0','b1','b2','z','sd.val')
b_mat 

#         b0       b1         b2    z sd.val
#[1,] 11.410  0.85390 -0.0130200 0.25   0.48
#[2,] 11.360  0.95650 -0.0165400 0.47   0.65
#[3,] 11.410  0.85390 -0.0130200 0.25   0.48
#[4,]  6.950  0.06752 -0.0026720 0.57   0.65
#[5,]  8.767 -0.03179 -0.0002822 0.48   0.48
#[6,] 11.360  0.95650 -0.0165400 0.47   0.65

pt_matrix_no_sd <- apply(b_mat, 1, function(x) x[1] + x[2] * fu.time + x[3] * fu.time^2 + x[4])
pt_matrix_pos_sd <- apply(b_mat, 1, function(x) x[1] + x[2] * fu.time + x[3] * fu.time^2 + x[4] + 2 * x[5])
pt_matrix_neg_sd <- apply(b_mat, 1, function(x) x[1] + x[2] * fu.time + x[3] * fu.time^2 + x[4] - 2 * x[5])

如果您注意到，最后 3 行有很多共同点。当我们只添加一个常量时，我​​们可以使用sweep。对于pt_matrix_no_sd 中的每一列，我们将使用以下语句添加2 * sd.val：

sweep(pt_matrix_no_sd, 2, 2*b_mat[,'sd.val'], FUN = '+')

identical(pt_matrix_pos_sd,
          sweep(pt_matrix_no_sd, 2, 2*b_mat[,'sd.val'], FUN = '+')
) #TRUE

然后要获取您的汇总统计信息，我们可以使用colMeans 或matrixStats 中的其他colXs：

library(matrixStats)
pt_summary = array(t(apply(pt_matrix_no_sd, #change as needed
                         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, ,]   

#        mean      LCL      UCL
#pat1 11.7450 11.70250 11.82575
#pat2  9.5900  8.55500 11.55650
#pat3 10.5385  9.89275 11.76543

pt_summary2 = array(t(apply(pt_matrix_pos_sd, 1, #change as needed
                           function(row) colMeans2(matrix(row, ncol = pat)))),
                    dim = c(length(fu.time), pat, 1),
                    dimnames = list(NULL, paste0('pat', seq_len(pat)), c('mean')))
pt_summary2[1,,]

#   pat1    pat2    pat3 
#12.8750 10.7200 11.6685 

#you should be able to do the negative sd