How to assess correlation when each variable is measured by independent replicates? I frequently measure multiple variables, in multiple replicates, at many sites. For instance, I might measure bacterial abundance and bacterial growth rates, each in 3 replicates, at many sites. Each replicate is independently sampled and each variable is measured in a different sample (i.e., I can't measure both bacterial abundance and growth rate in the same sample).
I'd like to test for correlation between those variables. The problem is, since the variables are measured independently, the variables are not paired. Replicate 1 of variable A is not related to replicate 1 of variable B, any more than it is to replicate 2 of variable B. 
I could test for correlation among the average of replicates at each site - but that seems blunt, since you lose information about variation among replicates for each parameter. I can imagine some kind of resampling approach, where I randomly select one replicate for each variable at each site. Is there a better way?
 A: From your description I think the only viable way to go is what you don't want to do: Use the bucket as the level of analyses. That is, aggregate the 3 measurements within each bucket and you have your pairings. With this approach you should effectively aggegate out the measurement error.
I made a small simulation and compared this approach with a second approach in which I used all possible pairings per bucket to estimate the correlation. The results show, the aggregation method is better in recovering the original correlation:
# I use R and the mvtnorm library to generate the data
library(mvtnorm)

set.seed(12345) # make reproducible

nbuckets <- 50  #number of buckets
r.buckets <- 0.5  # correlation across buckets

# generate data
Cor <- array(c(1, r.buckets, r.buckets, 1), dim=c(2,2))
d.bucket <- rmvnorm(nbuckets, sigma = Cor)
measurement.error = 0.5 # size of eror in relation to sd of the data
data <- vector("list", nbuckets)

for (bucket in seq_len(nbuckets)) {
    data[[bucket]] <- list(x = rep(d.bucket[bucket, 1], 3) + rnorm(3, measurement.error), y = rep(d.bucket[bucket, 2], 3) + rnorm(3, sd = measurement.error))
}
# Note that there are separate error terms for the two types of measurements 

# aggregating per bucket:
data.agg <- lapply(data, function(x) data.frame(x = mean(x[[1]]), y = mean(x[[2]])))
data.agg <- do.call("rbind", data.agg)
cor(data.agg$x, data.agg$y) # should give .408

# using all pairs:
all.pairs <- lapply(data, function(x) data.frame(x = x[[1]], y = x[[2]][c(1:3,3:1,2,1,3,2,3,1,1,3,2,3,1,2)]))
all.pairs <- do.call("rbind", all.pairs)
cor(all.pairs$x, all.pairs$y) # should give .321

If you allow for even larger measurement error (although it already is large) the difference remains. If you allow for a single error term within each bucket the results will be a lot nearer to the real value of r and the difference between the methods will decrease. However, aggregating remains the better tactic.
I recommend you play around with it a little with more realistic values. You may even implement a bootstrap approach as was your initial thought. 
