rsae.R
In rsae: Robust Small Area Estimation

### R code from vignette source 'rsae.Rnw'
### Encoding: UTF-8

###################################################
### code chunk number 1: rsae.Rnw:172-173 (eval = FALSE)
###################################################
## install.package("rsae")


###################################################
### code chunk number 2: rsae.Rnw:181-182
###################################################
library(rsae)


###################################################
### code chunk number 3: rsae.Rnw:204-205
###################################################
data(landsat)


###################################################
### code chunk number 4: rsae.Rnw:209-210
###################################################
names(landsat)


###################################################
### code chunk number 5: rsae.Rnw:214-215
###################################################
table(landsat$CountyName)


###################################################
### code chunk number 6: rsae.Rnw:222-223
###################################################
landsat[30:35, c("HACorn", "PixelsCorn", "PixelsSoybeans", "CountyName", "outlier")]


###################################################
### code chunk number 7: rsae.Rnw:228-230
###################################################
   linmodel <- lm(HACorn ~ PixelsCorn + PixelsSoybeans, data = landsat)
   plot(linmodel, 5)


###################################################
### code chunk number 8: rsae.Rnw:236-238 (eval = FALSE)
###################################################
##    linmodel <- lm(HACorn ~ PixelsCorn + PixelsSoybeans, data = landsat)
##    plot(linmodel, 5)


###################################################
### code chunk number 9: rsae.Rnw:257-260
###################################################
bhfmodel <- saemodel(formula=HACorn ~ PixelsCorn + PixelsSoybeans, 
                    area=~CountyName, 
                    data=subset(landsat, subset=(outlier==FALSE)))


###################################################
### code chunk number 10: rsae.Rnw:267-268
###################################################
bhfmodel


###################################################
### code chunk number 11: rsae.Rnw:273-274
###################################################
summary(bhfmodel)


###################################################
### code chunk number 12: rsae.Rnw:294-295
###################################################
mlfit <- fitsaemodel("ml", bhfmodel)


###################################################
### code chunk number 13: rsae.Rnw:299-300
###################################################
mlfit


###################################################
### code chunk number 14: rsae.Rnw:309-310
###################################################
summary(mlfit)


###################################################
### code chunk number 15: rsae.Rnw:319-320
###################################################
failconvg <- fitsaemodel("ml", bhfmodel, niter=1)


###################################################
### code chunk number 16: rsae.Rnw:325-326
###################################################
failconvg


###################################################
### code chunk number 17: rsae.Rnw:331-332
###################################################
convergence(failconvg)


###################################################
### code chunk number 18: rsae.Rnw:347-354
###################################################
set.seed(12345)
n <- 200; beta <- c(1, 1)
cst <- rep(1, n)
x <- rnorm(n)
y <- as.matrix(cbind(cst, x)) %*% beta + rnorm(n)
areaid <- rep(1:10, each=10)
df <- data.frame(y=y, x=x, areaid=areaid)


###################################################
### code chunk number 19: rsae.Rnw:358-359
###################################################
lm(y ~ x, data=df)


###################################################
### code chunk number 20: rsae.Rnw:363-364
###################################################
fakemlm <- saemodel(y ~ x, area=~areaid, data=df)


###################################################
### code chunk number 21: rsae.Rnw:369-370
###################################################
fitsaemodel("ml", fakemlm)


###################################################
### code chunk number 22: rsae.Rnw:390-391
###################################################
huberfit <- fitsaemodel("huberm", bhfmodel, k=1.5)


###################################################
### code chunk number 23: rsae.Rnw:398-399
###################################################
huberfit


###################################################
### code chunk number 24: rsae.Rnw:404-405
###################################################
summary(huberfit)


###################################################
### code chunk number 25: rsae.Rnw:419-420
###################################################
convergence(huberfit)


###################################################
### code chunk number 26: rsae.Rnw:457-458
###################################################
fitsaemodel("huberm", bhfmodel, k=1.2, init="s")


###################################################
### code chunk number 27: rsae.Rnw:489-493
###################################################
d <- unique(landsat[-33, c("MeanPixelsCorn", "MeanPixelsSoybeans", "CountyName")])
d <- cbind(rep(1,12), d)
rownames(d) <- d$CountyName
d <- d[,1:3]


###################################################
### code chunk number 28: rsae.Rnw:497-498
###################################################
d


###################################################
### code chunk number 29: rsae.Rnw:505-506
###################################################
pr <- robpredict(mlfit, areameans=d, reps=500)


###################################################
### code chunk number 30: rsae.Rnw:510-511
###################################################
pr 


###################################################
### code chunk number 31: rsae.Rnw:516-517
###################################################
plot(pr)


###################################################
### code chunk number 32: rsae.Rnw:522-523
###################################################
plot(pr, type="l", sort="means")


###################################################
### code chunk number 33: rsae.Rnw:537-540
###################################################
res <- residuals(pr)
qqnorm(res)
qqline(res)


###################################################
### code chunk number 34: rsae.Rnw:592-593
###################################################
mymodel <- makedata()


###################################################
### code chunk number 35: rsae.Rnw:600-601
###################################################
mymodel


###################################################
### code chunk number 36: rsae.Rnw:606-607
###################################################
summary(mymodel)


###################################################
### code chunk number 37: rsae.Rnw:615-616
###################################################
makedata(ve.epsilon=0.1)


###################################################
### code chunk number 38: rsae.Rnw:621-622
###################################################
makedata(vu.epsilon=0.1)


###################################################
### code chunk number 39: rsae.Rnw:627-628
###################################################
makedata(vu.epsilon=0.1, ve.epsilon=0.1)


###################################################
### code chunk number 40: rsae.Rnw:638-639
###################################################
my_area_id = c(1,1,1,1,1,2,2,3,3,3,4,4,4,4,4,5)


###################################################
### code chunk number 41: rsae.Rnw:646-647
###################################################
mymodelub <- makedata(n=NULL, g=NULL, areaID=my_area_id)


###################################################
### code chunk number 42: rsae.Rnw:658-659
###################################################
citation("rsae")


###################################################
### code chunk number 43: rsae.Rnw:667-668
###################################################
fitsaemodel("ml", bhfmodel)


###################################################
### code chunk number 44: rsae.Rnw:673-675
###################################################
require(nlme)
nlme::lme(HACorn ~ PixelsCorn + PixelsSoybeans, random=~1|CountyName, data=subset(landsat, subset=(outlier==FALSE)), method="ML")