I measure the response of the ADC for over 4000 voltages. I use a linear fit to these measurements. This gives very good results for almost all ATmega chips.
The biggest error is usually around change over in the MSB. Do you see evidence of that in the 4000 readings?
If a simple linear correction is to be made , it may be better to do a separate fit for each half of the range, or the same slope would suffice but with a different offset for each half of the scale.
WRT 1023 vs 1024, I realised after thinking about the succ. approx technique that 2^N is the correct factor not 2^N-1 . The approx is done by halving the interval to find each bit, so it is based on full scale of Vref. The implication is that the ADC will never measure Vref , only Vref-1LSB as its max value .
This ties in to your +0.5 suggestion. You are correct that would be statistically more accurate.
Considering a 1 bit S-A conversion , it will only return 0 or 1 ie 0% or 50% Vref. Clearly 25% or 75% would be a more balanced representation of the whole input range.