CODEWITHSUNDEEP
r
Switters
Switters

Reputation: 119

Analyzing multiple repetitions

I have a dataset with 1000 repetitions. Each repetition contains 50 entries. I need to run and output some statistics for each repetition. It seems like a for loop is the right way to go, but I can't make it work.

Example dataset (truncated to 3 samples per rep) and the code I would like to run on each set of reps

'''R
####DATASET
Rep Y_Star_T m
1 3 3 
1 8 7
1 9 6
2 13 2 
2 5 5
2 19 16
3 12 7 
3  7
3 9 6
####Global variables needed outside the loop
a <- .25
L <- 75
W <- 20
big_N <- (L*W)/a

####Begin appended calculations needed for the loop
w <- netdata$y_star_T/netdata$m
mu_hat <- (1/n *sum(w))/a
tau_hat <- (big_N*a)*mu_hat
var_mu_hat <- (1/(n*(n-1))*sum((w-mu_hat)^2))/a^2
var_tau_hat <- (big_N*a)^2*var_mu_hat

I would like a vector of length equal to the number of repetitions that contain the values calculated above (e.g., mu_hat). For 1000 repetitions, the vector length would be 1000 One vector per for mu_hat, tau_hat, etc is fine if it's easy or a list of lists would work as well.

Upvotes: 0

Views: 83

Answers (1)

Rui Barradas
Rui Barradas

Reputation: 76651

You have to define an auxiliary function to compute the statistics needed. Then *apply that function by groups defined by column Rep.

funStats <- function(DF, a = 0.25, L = 75, W = 20, big_N = L*W/a){
    n <- nrow(DF)
    w <- DF[["Y_Star_T"]]/DF["m"]
    mu_hat <- (1/n *sum(w))/a
    tau_hat <- (big_N*a)*mu_hat
    var_mu_hat <- (1/(n*(n-1))*sum((w-mu_hat)^2))/a^2
    var_tau_hat <- (big_N*a)^2*var_mu_hat
    res <- c(mu_hat = mu_hat,
      tau_hat = tau_hat,
      var_mu_hat = var_mu_hat,
      var_tau_hat = var_tau_hat
    )
    res
}

sapply(split(netdata, netdata$Rep), funStats)
#                       1            2            3
#mu_hat      4.857143e+00 1.158333e+01 6.428571e+00
#tau_hat     7.285714e+03 1.737500e+04 9.642857e+03
#var_mu_hat  1.065170e+02 6.557882e+02 3.721224e+02
#var_tau_hat 2.396633e+08 1.475523e+09 8.372755e+08


result <- lapply(split(netdata, netdata$Rep), funStats)
result <- do.call(rbind, result)

result
#     mu_hat   tau_hat var_mu_hat var_tau_hat
#1  4.857143  7285.714   106.5170   239663265
#2 11.583333 17375.000   655.7882  1475523437
#3  6.428571  9642.857   372.1224   837275510

Data.

Note that row 8 is commented out.

netdata <- read.table(text = "
Rep Y_Star_T m
1 3 3 
1 8 7
1 9 6
2 13 2 
2 5 5
2 19 16
3 12 7 
#3  7
3 9 6
", header = TRUE)

Upvotes: 1

Related Questions