H Kelsey,
Here is one way:
# some data
set.seed(123)
x <- rnorm(100)
e <- rexp(100, 2)
y <- 3 + 1.5*x - 1.3*x^2 + e
plot(x, y, las = 1)
d <- data.frame(x, y)
d
# model
fit <- lm(y ~ x + I(x^2) , data = d)
summary(fit)
rse <- summary(fit)$sigma
rse
# [1] 0.4646164
# function to calculate the RSE
rse.boot <- function(d){
di <- d[sample(NROW(d), replace = TRUE), ]
fit <- lm(y ~ x + I(x^2) , data = di)
rse <- summary(fit)$sigma
rse
}
# 5000 bootstrap samples for the RSE
out <- replicate(5000, rse.boot(d))
head(out)
hist(out, prob = TRUE, las = 1, main = 'RSE distribution', xlab = "RSE*")
# confidence interval
c(RSE = rse, quantile(out, probs = c(0.025, 0.975)))
#RSE 2.5% 97.5%
#0.4646164 0.3700671 0.5367165
HTH,
Jorge.-
On Wed, May 9, 2012 at 9:15 PM, missToo Quick <> wrote:
> Hello,
> I am trying to write a code to estimate the residual standard error and
> create a confidence interval using bootstrap, since it does not follow
> a normal distribution.
>
> So far I have found a linear model for the data (m1<-lm(y~x+I(x^2))), but
> I am not sure how to create the bootstrap code to estimate its residuals.
>
> Thank you for any input,
> Kelsey
> [[alternative HTML version deleted]]
>
>
> ______________________________________________
> [email protected] mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
>
[[alternative HTML version deleted]]
______________________________________________
[email protected] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
