7
$\begingroup$

Let $p$ be a matrix, where each row represents an observation of a 4-variate normally distributed random variable $\mathcal{N}_4(\mu,\Sigma)$.

  • Is there any Bootstrap methode to get a good estimation for Σ?

  • If not, is the number of the following sample enough to bootstrap the distribution of any statistic T which operates on the population where Y comes from?

Here's what I tried so far:

Y<-data.frame(response=c(10,19,27,28,9,13,25,29,4,10,20,18,5,6,12,17),
               treatment=factor(rep(1:4,4)),
               subject=factor(rep(1:4,each=4))
               )

p<-matrix(Y$response,4,4,byrow=T)
B<-1000
sampleB<-sample(1:4,4*B,replace=T)
fit<-lm(p[sampleB,]~1)
cov(residuals(fit))

I also tried

require(nlme)
require(mgcv)
nSubj <- 20
sampleB<-sample(1:4,nSubj,replace=T)
y<-data.frame(response=c(t(p[sampleB,])),
           treatment=factor(rep(1:4,nSubj)),
           subject=factor(rep(1:nSubj,each=4))
           )

fit <- lme(response~-1+treatment,y,random=~1|subject,correlation=corSymm())
extract.lme.cov(fit,y)[1:4,1:4]

but I get the error code:

Error in lme.formula(response ~ -1 + treatment, y, random = ~1 | subject,  : 
nlminb problem, convergence error code = 1
message = iteration limit reached without convergence (10)
$\endgroup$
5
  • $\begingroup$ :there is a boot package in R. You can see that. $\endgroup$ – Metrics Jul 15 '13 at 12:19
  • $\begingroup$ I am familiar with this package. But thats not the point of my question. For sure I could write a function to do what I implement in the code above and use the boot-function. What I am interested in is, if the bootstrap make sense in this case and when this bootstrap method is meaning full how many replications I need. $\endgroup$ – Klaus Jul 15 '13 at 12:31
  • 3
    $\begingroup$ I noticed a flag for migration to SO, but I think there really is a statistical question beyond R code. $\endgroup$ – chl Jul 15 '13 at 13:37
  • 3
    $\begingroup$ You are asking about a proposed solution to a problem. Let's back up: the problem appears to be to estimate the covariance matrix. The R command cov performs that estimate and it is known to be unbiased. Bootstrapping can accomplish two things: (1) it can assess the amount of bias (and let you correct for it) and (2) when the dataset is sufficiently large, it can display the samling distribution. You don't need (1) because your estimates are unbiased and (2) is out of the question. Why then look any further at bootstrapping, unless it's purely a programming exercise? $\endgroup$ – whuber Jul 16 '13 at 13:57
  • $\begingroup$ Thx for this answer, it approves my fears. 1) In this case of small sample size cov gives us not a good inference of $\Sigma$. You will see that cov(p) is not positive semidefinite, so it gives us either a singular or non singular estimation for a normal distributed population. And if I understand you right the (2) point is obsolet because the sample is to small. So what to do with data like this in the context of a mixed model? In my opinion this data is useless to do some statistic inference on it. $\endgroup$ – Klaus Jul 16 '13 at 14:12
5
$\begingroup$

Is the Bootstrap a good option to estimate $\Sigma$?

No, the boostrap will help you infer about the uncertainty of your sample estimate. Specifically, it might be used to get confidence intervals on the elements of $\widehat{\Sigma}$.

And how does the bootstrap work?

The approach is to create $R$ replicate datasets from the original dataset by resampling the observations with replacement. Then you compute the estimate of interest on each of the $R$ replicates, in your case the covariance matrix, for each of the $R$ replicates, obtaining $\widehat{\Sigma}^1, \ldots, \widehat{\Sigma}^R$. Confidence intervals for $\widehat{\Sigma}$ can then be computed empiricaly from $\widehat{\Sigma}^1, \ldots, \widehat{\Sigma}^R$.

For further information, you might want to have a look at the Wikipedia page.


Edit

Could I use $\widetilde{\Sigma}_R = R^{-1} \sum_{r=1}^R \widehat{\Sigma}^r$ instead of $\widehat{\Sigma}$?

Actually, the matrix $\widetilde{\Sigma}_R$ is an estimate of ${\rm E} (\widehat{\Sigma}^r)$, where $\widehat{\Sigma}^r$ is the estimate of the covariance matrix based on a bootstrap replication of the initial sample. Bootstrapping comes down to sample from the empirical distribution $\widehat{F}$. Therefore, I think, but I don't have a formal proof, that $\widetilde{\Sigma}_R = {\rm E} (\widehat{\Sigma}^r) \to \widehat{\Sigma}$ as $R \to \infty$.

So, I think you could use $\widetilde{\Sigma}_R$ instead of $\widehat{\Sigma}$, but that would be like using a sledgehammer to crack a nut.

The R code below is a numerical investigation, which is by no means of proof of the above assertion about the convergence. The structure of dependence is a Gumbel copula, and the margins are two standard normal distribution.

## Initialization
library(copula)
set.seed(531)
n <- 200            # Number of observations in the original sample
R <- 10000          # Number of replications
## Specification for the dependence structure (Gumbel copula)
spec.cl <- archmCopula("gumbel", 1.2)
## Create a fake original dataset
pseudo  <- rCopula(n, spec.cl)
obs     <- qnorm(pseudo)
cov.obs <- cov(obs)[1, 2]
## Get an idea of the "true" covariance
pseudo   <- rCopula(10000, spec.cl)
obs.big  <- qnorm(pseudo)
cov.true <- cov(obs.big)[1, 2]
## Get the bootstrap covariances
cov.sim <- sapply(1:R,
                  function(i, x, n){x.boot <- x[sample(1:n, size = n, replace = TRUE), ]
                                    cov(x.boot)[1, 2]},
                  x = obs, n = n)
## Visualization
plot(1:R, cov.sim, xlab = "Replication", ylab = "", pch = 16, cex = 0.7, col ="grey",
     ylim = quantile(cov.sim, probs = c(0.1, 0.9)))
lines(1:R, rep(cov.true, R), col = "green", lwd = 2)
lines(1:R, rep(cov.obs, R), col = "red", lwd = 2)
lines(1:R, cumsum(cov.sim)/(1:R), col = "blue", lwd = 2)
legend("topright", legend = c("Boot cov", "True", "Initial", "Boot average"),
       col = c("grey", "green", "red", "blue"),
       bg = "white", pch = c(16, NA, NA, NA), lwd = c(NA, 2, 2, 2))

The blue line corresponds to $\widetilde{\Sigma}_R$ as a function of the number of replications, and the red line is $\widehat{\Sigma}$.

enter image description here

$\endgroup$
8
  • 1
    $\begingroup$ I am not sure, but I dont think you could get only a confidence interval for my $\hat \Sigma$. I think it is also possible to estimate $\hat\Sigma$ through $\frac{1}{R}\sum_{i=1}^R\hat\Sigma^i$. If this is true and assume the theoretical backround for the bootstrap gives you the permission to take this estimation as a estimation for $\Sigma$, my question is which methode you would recommend, because I have some doubts that bootstrap would work on a sample size of 4 in this case. So I need a reason that this method will working. $\endgroup$ – Klaus Jul 15 '13 at 13:44
  • 1
    $\begingroup$ Thx that is what I am speaking about. Now I like to know is there are a bootstrap methode to get a blue line which comes closer the green line? $\endgroup$ – Klaus Jul 16 '13 at 11:08
  • 1
    $\begingroup$ This seems to be a new question, albeit somehow related to the initial question. It's probably best to ask a new question for this, and to close the current question if the answer provided is useful. $\endgroup$ – QuantIbex Jul 16 '13 at 11:20
  • 1
    $\begingroup$ But my first question under this subject is "Is the Bootstrap a good option to estimate Σ". If not, then I would like to discuss, if bootstrap is a good option to get the distribution of any statistic T which depends from this data. That means is the number of this sample enough to execute any bootstrap-methode to get the distribution of T. $\endgroup$ – Klaus Jul 16 '13 at 11:39
  • 1
    $\begingroup$ But my questions was is there any methode, to clearify that my examples wouldnt work was not the main question, I only post this, that the first comment is not the question to what I was trying so far. $\endgroup$ – Klaus Jul 16 '13 at 11:50

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.