``````## 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))))))
``````