# Calculate mean of one variable while controlling for another using regression

This is probably a very naive question... I'd like to estimate "adjusted" or "conditional" means for a variable (i'm unsure of the correct terminology). My data are on cortisol levels (dependent variable) in rabbits (n=56). I have many measurements at different times of the day, over many months. I'd like to calculate mean weekly values of cortisol for each individual rabbit so these can be used as a predictor in another model for which I only have weekly data. Rather than calculate the means from the raw data, i'd like to control for the time of day the samples were taken (this can influence the measurement). I thought i'd regress time of day (in minutes from 00:00 each day) on cortisol level and then extract the fitted values and calculate the weekly mean for each rabbit from these. Would this give me the estimated mean for cortisol, while controlling for time of day?

I can't share my data, but i've created a similar mock up using the iris data set. Here I fit a model, extract the fitted values and then calculate "adjusted" means for each species while controlling for the predictor. Am I right in thinking the difference between these means and the ones for the raw data (below) reflect the adjustment made when controlling for the independent variable?

data(iris)

fit <- lm(Sepal.Length ~ Petal.Length, data = iris)
summary(fit)

with(iris, plot(Sepal.Length ~ Petal.Length, col = as.numeric(Species), asp = 1))
abline(coef(fit))

iris\$fitted <- fitted(fit)

with(iris, aggregate(fitted, list(Species), mean))
#      Group.1       x
# 1     setosa  4.9044
# 2 versicolor  6.0486
# 3  virginica  6.5769

with(iris, aggregate(Sepal.Length, list(Species), mean))
#      Group.1      x
# 1     setosa  5.006
# 2 versicolor  5.936
# 3  virginica  6.588


Using fitted means instead of actual means is risky. You have to be able to really trust your model. If your model is mis-specified, then the fitted means may be badly biased.

If you feel like you really need to control for time of day, though, I would try fitting a non-linear effect for time of day. If you're fitting a non-parametric effect for the time of day, then the chance of model mis-specification may be lower, and so you may get less bias.

For example, you could try fitting a GAM:

library(mgcv)

gam.fit = gam(Weight ~ Rabbit + te(Time.Of.Day, bs="cc"), data=Rabbit.Data)

This command will fit a periodic spline effect to the time of day (you can look at it using plot). You can then treat the coefficients for each rabbit as means that have been controlled for time of day effects.

• Yes, and then estimate mean cortisol for each rabbit at the median value of time (or any other point you wish to choose, perhaps for subject matter/theoretical reasons). Apr 25, 2013 at 2:28
• @Stefan Wager: thanks for your advice. I've been thinking more about this. Would it not make sense to use the residuals from the model with just time of day as the predictor and calculate the mean of those for each rabbit by each week? Apr 25, 2013 at 23:20
• @Chris: For the model I wrote down, the model residuals will be mean-zero for each rabbit, because each rabbit has its own factor. So the two methods would be equivalent. Apr 26, 2013 at 2:54
• @Stefan Wager: ok, yes, I can see that. What about just fitting the model I originally proposed - an OLS model without rabbit as a predictor (cortisol ~ time.of.day) - and then calculating the mean of the residuals for each rabbit by each week? Would that still suffer from the bias issues you mentioned for the fitted values? Apr 26, 2013 at 4:24
• @Chris: Yes, you should include the rabbit predictors. Otherwise, you will have bias problems unless your groups are perfectly balanced. Leaving out the rabbit factor can only hurt you, but not help you. Apr 26, 2013 at 7:39