Good data example needed with covariate affected by treatments

Question

I have looked at a lot of R datasets, postings in DASL, and elsewhere, and am not finding very many good examples of interesting datasets illustrating analysis of covariance for experimental data. There are numerous "toy" datasets with contrived data in stat textbooks.

I'd like to have an example where:

• The data are real, with an interesting story
• There is at least one treatment factor and two covariates
• At least one covariate is affected by one or more of the treatment factors, and one is not affected by treatments.
• Experimental rather than observational, preferably

Background

My real goal is to find a good example to put in the vignette for my R package. But a larger goal is that people need to see good examples to illustrate some important concerns in covariance analysis. Consider the following made-up scenario (and please understand that my knowledge of agriculture is superficial at best).

• We do an experiment where fertilizers are randomized to plots, and a crop is planted. After a suitable growing period, we harvest the crop and measure some quality characteristic - that's the response variable. But we also record total rainfall during the growing period, and soil acidity at time of harvest -- and, of course, which fertilizer was used. Thus we have two covariates and a treatment.

The usual way to analyze the resulting data would be to fit a linear model with the treatment as a factor, and additive effects for the covariates. Then to summarize the results, one computes "adjusted means" (AKA least-squares means), which are predictions from the model for each fertilizer, at the average rainfall and the3 average soil acidity. This puts everything on an equal footing, because then when we compare these results, we are holding rainfall and acidity constant.

But this is probably the wrong thing to do -- because the fertilizer probably affects the soil acidity as well as the response. This makes the adjusted means misleading, because the treatment effect includes its effect on acidity. One way to handle this would be to take acidity out of the model, then the rainfall-adjusted means would provide a fair comparison. But if acidity is important, this fairness comes at great cost, in the increase in residual variation.

There are ways to work around this by using an adjusted version of acidity in the model instead of its original values. The upcoming update to my R package lsmeans will make this downright easy. But I want to have a good example to illustrate it. I will be very thankful to, and will duly acknowledge, anyone who can point me to some good illustrative datasets.

Answer 1

You may want to check out the mediation R package. It does include experimental data like jobs and framing where the treatment variable affects both a response variable and covariates (i.e., mediators of the treatment effect), along with covariates not affected by the treatment.

I looked into the mediation literature because I though you exactly described a mediation study: the fertilizer effect on the crop quality is mediated through its effect on soil acidity. Even if the datasets in the mediation package do not satisfy you, you may find one if you look into the mediation literature.

Answer 2

I thought I'd show how an analysis comes out with one of the datasets in the mediation package. In framing, an experiment is done where subjects have the opportunity to send a message to Congress regarding immigration. However, some subjects (treat=1) were first shown a news story that portrays Latinos in a negative way. Besides the binary response (whether or not they sent a message), we also measured emp, the subjects' emotional state after the treatment was applied. There are various demographic variables as well.

First, let's load the needed packages in R, and change the labels for educ to shorter strings.

> library("lsmeans")
> library("mediation")
> levels(framing$educ) = c("NA","Ref","< HS", "HS", "> HS","Coll +")

Now fit a logistic regression model

> framing.glm = glm(cong_mesg ~ age + income + educ + emo + gender * factor(treat),
+                   family = binomial, data = framing)

Here is a display of the conventional adjusted means, where predictions are made with the covariates age, income, and emo set at their mean values:

> lsmip(framing.glm, treat ~ educ | gender, type = "response")

This is a curious result because the displayed treatment effects are the opposite for females as for males, and the effect of education isn't monotone as one might expect.

Note, hHowever, emo is a post-treatment measurement. This means that the treatment could have affected it, i.e. emo is a mediating covariate; and so it may not be meaningful to compare predictions of the response variable while holding emo constant. Instead, let's look at the predictions where emo is set to its predicted values given treat and the demographic variables.

> lsmip(framing.glm, treat ~ educ | gender, type = "response",
+       cov.reduce = emo ~ treat*gender + age + educ + income)

This result is quite different, suggesting that emo plays a strong mediating role. (The mediation package has functions for estimating the strength of these effects.) The above predictions suggest that, taking emotional response into account, male subjects exposed to the negative news story are more likely to send the message than are females or those not seeing the negative news story. Also, the effect of educ is (almost) monotone.

Thanks again to @MasatoNakagawa for pointing me to this interesting example and attuning me to some recent research on causality.

Answer 3

Look up gene-environment interaction GWAS studies. The statistical analysis they perform in essence is what you have described. The question is does your environment matter to a phenotype (observable feature)? One school of thought generally ignores all environmental information and says your genetic makeup describes your phenotype. This is in complete contrast with ecological studies where the story is environment is everything and they ignore the genes. Since both parties are trying to understand the same problem, there have been recent attempts to coalesce the two.

Say we are studying BMI. We take the first few principal components of the genetic matrix as the fixed effects due to genes. We fit education with an index 1 for well educated and 0 for poorly educated as a fixed effect. There is a reasonably strong correlation between the education index and wealth of the community the person is from. So one would argue that the low income communities are more likely to have more fast food restaurants. The fast food acts as an obesogenic trigger.. "Triggers something in your genetic setup which encourages fat accumulation" so it will show up in the genetic makeup in some form.

Simulating such data is not a problem. Look up

http://pngu.mgh.harvard.edu/~purcell/plink/simulate.shtml

This lets you simulate GWAS (think of this as genetic units) data responsible for a symptom. If not instructed otherwise it will generate 1000 with the symptom and 1000 controls. The norm in these simulations which I use is 9990 SNPs do not cause the symptom and 10 SNPs do. Read the instructions on how these are simulated.

The output will be 1 if the person is obese and 0 if he isn't. Simulate education factors (finished college eduction/not finished college education )based on some reasonable correlation with obesity levels.

Hope this helps!!!

Answer 4

I'd recommend reading Freakonomics, and finding the papers their work is based on, and seeing if you can grab that data. They have some really interesting work on really interesting datasets, and in some cases they figure out very clever ways to test hypotheses despite limitations in the data.