CODEWITHSUNDEEP
rfunctionloopsmatrix
MNU
MNU

Reputation: 764

Multiple function for the same data in r

With the following information,

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)

I can calculate three different functions. The first function is b0 + b1*fu + b2*fu^2+z. and calculated as 

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]
    }
  }
}

Now finding the mean and quantiles as 

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)

The second function is b0 + b1*fu + b2*fu^2+z+2*sqrt(sd.val) and calculated as `

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)

The third function is b0 + b1*fu + b2*fu^2+z-2*sqrt(sd.val) and calculated as 

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)

`

All of the codes work well. My question is, Can I calculate all of these functions in one loop or using any other function?. Any help is appreciated.

Upvotes: 0

Views: 90

Answers (1)

Cole
Cole

Reputation: 11255

To continue on the previous answer:

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

We can simplify your loops to this:

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])

If you notice, the last 3 lines share a lot of elements in common. When we are only adding a constant, we can use sweep. For each column in pt_matrix_no_sd, we are going to add 2 * sd.val with the following statement:

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

Then to get your summary statistics, we can use colMeans or other colXs from matrixStats:

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

Upvotes: 1

Related Questions