# Combining sources of uncertainty/variation in a multi-layered linear model

I’m having trouble understanding how to combine parameter uncertainty and interannual variability from various levels in a bootstrapped linear model. Specifically, this model is designed to generate estimates of some phenomenon $Y$ for 0.5º grid cells in a spatial region (e.g, Central America).

It looks like this: $$\boldsymbol{Y_{m,y}} = \boldsymbol{\beta_{m,y}X_{m,y}}$$ where $\boldsymbol{\beta}$ is the vector of parameters being estimated (size 1x3), $\boldsymbol{X}$ is a matrix of observations of land cover fractions (size 3xN), $\boldsymbol{Y}$ is the vector of predicted values (size 1xN), and the subscripts indicate that the model is being parameterized for each month $m$ and year $y$ (e.g., May 2006, June 2006, and May 2007 can all have different parameter estimates). (No intercept, although I don’t think that matters for my question.)

If we think of the model as $Y_{m,y}=\sum\limits_{i=1}^3{\beta_{i,m,y}X_{i,m,y}}$, then each term $\beta_iX_i$ can be considered a subset of $Y$ — i.e., $Y_i$, the amount of $Y$ that happens on land cover type $i$ — in a real physical sense. Thus, there will likely be some correlation among the various parameter estimates. This is embodied/exacerbated by the fact that for each grid cell $g$, $\sum\limits_{i=1}^3{\boldsymbol{X_{g,i,m,y}}} = 1$ — that is, the land cover type fractions add up to 1.

To create uncertainty bounds for the parameters, I have done 10,000 bootstrapping runs. Let’s say there are 500 grid cells in the region. For the first bootstrapping run, I randomly choose 500 grid cells from that region with replacement. I then estimated the parameters for each month and year using that sample set. I then repeated this sampling-fitting procedure 9,999 more times.

My goal is to compare the measured/observed mean annual amount of phenomenon $Y$ with the estimated/modeled amount. I would like the uncertainty bars around the measured amount (which will take into account just interannual variability in the measured amount) to overlap with the uncertainty bars around the estimated amount (which will take into account both uncertainty in the model parameters as well as interannual variability in the annual estimated amount).

For the observations, it’s easy enough — sum up observed $Y$ for all the months in each year, then find the standard deviation across all the years. I get confused when thinking about how to do this for the estimates, though. The problem becomes even stickier if I want to compare observed vs. estimated annual GLOBAL $Y$ (i.e., across all regions).

I know that adding standard deviations in quadrature is going to be key here, but the fact that there are so many levels is confusing me. I'm working in Matlab now and am also fluent in R, but any help you can provide would be much appreciated.

• Could you please describe your grid and grid cells more in detail? When you write that there are 500 grid cells "in the region", what does that mean? Moreover, is the i you mention in your second formula running between 1 and 3? What does the sum on the grid cells (=1) means? Sorry for the avalanche of questions, but even though you explained a lot of details, the problem seems to be quite complex. Commented May 27, 2014 at 0:02
• No problem about the questions. I appreciate you trying to work through it! I've changed the question text to clarify the points you mentioned, but briefly, here are the answers. A region is something like Central America. I've overlaid a grid of 0.5º cells across the world, and so each region is made up of a number of these grid cells. You are correct about the i's in the summations, and they go from 1 to 3 because each X_i represents the fraction of each of three land cover types, which -- answering your third question -- sum to 1. Commented May 27, 2014 at 0:43

I have an idea for how to do this (this is Idea 1; see Idea 2 in a different answer). I won't mark this as the correct answer until I get some confirmation, so please let me know if this looks right to you (although make sure you understand what I'm trying to do!). My main point of concern is the way I include covariances in quadrature -- I'm afraid I am being redundant, and therefore the estimates of SD will be too large.

Note:

• $cov_b(i,j)$ means "the covariance of $i$ and $j$ across all bootstrapping runs $b$"

Find the mean and SD for a given cover ($c$) in a given month ($m$) of a given year ($y$) in a given region ($r$) across all B bootstrapping runs (e.g., forest in May 2006 in Central America):

$$\overline{Y_{c,m,y,r}} = mean_{b}(Y_{c,m,y,r,b}) = \frac{\sum\limits_{b=1}^B{Y_{c,m,y,r,b}}}{B}$$ $$\sigma_{c,m,y,r} = sd_{b}(Y_{c,m,y,r,b}) = \sqrt{\frac{\sum\limits_{b=1}^B{(Y_{c,m,y,r,b}-\overline{Y_{c,m,y,r}})^2}}{B}}$$ Find the mean and SD of total (i.e., across all 3 covers $c$) Y for a given month ($m$) of a given year ($y$) in a given region ($r$) (e.g., May 2006 in Central America): $$\overline{Y_{m,y,r}} = \sum\limits_{c=1}^3{\overline{Y_{c,m,y,r}}}$$ $$Y_{m,y,r,b} = \sum\limits_{c=1}^3{Y_{c,m,y,r,b}}$$ $$\sigma_{m,y,r} = \sqrt{\sum\limits_{c=1}^3{\sigma^2_{c,m,y,r}}+2\sum\limits_{c1=1}^2{\sum\limits_{c2=2\neq{c1}}^3{cov_b(Y_{c1,m,y,r,b},Y_{c2,m,y,r,b})}}}$$ Find the mean and SD of total Y for a given month ($m$) in a given region ($r$) (e.g., for all Mays in Central America). $$\overline{Y_{m,r}} = mean_y(\overline{Y_{m,y,r}}) = \frac{\sum\limits_{y=1}^N{\overline{Y_{m,y,r}}}}{N}$$ $$Y_{m,r,b} = \frac{\sum\limits_{y=1}^N{\sum\limits_{c=1}^3{Y_{c,m,y,r,b}}}}{N}$$ $$\sigma_{m,r} = \frac{1}{N}\sqrt{\sum\limits_{y=1}^N{\sigma^2_{m,y,r}}+2\sum\limits_{y1=1}^{N-1}{\sum\limits_{y2=2\neq{y1}}^N{cov_b(Y_{m,y1,r,b},Y_{m,y2,r,b})}}}$$ Find the mean and SD of annual total Y for a given region ($r$) (e.g., for Central America): $$\overline{Y_{r}} = \sum\limits_{m=1}^{12}{\overline{Y_{m,r}}}$$ $$\sigma_r = \sqrt{\sum\limits_{m=1}^{12}\sigma^2_{m,r}+2\sum\limits_{m1=1}^{11}{\sum\limits_{m2=2\neq{m1}}^{12}{cov_b(Y_{m1,r,b},Y_{m2,r,b})}}}$$

• I think what you did makes sense, but I'll let it to more experienced people to comment. Yet, I'll have to stress again the same issue, I am sorry, as I repeat myself and this is certainly heavy, but I think the way you bootstrap is not clear. You write that each grid has 500 regions; and then that you draw 500 regions with replacement from this. So this means you will for sure repeat some values but I am not sure this provides a correct assessment of the grid intrinsic variability, and this is central to the deduction that follows. Commented May 28, 2014 at 16:14
• Yet, I think the question is well-formed and I will vote it up (in the hope it will attracts others attention as well). Commented May 28, 2014 at 16:15
• Technically, I'd say the (fake) REGION in my example has 500 GRID CELLS. But yes, it's a good question as to whether the sampling procedure adequately assesses the variability within the region. Thanks for the upvote and for all your help in clarifying the problem. Commented May 28, 2014 at 16:26

I am not sure I understood the problem, so I'll star by trying to resume all important features for the discussion. Correct me if wrong, I'll edit the answer (or someone will) and we will work on it iteratively.

You have a linear model that predicts the occurrence some phenomenon Y of which you want to estimate variability. This phenomenon is calculated for a month $m$ and year $y$ and the input parameters of it are just 3: the fraction of terrain of each of the three types covered by the geographic grid. These fractions add up to one, naturally. Your objective is to estimate the inter-annual variability as measured, as estimated by the model, and to compare these with the variability that is intrinsically associated to the model.

The first thing I find hard to understand is the concept of "region". When you draw from a region, my understanding is that you are drawing groups of $\beta_i$ such that these are valid ($\Sigma_i \beta_i=1$) and that somehow are still equivalent (maybe in the error bars) of the original observed grid point. I do not understand if $Y_{m,y}$ only makes sense for the context of a wide-scale value calculated for all the grid points, resulting from the N observations on N grid points, or if it actually that runs on the gridpoints as well. In other words, I do not know what N refers to and what are the dependent variables of your model.

I think your method of calculating variability through bootstrap is justified, and you could go even for more "raw" methods, more conservative but that make less assumptions, like the jacknife, but I think before you should try to present your problem in a clearer way.

• Your understanding of the situation is mostly correct. The point where we diverge is the end target here. I want to compare the MEASURED/OBSERVED mean annual amount of phenomenon Y with the ESTIMATED/MODELED amount. I would like the uncertainty bars around the measured amount (which will take into account just interannual variability in the measured amount) to overlap with the uncertainty bars around the estimated amount (which will take into account both uncertainty in the model parameters as well as interannual variability in the annual estimated amount). I'll add another comment re:regions. Commented May 27, 2014 at 16:22
• N refers to the number of grid cells in the region (500 in my fake example). Note that by drawing (with replacement) a random subset of grid cells, I'm not drawing groups of Beta_i (the parameters), I'm drawing groups of X_i (the observed land cover fractions). The optimum Beta_i values (as determined by minimum sum of squared errors) for that sample are then calculated. Y_{m,y} does make sense in terms of an individual grid cell -- it's the amount of Y that happens on that grid cell. Let me know if this makes sense. Thanks again for your patience; it's always clearer in one's head, isn't it? Commented May 27, 2014 at 16:30
• It's getting clearer, no need to apologize. =) I would advise you to add this extra information of what you want to compare directly to the main text, paying attention to state the question clearly as you wrote in your first comment. I understand now my mistake on the interpretation of the $\beta_i$'s. Yet, I still have to say you did not answer me as to where you drew your $X_i$'s from. This is fundamental to understand the validity of the bootstrapping procedure. Commented May 27, 2014 at 16:54
• Let's say a region has $N$ grid cells. For each bootstrapping run, I randomly choose $N$ grid cells from the region with replacement. Then for each month and year, each grid cell $j$ has four measurements associated with it: The different land cover fractions present there ($X_{j1}$, $X_{j2}$, $X_{j3}$) and the amount of the phenomenon observed there ($Y_j$). I use these to estimate the regional $\beta_1$, $\beta_2$, and $\beta_3$ for the month+year in question. I do this for the same sample set for all other months and years. Then I start over, choosing a new sample set. Does that help? Commented May 27, 2014 at 17:38

I've come up with an alternative that is much simpler (this is Idea 2; see Idea 1 in a different answer). I still don't know if this is right, though! This explanation is for a given region. I will be using Matlab code to show what I'm doing.

Definitions:

• Z = number of years modeled
• B = number of bootstrapping runs
• estY_zbmc = $Z \times B \times 12 \times 3$ array, with model-estimated $Y$ on each land cover type $c$ in each month $m$ (e.g., July or August) in each year $z$, for each bootstrapping run $b$.
• obsY_z = $Z \times 1$ array of total $Y$ (i.e., for all land cover types combined) for each year.

Here's the code:

estY_zb = sum(sum(estY_zbmc,4),3) ; % Find total estimated Y for each year in each bootstrapping run
estY_z_mean = mean(estY_zb,2) ; % Across all bootstrapping runs, find mean total estimated Y for each year
estY_z_sd = std(estY_zb,2) ; % Across all bootstrapping runs, find st. dev. of total estimated Y for each year
estY_meanAnn = mean(estY_z_mean) ; % Find mean annual total estimated Y


That gives us the mean annual Y easily enough. What we want next is to find the error bars around that number. We'll have to do something like this:

$\sigma = \frac{1}{Z} \sqrt{\sum\limits_{z=1}^Y{\sigma_z}+2\sum\limits_{z1=1}^{Z-1}{\sum\limits_{z2=2 \neq z1}^Z{\sigma_{z1,z2}}}}$

The expression under the radical is just the sum of all the elements in the covariance matrix, so I can do the following (transp used because cov(X) treats rows of matrix X as observations, and columns as variables):

cov_z1z2 = cov(transp(estY_zb)) ; % transp() used because cov(X) treats rows of matrix X as observations, and columns as variables
sum_cov_z1z2 = sum(cov_z1z2(:)) ;


And now:

estY_sd = sum_cov_z1z2 / Y ;


I can then compare estY_meanAnn and estY_sd with the mean and SD of mean annual observed $Y$:

obsY_meanAnn = mean(obsY_z) ;
obsY_sd = std(obsY_z) ;