I have built various different types of regression model (linear model, non-linear model, generalized linear model), and wish to determine the error/uncertainty of each one in order to compare them.
I have built the three models in R, and understand that I am able to use the predict
function to obtain a confidence interval around each of the fitted values. A bit of background about my models. Each output gives a prediction of a households electricity usage when it is given a predefined continuous variable. The way in which these models will be used, is to sum the predicted values, in order to obtain a total electricity consumption for an area. Therefore, i would like to quantify the error surrounding this summed value to determine which model has the least uncertainty in it. The models will then be used in other areas for which we do not have the actual values, so the model with the smallest error will be desired.
My understanding, is that I can do the following, to work out the model error:
- Compute the width of the confidence interval for each of the fitted values.
- Divide the width by two, for each fitted value.
- Each 'half-width' is equivalent to the error for that reading. Square each result.
- Sum all of the squared errors, and square root the result.
- This should give the error for the sum of the fitted values.
Firstly, is this correct? Does it give the error for the model?
Some R code (to demonstrate) is given at the end.
If this is correct, then my next question asks whether it is possible to use this same method for linear models which have manually computed coefficients?
Since I have created a linear model which uses the log of y and x, I have to take exponentials for my results to have meaning. Therefore, I have introduced Duan's Smearing factor, which applies a correction factor to the result. Therefore, my alpha (intercept) coefficient is different for one of my models, and i wish to determine the uncertainty of this type of model to compare.
R Code:
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%