1
$\begingroup$

I'm doing a simulation. I first simulate a mean vector and a semi - positive definite covariance matrix.

I then go through and generate 5000 samples of 10 length. I then construct the 95% Confidence interval for each sample. Finally, I see how many samples overlap the true means (which I know). However it's all of them. Shouldn't there be a few that fall out of the range?

The code is attached.

library(MASS)
library(clusterGeneration)

#set.seed(314159) # For ease of replication

######################
#Generate mean vector and sigma information
######################
trueMu<- runif(n = 2, min = 0, max = 2)
matrixData <- genPositiveDefMat(dim = 2, covMethod = "eigen", lambdaLow = 8)
trueSigma <- matrixData$Sigma
trueMu
trueSigma

######################
#Generate random data 5000 times @ sizes of 10,20, 150
######################
TenSampleList <- lapply(1:5000,function(x){
  mvrnorm(n = 10, mu = trueMu, Sigma = trueSigma)  
})

######################
#Function to calculate confidence interval
######################
calcConfInterval <- function(currentData, corrected= FALSE){
  #browser()
  meanCol1 <- mean(currentData[,1])
  meanCol2 <- mean(currentData[,2])

  sdCol1  <- sqrt(var(currentData[,1]))
  sdCol2  <- sqrt(var(currentData[,2]))

  lengthObs <- nrow(currentData)
  confidenceLevel = 1 - .05/(2 * ifelse(corrected ==TRUE, 2 ,1))

  col1Conf <- c(meanCol1 - qt(confidenceLevel, lengthObs - 1)*sdCol1,meanCol1 + qt(confidenceLevel, lengthObs - 1)*sdCol1)
  col2Conf <- c(meanCol2 - qt(confidenceLevel, lengthObs - 1)*sdCol2,meanCol2 + qt(confidenceLevel, lengthObs - 1)*sdCol2)

  return(list(col1Conf = col1Conf, col2Conf = col2Conf))
}

######################
#Get ConfidenceIntervals for each sampel set
######################

TenSampleConfIntUncorrected <- lapply(TenSampleList, calcConfInterval, corrected = FALSE)

checkIfInside <- function(data, trueMu1, trueMu2){
  #browser()
  result1 <- if(trueMu1 >= data$col1Conf[1] && trueMu1 <= data$col1Conf[2]){1}else{0}
  result2 <- if(trueMu2 >= data$col2Conf[1] && trueMu2 <= data$col2Conf[2]){1}else{0}
  return(result1 * result2)
}


TenSampleUncorrectedInsideProportion <- sum(unlist(lapply(TenSampleConfIntUncorrected, checkIfInside, trueMu1 = trueMu[1], trueMu2 = trueMu[2])))

TenSampleUncorrectedInsideProportion
$\endgroup$
  • 2
    $\begingroup$ The standard error on the mean is the one for the population divided by $\sqrt{n}$. Moreover you take 0.95 confidence on each dimension so you have a 0 9025 CI ( if I understand your code well) $\endgroup$ – user83346 Nov 14 '16 at 5:31
  • 1
    $\begingroup$ Ack! Such a simple thing. Thanks for the reminder. And yes, I realize that it's 0.9025 CI for both together. That's why I have the ifelse in the "confidenceLevel" line. To show that after I do bonferroni's correction it's more on target to 0.95 $\endgroup$ – user1357015 Nov 14 '16 at 5:49
  • 1
    $\begingroup$ Fine :-) so if you use the $\sqrt{n}$ then you should have some non-covering intervals? $\endgroup$ – user83346 Nov 14 '16 at 6:17
  • 1
    $\begingroup$ Yep, I do. Using the sqrt(n) it's about 90.something % and after I use the correction it's about 95%-- exactly as I expected. The mistake to forget sqrt(n) was silly. Thank you for pointing it out. Much appreciated. If you put it as an answer I'll mark it correct and give you credit. $\endgroup$ – user1357015 Nov 14 '16 at 6:20
  • 1
    $\begingroup$ Your title gives misleading impression of the question. Maybe you could change it to: "Why my simulation does not result in correct coverage for confidence interval?" or something similar. $\endgroup$ – Tim Nov 14 '16 at 8:18
4
$\begingroup$

As indicated in my comment, the standard error on the mean should be the one of the random variable divided by $\sqrt{n}$.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.