saipe_eb

Reading in the data

### read the saipe state-level data for the year of interest

rm(list = ls(all = TRUE))

saipedat <- read.csv('data93.csv',header=TRUE)
ve<-saipedat$fnlse^2

### extract the variables of interest

alldat <- saipedat[,c('cps','irspr','irsnf','fs','fnlse','smpsize','cenres')]
alldat$ve=ve
###

Creat plots

par(mfrow=c(2,2))
plot(alldat$cps~alldat$irspr,ylab='1993 CPS poverty rates', xlab='IRS pseudo poverty rates', main='(a)',cex.main=0.8,cex.lab=0.8)

plot(alldat$cps~alldat$irsnf,ylab='1993 CPS poverty rates', xlab='IRS non-filer rates', main='(b)',cex.main=0.8,cex.lab=0.8)

plot(alldat$cps~alldat$fs,ylab='1993 CPS poverty rates', xlab='Food stamp participation rates', main='(c)',cex.main=0.8,cex.lab=0.8)

plot(alldat$cps~alldat$cenres,ylab='1993 CPS poverty rates', xlab='Census residuals', main='(d)',cex.main=0.8,cex.lab=0.8)

model

### Fay-Herriot models

library('sae')

## Warning: package 'sae' was built under R version 3.6.3

## Loading required package: MASS

## Loading required package: lme4

## Warning: package 'lme4' was built under R version 3.6.3

## Loading required package: Matrix

## one variable

m1 <- mseFH(alldat$cps ~ alldat$irspr, (alldat$ve),'REML')

m2 <- mseFH(alldat$cps ~ alldat$irsnf, alldat$ve,'REML')

m3 <- mseFH(alldat$cps ~ alldat$fs, (alldat$ve),'REML')

m4 <- mseFH(alldat$cps ~ alldat$cenres, (alldat$ve),'REML')

## two variables

m12 <-mseFH(alldat$cps ~ alldat$irsnf+ alldat$irspr, (alldat$ve),'REML')
m13 <-mseFH(alldat$cps ~ alldat$irspr + alldat$fs, (alldat$ve),'REML')
m14 <-mseFH(alldat$cps ~ alldat$irspr + alldat$cenres, (alldat$ve),'REML')
m23 <-mseFH(alldat$cps ~ alldat$irsnf + alldat$fs, (alldat$ve),'REML')
m24 <-mseFH(alldat$cps ~ alldat$irsnf + alldat$cenres, (alldat$ve),'REML')
m34 <-mseFH(alldat$cps ~ alldat$fs +alldat$cenres, (alldat$ve),'REML')

## three variables

m123 <-mseFH(alldat$cps ~ alldat$irsnf+ alldat$fs+alldat$irspr, (alldat$ve),'REML')
m124 <-mseFH(alldat$cps ~ alldat$irsnf+alldat$irspr +alldat$cenres, (alldat$ve),'REML')
m134 <-mseFH(alldat$cps ~ alldat$cenres+ alldat$fs+alldat$irspr, (alldat$ve),'REML')
m234 <-mseFH(alldat$cps ~ alldat$irsnf+ alldat$fs+alldat$cenres, (alldat$ve),'REML')


## four variables

m1234 <-mseFH(alldat$cps ~ alldat$irspr +alldat$irsnf+ alldat$fs+alldat$cenres, (alldat$ve),'REML')

#m1234$est$fit

Output

### output

d1 <- round(m1234$est$fit$estcoef,3)

out <- data.frame(modelvar=c(m1$est$fit$refvar , m2$est$fit$refvar ,m3$est$fit$refvar ,m4$est$fit$refvar ,
                             m12$est$fit$refvar ,m13$est$fit$refvar ,m14$est$fit$refvar ,m23$est$fit$refvar,m24$est$fit$refvar ,m34$est$fit$refvar ,
                             m123$est$fit$refvar,m124$est$fit$refvar,m134$est$fit$refvar,m234$est$fit$refvar,m1234$est$fit$refvar ),
                  AICs = c( m1$est$fit$goodness['AIC'], m2$est$fit$goodness['AIC'], m3$est$fit$goodness['AIC'],m4$est$fit$goodness['AIC'],
                            m12$est$fit$goodness['AIC'], m13$est$fit$goodness['AIC'], m14$est$fit$goodness['AIC'],m23$est$fit$goodness['AIC'],m24$est$fit$goodness['AIC'],m34$est$fit$goodness['AIC'],
                            m123$est$fit$goodness['AIC'],m124$est$fit$goodness['AIC'],m134$est$fit$goodness['AIC'],m234$est$fit$goodness['AIC'],m1234$est$fit$goodness['AIC']))

# library('xtable')
# xtable(d1)
# xtable(out)

Results

#### results

### model estimates versus direct estimates

par(mfrow=c(1,1))

plot(alldat$cps~m1234$est$eblup,xlab='1993 CPS poverty rate model estimates', ylab='1993 CPS poverty rate direct estimates', main='',cex.main=0.8,cex.lab=0.8)
abline(0,1)

### standardized residuals for full model

par(mfrow=c(1,3))

r1234 <- (alldat$cps-m1234$est$fit$estcoef[1,1]-m1234$est$fit$estcoef[2,1]*alldat$irspr - m1234$est$fit$estcoef[3,1]*alldat$irsnf- m1234$est$fit$estcoef[4,1]*alldat$fs-m1234$est$fit$estcoef[5,1]*alldat$cenres)/sqrt(alldat[,'ve']+m1234$est$fit$refvar)

plot(r1234~m1234$est$eblup,xlab='1993 CPS poverty rate model estimates', ylab='Standardized residuals', main='(a)',cex.main=0.8,cex.lab=0.8)
abline(h=0)

# fittedvals <- m1234$est$fit$estcoef[1,1]+m1234$est$fit$estcoef[2,1]*alldat$irsprchld + m1234$est$fit$estcoef[3,1]*alldat$irsnf0.64+ m1234$est$fit$estcoef[4,1]*alldat$fsrate+m1234$est$fit$estcoef[5,1]*alldat$cenrsd
# plot(r1234~fittedvals,xlab='1993 CPS poverty rate model fitted values', ylab='Standardized residuals', main='(a)',cex.main=0.8,cex.lab=0.8)
# abline(h=0)

qqnorm(r1234,main='(b)',cex.main=0.8,cex.lab=0.8)
abline(0,1)

hist(r1234, main='(c)',xlab='',cex.main=0.8,cex.lab=0.8)

### estimated CVs

par(mfrow=c(1,2))

cvdir <- sqrt(alldat$ve)/alldat$cps
cveblup <- sqrt(m1234$mse)/m1234$est$eblup

plot((sqrt(m1234$mse)/sqrt(alldat$ve))~alldat$smpsize, xlab='CPS sample size (log scale)',ylab='SE model estimates over SE direct estimates',main='(a)',cex.main=0.8,cex.lab=0.8,log='x')
plot( (cveblup/cvdir) ~ alldat$smpsize, xlab='CPS sample size (log scale)',ylab='CV model estimates over CV direct estimates', main='(b)',cex.main=0.8,cex.lab=0.8  ,log='x')