JoSAE.R
In JoSAE: Unit-Level and Area-Level Small Area Estimation

### R code from vignette source 'JoSAE.Rnw'
### Encoding: ISO8859-1

###################################################
### code chunk number 1: JoSAE.Rnw:68-70 (eval = FALSE)
###################################################
## install.packages("JoSAE")
## 


###################################################
### code chunk number 2: JoSAE.Rnw:76-77
###################################################
library(JoSAE)


###################################################
### code chunk number 3: JoSAE.Rnw:81-82 (eval = FALSE)
###################################################
## ?JoSAE


###################################################
### code chunk number 4: JoSAE.Rnw:101-111
###################################################
	#mean auxiliary variables for the populations in the domains
data(JoSAE.domain.data)
	#data for the sampled elements
data(JoSAE.sample.data)

#plot(biomass.ha~mean.canopy.ht,JoSAE.sample.data)

library(lattice)

print(xyplot(biomass.ha ~ mean.canopy.ht | domain.ID, data = JoSAE.sample.data))


###################################################
### code chunk number 5: JoSAE.Rnw:132-135
###################################################
    #lme model
summary(fit.lme <- lme(biomass.ha ~ mean.canopy.ht, data=JoSAE.sample.data
                       , random=~1|domain.ID))


###################################################
### code chunk number 6: JoSAE.Rnw:144-147
###################################################
    #domain data need to have the same column names as sample data or vice versa
d.data <- JoSAE.domain.data
names(d.data)[3] <- "mean.canopy.ht"


###################################################
### code chunk number 7: JoSAE.Rnw:155-156
###################################################
result <- eblup.mse.f.wrap(domain.data = d.data, lme.obj = fit.lme)


###################################################
### code chunk number 8: JoSAE.Rnw:166-175
###################################################
library(xtable)
#tmp <- result[,c(1:2,6,4,9:11)]
tmp <- result[,c("domain.ID","N.i.domain","n.i.sample", "biomass.ha.sample.mean", "GREG"
                 , "EBLUP", "Synth")]
names(tmp)[4] <- c("sample.mean")
print(xtable(tmp, label = "tab:one", caption="Number of population and sampled elements as well as simple random sample, synthetic, GREG and EBLUP estimates of the mean above-ground forest biomass within 14 Norwegian municipalities."), include.rownames=FALSE)

tmp <- result[,c("domain.ID","n.i.sample","sample.se", "GREG.se", "EBLUP.se.1", "EBLUP.se.2")]#result[,c(1:2,6,23:26)]
print(xtable(tmp, label = "tab:two", caption="Number of population and sampled elements as well as standard errors of the simple random sample, GREG and EBLUP estimates of the mean above-ground forest biomass within 14 Norwegian municipalities."), include.rownames=FALSE)


###################################################
### code chunk number 9: JoSAE.Rnw:193-219
###################################################

tmp <- result[,c("biomass.ha.sample.mean", "Synth", "GREG", "EBLUP")]
    #actual plot
tmp1 <- barplot(t(as.matrix(tmp)), beside=T
            , names.arg=result$domain.ID
            , xlab="Municipalities"
                , ylab=expression(paste("Estimated biomass (Mg ", ha^{-1}, ")" ))
                , ylim=c(0,200))
    #print n.sample plots
text(tmp1[2,]+.5, y = 50, labels = result$n.i.sample,cex=1.5)
    #error bars
tmp2<- result[,c("sample.se", "sample.se", "GREG.se", "EBLUP.se.2")]#sample.se twice to fill the column, only used once.
tmp2[is.na(tmp2)] <- 0
        #plot error bars
            #sample mean
arrows(x0=tmp1[1,], y0=tmp[,1]+tmp2[,1], x1=tmp1[1,], y1 = tmp[,1]-tmp2[,1]
           , length = 0.01, angle = 90, code = 3)
            #GREG
arrows(x0=tmp1[3,], y0=tmp[,3]+tmp2[,3], x1=tmp1[3,], y1 = tmp[,3]-tmp2[,3]
           , length = 0.01, angle = 90, code = 3)
            #EBLUP
arrows(x0=tmp1[4,], y0=tmp[,4]+tmp2[,4], x1=tmp1[4,], y1 = tmp[,4]-tmp2[,4]
           , length = 0.01, angle = 90, code = 3)
    #legend
legend(13,200, fill=grey(c(.3, .6, .8, .9)), legend=c("SRS", "Synth", "GREG", "EBLUP"), bty="n")



###################################################
### code chunk number 10: JoSAE.Rnw:237-253
###################################################


data(landsat)

    #prepare the domain data - exclude "outlying" domain
landsat.domains <- unique(landsat[-33,c(1, 7:8,10)])
        #add a numeric domain ID
landsat.domains$domain.ID <- 1:nrow(landsat.domains)
        #change names to the names in the sample data
names(landsat.domains)[2:3] <- c("PixelsCorn", "PixelsSoybeans")

    #prepare the unit-level sample data
tmp <- landsat[-33,c(2:6, 10)]
        #add numeric domain ID
landsat.sample <- merge(landsat.domains[4:5], tmp, by="CountyName")



###################################################
### code chunk number 11: JoSAE.Rnw:259-266
###################################################
summary(landsat.lme <- lme(HACorn ~ PixelsCorn + PixelsSoybeans
                           , data=landsat.sample
                           , random=~1|domain.ID))
    #obtain EBLUP estimates and MSE
result <- eblup.mse.f.wrap(domain.data = landsat.domains
             , lme.obj = landsat.lme)



###################################################
### code chunk number 12: JoSAE.Rnw:269-279
###################################################

tmp <- result[,c("CountyName.domain","n.i.sample","EBLUP", "EBLUP.se.1", "EBLUP.se.2", "GREG.se")]
#tmp <- result[,c(5,9,13,28,29,27)]
names(tmp)[1:2] <- c("County.name","n_i")#, "EBLUP.mean", "EBLUP.se.1", "EBLUP.se.2", "GREG.se")
#names(tmp) <- c("County.name","n_i", "EBLUP.mean", "EBLUP.se.1", "EBLUP.se.2", "GREG.se")
print(xtable(tmp, label = "tab:landsat"
             , caption="EBLUP estimates of county means of hectares under corn and estimated standard errors of the EBLUP and GREG estimates.")
      , include.rownames=FALSE)




###################################################
### code chunk number 13: JoSAE.Rnw:308-316
###################################################
data(nfi)
    #fit the model
fit.nfi.iw <- lm(vol.2011~Elev.Mean, nfi, weights=1/Elev.Mean)
  #data (model matrix, X) of one validation stand
stand <- cbind(Intercept=1, Elev.Mean=c(147.41,127.48,98.66,118.85,124,120.81,119.7),
               N=c(0,23,0,55,27,80,56), E=c(73,77,0,39,37,54,54) )
  #aggregate to obtain X-bar
stand.agg <- apply(stand[,1:2], 2, mean)


###################################################
### code chunk number 14: JoSAE.Rnw:325-330
###################################################
    #obtain covariance matrix 
Sigma <- vcov(fit.nfi.iw)

    #residual variance
sig <- summary(fit.nfi.iw)$sigma


###################################################
### code chunk number 15: JoSAE.Rnw:335-336
###################################################
var.p <- t(stand.agg) %*% Sigma %*% stand.agg


###################################################
### code chunk number 16: JoSAE.Rnw:342-343
###################################################
var.prh <- var.p + sum(sig^2 * stand[,2])/nrow(stand)^2


###################################################
### code chunk number 17: JoSAE.Rnw:348-350
###################################################
library(nlme)
spG <- corGaus(23, form = ~N+E)


###################################################
### code chunk number 18: JoSAE.Rnw:356-359
###################################################
cormat <- corMatrix(Initialize(spG, data.frame(stand)))

varmat <- (sig * sqrt(stand[,2])) %o% (sig * sqrt(stand[,2]))


###################################################
### code chunk number 19: JoSAE.Rnw:364-365
###################################################
var.prhs <- var.p + sum(cormat * varmat)/nrow(stand)^2


###################################################
### code chunk number 20: JoSAE.Rnw:371-372
###################################################
sqrt(c(var.p, var.prh, var.prhs))