Consider the following ridiculous example. There is a real, non-ridiculous research question here, I promise -- I'm just a little bit uncomfortable posting what I'm working on on the internet in great detail. (And the below is hopefully a bit more entertaining):

A race of intelligent space monkeys has come to earth, liberated regular earth-monkeys from zoos, and taken their place. These intelligent space-monkeys are actually graduate students at a space-monkey university, and are conducting an experiment on behalf of a space-monkey professor: they wish to estimate how the happiness of human zoogoers changes when excrement is thrown at them. (This is an important and unsettled problem in extraplanetary psychology, and Prof. Spacemonkey is aiming for a pub in the prestigious Ann. Improb. Res.) They are also interested in what moderates these observed changes, in the interests of external validity and hopefully better out-of-sample prediction. Importantly, these space-monkeys are able to create forks in space-time, whereby they may both apply and not apply the treatment (flung excrement) to their sample onlookers. Furthermore, each cohort of zoo-space-monekys creates $m \geq 2$ parallel universes, one of which always has zero excrement, and in which the remainder have excrement flung at different intensities.

Thus, the dataset is comprised of experiments -- comparing treatments at varying intensities against identical untreated controls. We observe the average level of happiness $H^{\text{excrement}}$ of those subjected to the treatment at a given intensity $t \in T > 0$, and the average level of happiness $H^{\text{No excrement}}$ from the same population in the control universe, where $T = 0$. In addition, he observes audience characteristics, as well as summary statistics of audiences $\mathbf{X}$.

Having been trained as an ecologist before moving into extraplanetary psychology, he is familiar with the representation of this type of data in the form of response ratios: $RR \equiv log\left(\frac{H^{\text{Excrement}}}{H^{\text{No excrement}}} \right)$.

Collecting the data, he immediately specifies a meta-regression model: $$ RR_{iz} = \alpha + \tau_z + f\left[T_{iz}, \mathbf{X}_z \right] +\epsilon_{iz} $$ where $\alpha$ is an intercept, $\tau$ is zoo-level random effect, $i$ indexes a particular treatment control comparison, and $z$ indexes zoos.

Prof. Space Monkey is using mgcv, so he codes

mod = gam(RR~ 
     s(zoo, bs = "re")
+    s(T)
+    s(x1)
+    s(x2)
+    te(T,x1,x2)

Let's simulate some data for Prof. Space Monkey:

zoo = ceiling(runif(1000)*10)
T = runif(1000)*10
x1 = rnorm(1000)
x2 = rnorm(1000)
RR = -.4*T + .0005*T^2 + .69*sin(x1) + .042*exp(x2) + T*x1*x2*.314 + zoo*.5 +rnorm(1000)*5

After fitting the model, Prof Space Monkey is dismayed to discover that the expectation of $RR$ is not zero when $T=0$!

vis.gam(mod,view=c("T","x2"), plot.type="contour")

enter image description here

This doesn't make any sense! By construction, the dependent variable should be zero when there is no excrement throwing! Prof Spacemonkey is distraught. He can't just remove the intercept on the regression -- this would bias things, and wouldn't even solve the problem, because x1 and x2 would still cause $\hat {RR}$ to be nonzero even at T=0.

Tossing and turning in his space-bed one night, Our Hero has an idea! He'll simply interact the treatment (excrement-throwing) with the entire right-hand-side of the model:

$$ RR_{iz} = f(T_{iz}) + f(T_{iz}) \times \tau_z + f\left[T_{iz} \times\mathbf{X}_z \right] +T_{iz} \times \epsilon_{iz} $$

because the constant = 1.

Great! The intuition of the model is now a little different: Response ratio is a function of excrement intensity, the effect of which is mediated by characteristics of individual zoo-experients. These zoo-level characteristics do not have main effects because of the nature of the dependent variable as a ratio. Some of these zoo-experiment characteristcs are observable and modeled as fixed effects, while some aren't, and enter as random effects.

Importantly, the expectation of $RR$ is now zero when $T$ is zero. The model is also now heteroskedastic, so maybe a wild bootstrap is needed for confidence intervals. (?). But he has four worries:

  • How does he code this in R? The random effect now has a level of "intensity" that is dictated by the intensity of thrown excrement.
  • Generalized additive models constrain smooth terms to sum to zero -- but that is incompatible with the constraint that $E[RR]\rightarrow 0$ as $T\rightarrow 0$. What to do?
  • Is a wild bootstrap indeed the way to go to capture the heteroskedastic nature of the model?
  • As smart as this idea seems, Prof Space Monkey is really worried that he might be making some error somewhere. Unknown unknowns and all that.

Suggestions very much appreciated.

  • 6
    $\begingroup$ It is gratifying that superintelligent space monkeys also sometimes have trouble with R. $\endgroup$
    – Glen_b
    Apr 15, 2013 at 4:38
  • $\begingroup$ This is a bit late, but are you sure $E(RR|T=0)$ is zero? In your code, you specify RR = ... + zoo*.5 and zoo is an integer ranging from 1 to 10. $\endgroup$
    – Hong Ooi
    Mar 16, 2015 at 4:36

1 Answer 1


A solution has occurred to me: rather than interacting the entire RHS with the intensity variable, I could simply specify the model as $RR/f(I) = whatever$. This is to say that RR per unit $I$ is a function of my explanatory variables. Of course, I don't know RR per unit $I$. But if I did, I could set them all equal to the same $I$, and have $RR/f(I=c)$. So I need an estimate of the response ratio at an intensity defined by a common constant.

I estimate this using a random coefficients model with no intercept, and all terms interacted with the intensity variable. I then use the estimated random effect predictions and the coefficient estimates to get $RR_{I=c}$. I then fit the model without $I$ in it at all.

Some sort of multiple imputation approach is probably appropriate given the uncertainty in estimating $RR$ as a function of $I$


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.