DF=data.frame(y=c(182.00455,606.27273,2401.03918,1597.27500,179.68846,76.82728,313.85000,1431.00000,16.43620,22887.81818,1010.00000,20.65909,184.17273,483.35000,21.45000,291.01500,359.10000,602.75000,253.18636,35.74286),
              x=c(133.955464,1.913142,88.887131,95.512793,25.247257,11.938203,51.246909,96.265030,42.701863,9.082072,42.466148,86.741979,15.908011,55.756779,79.432585,61.395584,22.822762,22.853197,30.154734,96.494249))

fit <- lm(log(y) ~ log(x), data=DF)
summary(fit)

prediction=predict.lm(fit, interval="confidence")

prediction=data.frame(prediction)
prediction=exp(prediction) # Since the model was built using log, the values are exponentiated before use.
prediction$halfwidth=(prediction$upr-prediction$lwr)/2
prediction$error_squared=prediction$halfwidth^2

sqrt(sum(prediction$error_squared, na.rm=T))

sum(prediction$fit)

## uncertainty percentage:
(sqrt(sum(prediction$error_squared, na.rm=T))/sum(prediction$fit))*100
# [1] 108.656%