# Generate data for the running example
# X: Gaussian but with more-than-standard variance
x = rnorm(100,0,3)
# Y: Linearly dependent on X but with variance that grows with X
# (heteroskedastic)
y = 3-2*x + rnorm(100,0,sapply(x,function(x){1+0.5*x^2}))
# Plot the data
plot(x,y)
# Plot the true regression line
abline(a=3,b=-2,col="grey")
# Fit by ordinary least squares
fit.ols = lm(y~x)
# Plot that line
abline(fit.ols$coefficients,lty=2)
#### Save first figure at this point
# Fit a weighted-least-squares regression (with weights inversely proportional
# to actual variances), add that line
fit.wls = lm(y~x, weights=1/(1+0.5*x^2))
abline(fit.wls$coefficients,lty=3)
##### Save figure again
# Generate more random samples from the same model and the same x values,
# but different y values
# Inputs: number of samples to generate
# Presumes: x exists and is defined outside this function
# Outputs: errors in linear regression estimates
ols.heterosked.example = function(n) {
y = 3-2*x + rnorm(n,0,sapply(x,function(x){1+0.5*x^2}))
fit.ols = lm(y~x)
# Return the errors
return(fit.ols$coefficients - c(3,-2))
}
# Calculate average-case errors in linear regression estimates (SD of
# slope and intercept)
# Inputs: number of samples per replication, number of replications (defaults
# to 10,000)
# Calls: ols.heterosked.example
# Outputs: standard deviation of intercept and slope
ols.heterosked.error.stats = function(n,m=10000) {
ols.errors.raw = t(replicate(m,ols.heterosked.example(n)))
# transpose gives us a matrix with named columns
intercept.sd = sd(ols.errors.raw[,"(Intercept)"])
slope.sd = sd(ols.errors.raw[,"x"])
return(list(intercept.sd=intercept.sd,slope.sd=slope.sd))
}
### As previous two functions, but with weighted regression
# Generate random sample from model (with fixed x), fit by weighted least
# squares
# Inputs: number of samples
# Presumes: x fixed outside function
# Outputs: errors in parameter estimates
wls.heterosked.example = function(n) {
y = 3-2*x + rnorm(n,0,sapply(x,function(x){1+0.5*x^2}))
fit.wls = lm(y~x,weights=1/(1+0.5*x^2))
# Return the errors
return(fit.wls$coefficients - c(3,-2))
}
# Calculate standard errors in parameter estiamtes over many replications
# Inputs: number of samples per replication, number of replications (defaults
# to 10,000)
# Calls: wls.heterosked.example
# Outputs: standard deviation of estimated intercept and slope
wls.heterosked.error.stats = function(n,m=10000) {
wls.errors.raw = t(replicate(m,wls.heterosked.example(n)))
# transpose gives us a matrix with named columns
intercept.sd = sd(wls.errors.raw[,"(Intercept)"])
slope.sd = sd(wls.errors.raw[,"x"])
return(list(intercept.sd=intercept.sd,slope.sd=slope.sd))
}