When it is a good idea to average data to reduce variability? I have 4 responses per day from animals, at different points of the day (morning, afternoon, etc).
My data is marred by a lot of variability, since the nature of the response is such that while some of it is systemic, a lot of it is just random behavior by the animals.
I fit an initial regression model which was very flexible (all continuous covariates were fitted with smoothers), and I only manage to explain at most 35-40 % of the deviance. 
However, if I make a boxplot of the response dependent on the time of day, I get 

So, it does not seem to have an affect. So, is all of the above a good enough reason to say: you know what, let's just average the per day data to reduce some of the variability, and thus ignore the "time"-factor, since it does not have an affect anyways?
 A: Averaging your data can only hurt in this case.
Although you could average the samples across each day, this won't provide any benefit. A standard linear regression model aggregates the data automatically. Averaging the samples for each day would decrease the variance but would also reduce the total number of observations, resulting in no net gain in power.
Furthermore, averaging time points might actually lose data.
I've wondered about this before myself so I made a toy example in R, and simulated data similar to your problem's setup.
Example code
noise = 2

sim1 <- function(seed, noise, control=FALSE, aggregate=FALSE)
{
  # Set up data
  day = c(1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4)
  set.seed(seed) # For reproducibility

  # Simulate noisy observations with mean 50
  obs = rep(50, length(day)) + rnorm(length(day), sd=noise)

  # Don't add true signal if it's a control simulation
  if (!control)
  {
    obs = obs + day  
  }

  # If FALSE, run method 1. If TRUE, run method 2.
  if (aggregate)
  {
    day = c(1, 2, 3, 4)
    obs = c(mean(obs[1:4]), mean(obs[5:8]), mean(obs[9:12]), mean(obs[13:16]))
  }

  fit = lm(obs ~ day)

  # Return TRUE if slope of line differs from 0 with p-value < 0.05
  return(summary(fit)$coefficients[2,4] < 0.05)
}

# Positive control: aggregation has less power to detect true positives
print(mean(sapply(1:100000,sim1, noise=noise, aggregate=FALSE)))
print(mean(sapply(1:100000,sim1, noise=noise, aggregate=TRUE)))

# Negative control: both methods maintain power of ~ 0.05.
print(mean(sapply(1:100000,sim1, noise=noise, control=TRUE, aggregate=FALSE)))
print(mean(sapply(1:100000,sim1, noise=noise, control=TRUE, aggregate=TRUE)))

Method 1: Model each data point separately.
We get a power (true positive rate) of 0.54882, and a false positive rate of 0.04961 as desired.
Method 2: Average points for each day.
We get a power (true positive rate) of 0.25541, and a false positive rate of 0.05038 as desired.
Method 1 significantly outperforms method 2 on these simulated data. Maybe someone else can give a more rigorous mathematical explanation of why this is true.
You can probably exclude the time of day from your regressors, as long as it's not interacting with any of the other variables to produce effects that are concealed right now. However, it will be impossible to explain more of the variance by excluding the time of day; removing variables can only decrease a model's performance on the data used to train it.
Hope this helps :)
