Removing outliers is common practice in statistical modeling and perfectly acceptable. However, with regards 1.5 IQR I am far from certain about this approach. Normally, if you want to be conservative then 3 standard deviations (SD) denote an outlier, which is more stringent than IQR. Some use 2 SD. If the value is lower than 2 SD from the group mean it isn't an outlier, that's because 1.96 SD is the normal distribution.
Normally you would delete the outlier and it takes no further part in the calculation, assigning it to 'na' should work. The subsequent calculation will recognise an 'na' value, I think in python it is np('na') ... Np = numpy.
In this case you appear to be deleting the value and treating it as missing data. Treatment of missing data is complicated I agree 100%. Mean, median are possibilities for missing data I can't see this working for binary data. One approach is regression and in python there is a command specifically to do this (which I forget), R will have the same command and it will assign a value for the missing data. However, you can't do this for outliers.
The other routes to remove outliers are via transformation, log transformation is the most common. Transformation is a big subject.