# Sampling and resampling data in R

My problem this time concerns sampling size-related errors, resample-based confidence intervals and a possible way to control for this error. My dataset consists of 50 measurements of certain cranial dimensions, from several groups (populations). My aim is to calculate a certain index based on a correlation matrix calculated for each of the populations separately. Now the problem is that sample sizes vary from population to population, and they are smaller in some, around 20 and large, around 100 in others. Some reviewers suggested that the sample size inequality surely leads to sampling error, such that as calculated index gets larger, sampling error also gets larger while the sample size is smaller. This is bad for me, as my work is based on finding the differences in this index.

This problem was also the subject of one of my earlier posts, and it shows hof far have I got with this dataset, and I think is more suitable for Cross Validated.

I will try to provide all the code and a possible reproducible example, although correlations between measurements are higher than in generated data:

First I generate some data:

rex <- rnorm(25000, mean = 100, sd = 25)
testMat <- matrix(rex, nrow = 500, ncol = 50)
populations <- factor(c(rep("a",25), rep("b",45), rep("c",100), rep("d",50), rep("e",30), rep("f", 50), rep("g", 50), rep("h",35), rep("i",45), rep("a",70)))
testMat <- data.frame(populations, testMat)


This is the function I use for the said index (standardized variance of eigenvalues)

IIvar <- function(R) {
R <- cor(R)
d <- eigen(R)$values p <- length(d) sum((d-1)^2)/(p*(p-1)) }  And I need to provide some measure of uncertainty so I wrote this function that is based on resampling from the whole matrix and calculating the index for every population (a, b, c, ...) separately. omat is the separate matrix for just one population, mat1 is the whole data matrix (testMat), freq is the number of individuals in a population, and numR is the number of resamples. This function gives the index and two types of confidence intervals. These intervals should be ok since I tested the same data with the boot package as well. ciint <- function(omat, mat1, freq, numR) { II <- IIvar(omat) n <- dim(mat1)[1] b <- numeric(numR) for (i in 1:numR) { b[i] <- IIvar(cor(mat1[sample(c(1:n),freq, replace = TRUE),]))} cinfmin <- II + qnorm(0.025)*sd(b) #normal intervals cinfmax <- II + qnorm(0.975)*sd(b) cinfminT <- (sqrt(II) + qnorm(0.025)*sd(sqrt(b)))^2 #transformed normal intervals cinfmaxT <- (sqrt(II) + qnorm(0.975)*sd(sqrt(b)))^2 list(integ = II, c25 = cinfmin, c95 = cinfmax, c25t = cinfminT, c95t = cinfmaxT) }  Unfortunately, when I use these functions I encounter a problem, like smaller the population sample, smaller the index and larger the confidence interval. I wanted to check how sampling size affects these two values so I gone further with this, by writing a random sampling function that will take a subset of individuals from every population, possibly with smaller number of measurements (sample rows and columns as well). These two functions perform this task, and the matrix I use is testMat, which must have a grouping variable in its first column. sampleria <- function(data, groups = 10, num = 20, num2 = 51, sampleCol = FALSE) { pops <- sample(unique(data[,1]),groups, replace = FALSE) pops <- sort(pops) b <- subset(data, data[,1] == pops[1]) b <- b[sample(1:nrow(b), num, replace = FALSE),] for(i in 2:groups) { a <- subset(data, data[,1] == pops[i]) a <- a[sample(1:nrow(a), num, replace = FALSE),] b <- rbind(b,a) } if (sampleCol == TRUE) { nd <- b[,2:ncol(b)] nd <- nd[,sample(1:ncol(nd), num2, replace = FALSE)] nd <- cbind(b[,1],nd) names(nd)[1] <- "groups" return(nd) } else { return(b) } }  sampleria returns the new dataset, with i.e. 20 individuals from every population so that ciintSampler <- function(mat1, numR, numy = 20, num2 = 51, sampleCol = FALSE) { matSample <- sampleria(mat1, groups = 10, num = numy, num2 = num2, sampleCol = sampleCol) pops <- unique(matSample[,1]) integration <- numeric(length(pops)) integration.ci <- numeric(length(pops)) for (i in 1:length(pops)) { omat <- subset(matSample, matSample[,1] == pops[i]) a <- ciint(omat[,2:num2], matSample[,2:num2], numy, numR) integration[i] <- a$integ
integration.ci[i] <- a$c95t - a$c25t
}
integrationMat <- data.frame(pops, integration, integration.ci)
return(integrationMat)
}


finally returns the index and its confidence intervals which could be compared. When I repeat this procedure, i.e., 10 times, changing the number of sampled individuals or the number of characters, I expected that the relationship between the index value and its CI would not be strict, linear and definite, with r2 value around 0.9 almost every time. I must be doing somethinga wrong since it is known that if sample size is less than the number of measurements the uncertainty in estimates rises.

I know this is a long post, but it is very important for me and I've been struggling with this a long time, so any idea would be greatly appreciated.

• Let's focus on the most fundamental issue: why have you chosen this formula for your index? What feature of the population is it intended to represent? How do you intend to interpret it? – whuber Jan 4 '14 at 20:11
• This is the index of morphological integration that represents the amount of connectedness between the measured cranial dimensions. It is calculated as standardized variance of correlation matrix eigenvalues, and if it is larger it means that the cranial elements (represented by measurements) are more tightly integrated. My hypothesis is that populations (biological) should have different degree of integration as a consequence of environmental differences. – Fedja Blagojevic Jan 5 '14 at 18:49
• I am confused because that is not what your IIvar function appears to be computing. That variance of the eigenvalues would correctly be computed using a denominator of p rather than p*(p-1). However, that's unimportant because this is just a constant factor applied regardless of the sample size. Are you asking how to estimate the standard error of this index? If so, I would expect the SE to be large, because the index ought to be unstable when sample sizes are smaller than or on the order of the number of variables, as is the case here. – whuber Jan 5 '14 at 22:21
• Yes yes, this is the exact issue, the inequality of sample sizes. But the goal of all these functions is to try to sample equal portions of individuals from every population (i.e. 20 is the maximum available), calculate the index again and see if it changes. Additionally, I would like to provide some CIs around the observed value of IIvar, but they are keep getting bigger (wider) as the observed value is higher. Finally, I am unsure how a standard error would be calculated for this index, and that is why I wanted to use bootstrap confidence intervals. – Fedja Blagojevic Jan 6 '14 at 19:15