What are good methods to deal with outliers when calculating mean of data? I have a dataframe with yearly energy uses of buildings over 5 years. In order to have a representative yearly energy use for data modelling, I'll have to take the mean of those data. As the data can contain outliers, I want to deal with outliers correctly (but keeping as much proper data as possible). (The df can contain empty cells (and older years are a little more important than recent years if a weighted decision has to be made somewhere).)
What are good methods to deal with outliers when calculating the mean of data?
I thought of:
-calculating the mean of the 5 datacolumns (y_2010 - y_2014) and then comparing all 5 datapoints with that mean. If there's a difference of e.g. >20%, this case is deleted and cannot be used for further analysis as there's too much variability in the data for that ID. (Normally the energy data over 5 years should be more or less the same unless retrofits have been done to the building, but most buildings haven't.)
- doing something with a rolling mean to come to a proper mean building energy use column
- ...
Example of the df:
   ID  y_2010   y_2011   y_2012  y_2013  y_2014  mean
21524   22631    21954    22314   22032   21843   ...
28965   27456    29654    28159   28654   27345   ...
10236   32165      NaN    31678   31895   32459   ...
89754   87621    86542    87542   88456   86961   ...
56457   58951    57486     2000       0       0   ...
25984   24587    25478      NaN   24896   25461   ...

 A: At the start, you have a fundamental decision to make:
Are you trying to learn something from the data? Or are you
trying to teach the data to behave as you suppose they should? 
This answer is mainly oriented toward the former approach.
Usually, it is a mistake to remove an 'outlier' from a dataset
unless you can establish that the observation in question arose from a documentable error (equipment failure, data entry error, etc.)
or you know for sure that its value is impossible (person's age above 140, negative height, etc.)

Anecdote: Where I live, the highest heating bill arises in December and January. For personal and family reasons unlikely ever to be
  repeated, my energy bills for 12/2019 and 1/2020 were extraordinarily high. By your criterion, I might be removed from your list, which I think would be a mistake. I can't foresee the same circumstances
  ever coming around for me again, but others may well be surprised by similar temporary periods of high energy usage in the future for very similar reasons.

One way to stabilize averages without tampering with data is to use trimmed means. To find a trimmed mean, data are sorted, a certain percentage of the very lowest and of the very highest
observations are ignored, and the mean of the more central remaining observations is taken. Depending on circumstances, typical trimming percentages may be 2% to 20% (sometimes higher), leaving the central 96% to 60% (sometimes fewer) to be averaged.
Consider data with $n=1000$ observations from
a gamma distribution with shape parameter 10 (perhaps waiting
times for finishing multi-phase projects.) Here is a graph of
its density curve--made in R.
curve(dgamma(x,10,1), 0, 25, lwd=2, ylab="PDF", 
      main="Density of GAMMA(10, 1)")
abline(v=0, col="green2");  abline(h=0, col="green2")


Almost all such samples have at least one boxplot outlier and the average number of outliers in a sample of 1000 is about 14.
set.seed(530)
nr.out = replicate(10^5, 
           length(boxplot.stats(rgamma(1000,10,1))$out) )
mean(nr.out);  mean(nr.out>0)
[1] 13.97049
[1] 1

Let's take a look at boxplots of 20 samples of size 1000 from this distribution in order to see the outliers.
set.seed(1234)
m = 20;  n=1000
x = rgamma(m*n,10,1);  g = rep(1:m, n)
boxplot(x~g, col="skyblue2", 
      main="GAMMA(10,1) Population: Boxplots of 20 Samples of 1000")

It seems that 2% trimmed means of the 1000 observations in each sample should allow us to ignore boxplot outliers in finding the means. (But the trimmed values are are not removed, so the ordinary mean
and quartiles are not affected.)

My simulated gamma observations have $\mu = 10, \sigma^2 = 10,$ so samples of
1000 have ordinary means average about $10$ with variances about $0.01$ (from theory). By contrast
 2% trimmed means of samples average about $9.93$ with variances about $0.01$ (from simulation).
set.seed(530)
a.02 = replicate(10^5, mean(rgamma(1000,10,1),trim=.02))
mean(a.02);  var(a.02)
[1] 9.932821
[1] 0.009988345

By using trimmed means we have retained all of the
data. In a fair and systematic way, we have mainly avoided using boxplot outliers to estimate means. Very roughly speaking we have
computed trimmed means by ignoring values that are more than double the
ordinary mean. And at the same time we have ignored values that are less than half the ordinary mean. Perhaps we find that
temporarily ignoring values that are proportionately far from
the ordinary mean (still the best estimate of the population mean), we can make better judgments from our data. 
qgamma(c(.02,.98), 10, 1)
[1]  4.618349 17.509813

However, over time we may come to realize that all of the observations have
a legitimate role to play in understanding how to use the data
to best advantage. In that case, the data are intact and we can do so.
Note: There are distributions with such heavy tails that
a trimmed sample mean is a better estimate of the location of the population than an ordinary sample mean.
The Cauchy is one such distribution. In that case tails are so heavy that
a 38% trimmed mean seems optimal. See a brief discussion 
here
and further information at its links.
