The PCB in Fish Example: Simple Linear Regression Model

The PCB in Fish Example: Simple Linear
Regression Model

1 Data

Data used here are PCB concentrations in lake trout collected by the Wisconsin Department of Natural Resources from 1974 to 2003 (Figure 1). The PCB concentration – fish size relationship (Figure 2)
represents the biological accumulation of PCB over time, as a larger fish is likely to be older.

Figure 1: Temporal trend of fish tissue PCB concentrations – PCB concentrations in lake trout from Lake Michigan decline over time, but shown a stabilizing trend in the last few years.

Figure 2: Fish tissue PCB concentrations vs. fish length – Large fish tend to have higher PCB concentrations in lake trout from Lake Michigan.

2 Regression with One Predictor

The first order rate model suggests that a simple
linear regression be used for assessing a temporal trend is a log linear model:

                                            (1)

The model coefficients are estimated using the least squares method, implemented in R function lm():

lake.lm1 <- lm(log(pcb) ~ year, data=laketrout)
display(lake.lm1, 3)
 
lm(formula = log(pcb) ~ year, data = laketrout)
(Intercept) 119.8467  10.9689
year         -0.0599   0.0055
---
n = 631, k = 2
residual sd = 0.8784, R-Squared = 0.16

The estimated (the intercept) is 119.85 and the estimated (the slope) is -0.06. With these two coefficients, we can calculate the mean log PCB concentration for a given year: . The estimated residual standard deviation of 0.8784 describes the variability or uncertainty. When putting the two parts together, the fitted model can be seen as a conditional normal distribution describing the probability distribution of log PCB concentrations. For example, the estimated log PCB distribution for year 1974 is or .

3 Model Interpretation

3.1 Centering the Predictor

The intercept of a simple regression model is the expected value of the response variable when the predictor is 0. For this model, we don’t believe that the model can be extrapolated to year 0. Consequently, the intercept cannot be interpreted to have any physical meaning. However, if the model is refit with using as the new predictor, the new intercept is 1.66, the mean log PCB concentration of 1974. The transformation , a linear transformation, does not change the fitted model, but the resulting intercept is easier to interpret.

3.2 Slope

The slope is the change in log PCB for a unit change in year. Because the response variable is log PCB concentration, a change of in the logarithm scale is a change of factor of in the original scale. That is, the initial year (1974) concentration is . The second year (1975) PCB concentration is , or . Given , the 1975 concentration is approximately 6% less than the 1974 concentration. The slope is the annual rate of reduction.

3.3 Residuals

The residual or model error term describes the variability of individuals. For this model, the estimated residual standard deviation is 0.87. When interpreting the fitted model in the original scale of PCB concentration, the predicted PCB concentration has a log normal distribution with log mean and log standard deviation 0.88. This model suggests that the middle 50% of the PCB concentrations in 1974 will be bounded between or (2.74, 8.97) mg/kg, and the middle 95% of the concentration values are bounded by (0.88, 27.79) mg/kg. The estimated mean concentration in 1974 is mg/kg, and the estimated standard deviation is, or , 1.081 times of the mean (i.e., the coefficient of variation cv = 1.081).

The model can be summarized graphically as in Figure 3.

Figure 3: Simple linear regression of the PCB example – PCB concentration
data are plotted against year. The simple linear regression resulted in highly
uncertain predictions. The solid line is the predicted mean PCB concentration
and the dashed lines are the middle 95% intervals.

The PCB in Fish Example: Multiple Regression -- Interaction

When fitting the multiple regression model with yr and len.c as the predictors, an important assumption is that the effect of year (the slope of year) is not affected by the size of the fish and the effect of fish size (the slope of length) is the same throughout the study period. This is the additive-effect assumption imposed on a multiple regression model. Is this assumption reasonable? Madenjian et al. [1998] reported that small lake trout ( cm) eat small alewives (Alosa pseudoharengus, which have an average PCB concentration of 0.2 mg/kg), intermediate-size lake trout (40 ~ 60 cm) eat alewives and rainbow smelt (Osmerus mordax, whose PCB concentrations ranged from 0.2 to 0.45 mg/kg) and large lake trout (≥ 60 cm) eat large alewives (with an average PCB concentration of 0.6 mg/kg). On the one hand, because larger fish tend to consume food with higher concentrations of PCB, its reduction over time should be slower than the rate of reduction of small fish. On the other hand, because PCB was banned in the 1970s, the natural reduction of PCB through microbiological metabolism resulted in the overall reduction of PCB concentration in the environment and in fish. We expect that the PCB – length relationship will change over time. In other words, the slope of year in the multiple regression model is expected to change with the size of a fish and the slope of length is expected to change over time. To model this “interaction” effect, we add a third predictor, the product of yr and len.c in the model:

#### R code ####
lake.lm4 <- lm(log(pcb) ~ I(year-1974)*len.c, data=laketrout)
display(lake.lm4, 4)
 
#### R output ####
lm(formula = log(pcb) ~ I(year - 1974)*len.c, data = laketrout)
                     coef.est coef.se
(Intercept)           1.8967   0.0465
I(year - 1974)       -0.0873   0.0036
len.c                 0.0510   0.0038
len.c:I(year - 1974)  0.0008   0.0003
---
n = 631, k = 4
residual sd = 0.5520, R-Squared = 0.67

When the interaction term is included, the model is expressed as:

(1)

Because of the product term, the model is no longer a linear model. The slopes of centered length (len.c) and year (yr) are no longer constant. We can rearrange the model to understand the interaction effect. First, the interaction term is grouped with yr:

That is, the effect (or slope) of is now a function of . The slope shown (-0.087) is the slope when [;Len.c = 0;] or the year effect for an average sized fish. When the fish size is 10 cm above average, the yr effect is -0.087 + 0.00085 ⋅ 10 = -0.0785. In other words, not only a larger fish has a higher PCB concentration on average, PCB in a larger fish tend to dissipate at a lower rate. This interpretation is true only when we are comparing same-sized fish over time. So, when comparing fish of the average length (Len.c = 0), the annual rate of dissipation is 8.7%. The annual dissipation rate is 7.6% for fish with a size 10 cm above average. When examining the log(PCB) fish length relationship, the model can be rearranged to be:

The relationship is still linear for any given year. But the slope changes over time. Initially, (yr = 0 or 1974), the size effect is 0.051. Each unit (1 cm) increase in size will result in a 5.1% increase in PCB concentration.  Ten years later (1984), the slope was 0.051 + 0.00085 ⋅ 10 = 0.0595. The size effect is stronger. This is reasonable because the rate of concentration decreasing for a large fish is smaller than the rate for a small fish. Consequently, the difference in concentration between the same two fish increases over time.

The interaction effect is small (albeit statistically significant). Can this small interaction effect be practically significant? Because the response variable is in logarithmic scale, we need to be careful in interpreting a small effect. For the slope of yr, the slope value for a small fish (-6.7 cm below average, or the first quartile) is 0.09 - 0.00085 × (-6.7) = 0.095 and the slope is 0.09 - 0.0008 × (8.5) = 0.083 for a large fish (8.5 cm above average, the third quartile). PCB concentration reduction is at a lower rate (~ 8%) for a large fish and a higher rate (~ 10%) for small fish. The slope of len.c increases from 0.05 in 1974 to 0.074 in 2004, indicating a much larger difference in PCB concentration between a large and a small fish.

References

C.P. Madenjian, R.J. Hesselberg, T.J. Desorcie, L.J. Schmidt, Stedman. R.M., L.J. Begnoche, and D.R. Passino-Reader. Estimate of net trophic transfer efficiency of PCBs to Lake Michigan lake trout from their prey. Environmental Science and Technology, 32:886–891, 1998.

Bootstrapping for threshold standard error

In my 2003 Ecological Modelling paper on ecological threshold, I used the bootstrapping method for estimating the standard error of a change point variable.  The change point was estimated as the binary break point that resulted in the greatest reduction in deviance.  Ian McKeague of FSU (now at Columbia) told me about Buhlmann and Yu (2002), which suggested that bootstrapping is inappropriate for a change point problem.  I was convinced that the bootstrap confidence interval is wrong, but was never able to explain the reasons convincingly until recently.

A change point problem limits the potential break points to be the gaps (intervals between distinct values) in the given data set.  These potential break points do not change from bootstrap sample to bootstrap sample, only that some gaps may become wider.  Should the population were sampled, the potential number of break points is infinite.  As a result, bootstrap estimated break point is of an artificially smaller sampling variance than the variance of the true sampling distribution of the break point. This reduced sampling variation is likely the cause of a much smaller standard error or narrower confidence interval when estimated using bootstrapping.  In Banerjeer and McKeague (2007), a simulation is performed to show that a 95% confidence interval based on bootstrapping covers the underlying true change point far less than 95% of the time.

Causal inference based on observational data

I have been reading causal inference using observational data for almost three years.  Three years ago, I read a paper on the effectiveness of agricultural best management practices on reducing nutrient loading from fields. The data were observational, and the estimated effects are affected by confounding factors.  Using that data set, I learned how to use multilevel model and made a few homework assignments for my class.  The question in that data set is clear and the approach of using the propensity score match is obvious.  I did a comparison of using the propensity score matching and multilevel modeling for estimating the effect of conservation practices on nutrient loading from agriculture fields. The results were not surprising -- on average, any conservation practice would lead to a reduction of nitrogen/phosphorus loading by 1/3-2/3.  I have since read Rosenbuam and Rubin (1983) and other classical literatures.  To summarize the main idea, the problem of causal inference of a treatment effect is always the problem of counter factual.  That is, we can never apply a treatment and control to the same subject to assess the real treatment effect on a single subject. When treatment and control were applied to two different subjects, the result can be due to confounding factors. Statistical solution to this problem is Fisher's randomization. Subjects are randomly applied control and treatment, effectively forming two groups that are otherwise similar and one applied treatment and the other control.  The difference in response can be confidently attributed to as the treatment effect.   When using observational data, treatment and control are not applied randomly.  Consequently, no causal inference should be made without careful adjustment. The propensity score matching approach is one of many methods used for causal inference based on observational data.  The basic idea is that data can be subset such that the treatment and control groups are similar in terms of confounding factors. The similarity is achieved though the propensity score, or the likelihood a subject receiving treatment given its confounding factor levels. Operationally, a propensity score is estimated as the probability of receiving treatment through a logistic regression using all available covariates.  Each observation in the treatment (or control) group is matched with an observation in the other group with a similar propensity score.  Not all observations will be included in the final data set.  After matching, we essentially have a treatment group and a control group that are similar with respect to confounding factors.  The difference in response is likely due to the treatment. 

As central idea of using the propensity score is matching and subset to balance out the effects of confounding factor. In environmental studies, the treatment is often a continuous variable. For example, the dose-response curve, the effect of nutrient on lake eutrophication.  I read Imai and van Dyk (2004) and Hirano and Imbens (2004) recently and learned the generalized propensity score approach for continuous and other forms of treatments.  The generalized propensity score is actually a natural extension of the binary treatment one. The propensity score of a binary treatment variable is the probability of an observation receiving treatment, or the probability distribution function of the binomial distribution of treatment.  The generalized propensity score is the probability density of treatment conditional on confounding factors. One can estimate this probability distribution by fitting a regression model if the treatment over confounding factors. The resulting regression model forms a conditional normal distribution of the treatment variable, from which the propensity scores can be calculated (as the density of the model residuals (~N(0, sigma)). Imai and van Dyk (2004) suggested that the resulting propensity scores be used to divide the observations into J groups of similar propensity scores, and within each group a regression model is fit for the treatment effect. The average causal effect is a weighted mean of the J effects.  Hirano and Imbens (2004) suggested a two step process: (1) the response is regressed to the polynomial of treatment and the propensity scores, and (2) the causal effect is estimated as the expected response of a given treatment value over all observations.

In any case, the propensity score approach still aims at achieving a balance of the confounding factors. When observations are divided into groups based on the estimated propensity scores (Imai and van Dyk, 2004), distributions of confounding factors (coveriates used for estimating propensity scores) are more less balanced.

A very odd application of the generalized propensity score is published recently in Ecological Applications (Yuan 2010).  The paper misinterpreted the definition of the propensity score.  Instead of the probability density of an observation, Yuan (2010) used the predicted mean treatment value as the propensity score. The resulting subgrouping further segregated the distributions of confounding factors, leading to a potentially more biased estimate of the causal effect.