Nothing
### 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")
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.