## ols2.r
## PJ Feb. 8, 2002
## This is a little fancier version of ols1. This one
## highlights the fact that the standard error of b estimated
## in an individual run are close to the true standard deviation
## of b in its sampling distribution.
nOfRuns <- 100
a <- 2 #intercept parameter
b <- 5 #slope parameter
stderror <- 3 #standard deviation of error term
getPhonyData <- function (x){
e<-rnorm(1000,sd=stderror)
y<-b + b*x + e
}
x<- rnorm(1000,mean=10,sd=7)
conductSim <- function(i){
y<-getPhonyData(x) #x is grabbed from environment
try(lm(y~x)) #output of last command is returned
}
## use lapply to run the sim over and over, stuff
## estimate objects into a list called result
result <- lapply(1:nOfRuns, conductSim)
## Get the estimates from run i
## Note you get the i'th one by result[[i]].
## here I take all coefficients
getCoeff <- function(i){result[[i]]$coefficients}
## collect up parameter estimates from all runs
paramvector <- sapply(1:nOfRuns,getCoeff)
## get standard errors of estimates of "a" and "b" from run i
## second column of coefficients in summary is the standard errors
## so grab them. There must be an easier way!
getSE <- function(i){summary(result[[i]])$coefficients[,2] }
## collect up standard errors of "a" and "b" from runs
sevector <- sapply(1:nOfRuns,getSE)
## print those out for fun
sevector
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## (Intercept) 0.1699 0.16959 0.16426 0.16515 0.16363 0.16486 0.16724 0.17167
## x 0.0138 0.01378 0.01334 0.01341 0.01329 0.01339 0.01358 0.01394
## [,9] [,10] [,11] [,12] [,13] [,14] [,15]
## (Intercept) 0.17628 0.16695 0.16426 0.17255 0.16657 0.16690 0.17518
## x 0.01432 0.01356 0.01334 0.01402 0.01353 0.01356 0.01423
## [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23]
## (Intercept) 0.16380 0.16312 0.1723 0.1637 0.16941 0.16643 0.16847 0.16779
## x 0.01331 0.01325 0.0140 0.0133 0.01376 0.01352 0.01368 0.01363
## [,24] [,25] [,26] [,27] [,28] [,29] [,30] [,31]
## (Intercept) 0.17079 0.16907 0.17002 0.16463 0.16549 0.17070 0.1662 0.17069
## x 0.01387 0.01373 0.01381 0.01337 0.01344 0.01387 0.0135 0.01386
## [,32] [,33] [,34] [,35] [,36] [,37] [,38]
## (Intercept) 0.16724 0.16170 0.16468 0.17047 0.16965 0.16813 0.16968
## x 0.01359 0.01314 0.01338 0.01385 0.01378 0.01366 0.01378
## [,39] [,40] [,41] [,42] [,43] [,44] [,45] [,46]
## (Intercept) 0.16474 0.16279 0.1674 0.16283 0.17147 0.17057 0.16914 0.16363
## x 0.01338 0.01322 0.0136 0.01323 0.01393 0.01386 0.01374 0.01329
## [,47] [,48] [,49] [,50] [,51] [,52] [,53] [,54]
## (Intercept) 0.16361 0.16431 0.16389 0.16446 0.1638 0.17051 0.17460 0.16207
## x 0.01329 0.01335 0.01331 0.01336 0.0133 0.01385 0.01418 0.01317
## [,55] [,56] [,57] [,58] [,59] [,60] [,61]
## (Intercept) 0.16487 0.16565 0.16669 0.17049 0.16639 0.16826 0.16539
## x 0.01339 0.01346 0.01354 0.01385 0.01352 0.01367 0.01343
## [,62] [,63] [,64] [,65] [,66] [,67] [,68] [,69]
## (Intercept) 0.16824 0.16674 0.16357 0.1674 0.16546 0.16760 0.15895 0.16715
## x 0.01367 0.01354 0.01329 0.0136 0.01344 0.01361 0.01291 0.01358
## [,70] [,71] [,72] [,73] [,74] [,75] [,76] [,77]
## (Intercept) 0.17099 0.16764 0.16705 0.16763 0.16913 0.16404 0.17705 0.1711
## x 0.01389 0.01362 0.01357 0.01362 0.01374 0.01332 0.01438 0.0139
## [,78] [,79] [,80] [,81] [,82] [,83] [,84] [,85]
## (Intercept) 0.16170 0.16773 0.16686 0.16779 0.16857 0.1650 0.17054 0.16722
## x 0.01313 0.01362 0.01355 0.01363 0.01369 0.0134 0.01385 0.01358
## [,86] [,87] [,88] [,89] [,90] [,91] [,92] [,93]
## (Intercept) 0.16964 0.16626 0.17274 0.16556 0.16297 0.1711 0.16604 0.16562
## x 0.01378 0.01351 0.01403 0.01345 0.01324 0.0139 0.01349 0.01345
## [,94] [,95] [,96] [,97] [,98] [,99] [,100]
## (Intercept) 0.16963 0.17329 0.16067 0.17268 0.17248 0.16883 0.16592
## x 0.01378 0.01408 0.01305 0.01403 0.01401 0.01371 0.01348
## the average of estimates of b across runs
meanb <- mean(paramvector[2,])
meanb
## [1] 4.999
## the standard deviation of estimates of b across runs
stdb <- sd(paramvector[2,])
stdb
## [1] 0.01291
hist(paramvector[2,])
text(x = 4.96, y = 20, paste("The true value of b is ", b))
text(x = 4.96, y = 15, paste("The average of the estimates \n of b is ", meanb))
## average of estimated standard errors of b
meanseb <- mean(sevector[2,])
## new histogram!
hist(sevector[2,], main="Sampling Distribution of Estimates of the Standard Error of b")
text(x = 0.0128, y = 24, adj = 0, paste("Std.Dev. of sampling dist of b= ", stdb))
text(x = 0.0128, y = 20, adj = 0, paste("The of mean se(b) was ", meanseb))
text(x = 0.0128, y = 15, adj = 0, paste ("The true sd(b) is ", stderror/sqrt(sum((x-mean(x))*(x-mean(x))))))
