.Rhistory

}
#define a seriesid(s) based on user input
seriesids <- paste0("LAUCN"
,fips_sub$State.ANSI
,sprintf("%03d",fips_sub$County.ANSI)
,"00000000"
,measure_sub)
#create a JSON file of all parameters specified
#this registrationKey is limited to 500 queries per day.
jsonbody <- toJSON(list('seriesid'=seriesids,
'startyear'=startyear,
'endyear'=endyear,
'registrationKey' = APIkey))
#the BLS API will fail if a single series is not in brackets.
#add brackets if we only have one seriesid
if (length(seriesids) == 1) {
jsonbody <- gsub(x=jsonbody
,replacement=paste0('[\"', seriesids, '\"]')
,paste0('\"', seriesids, '\"'))
}
#set the header for the http request
httpheader <- c(Accept="application/json; charset=UTF-8",
"Content-Type"="application/json")
#get BLS data with postForm method
bls_content <- postForm("http://api.bls.gov/publicAPI/v2/timeseries/data/"
,.opts=list(httpheader=httpheader
,postfields=jsonbody))
#convert the returned data to a List
bls_content <- fromJSON(bls_content[1])
#check count of series returned
series_cnt <- length(bls_content$Results$series)
#check count of data within the series (should be the same for every series)
data_cnt <- length(bls_content$Results$series[[1]]$data)
#initialize an object for storing data from the JSON object
results <- NULL
#loop through each series and data element
for (i in 1:series_cnt){
for (j in 1:data_cnt) {
nextrow <- c(bls_content$Results$series[[i]]$seriesID
,bls_content$Results$series[[i]]$data[[j]]$year
,bls_content$Results$series[[i]]$data[[j]]$period
,bls_content$Results$series[[i]]$data[[j]]$periodName
,bls_content$Results$series[[i]]$data[[j]]$value)
results <- rbind(results, nextrow)
}
}
#clean data and add some useful fields of information
rownames(results) <- NULL
results <- as.data.frame(results)
colnames(results) <- c("seriesID", "year", "period", "periodName", "value")
results$year <- as.numeric(as.character(results$year))
results$period <- as.numeric(substring(as.character(results$period), 2))
results$value <- as.numeric(as.character(results$value))
results$cd_measure <- substring(as.character(results$seriesID), 19, 20)
results$cd_fips_state <- as.numeric(substring(as.character(results$seriesID), 6, 7))
results$cd_fips_county <- as.numeric(substring(as.character(results$seriesID), 8, 10))
results$month <- as.Date(paste0(results$year, "-", results$period, "-01"))
require(sqldf)
results <- sqldf("select distinct
r.*
,rlf.State tx_state
,rlf.County_Name tx_county
,rlm.tx_measure
from results r
join ref_lookup_fips rlf
on rlf.State_ANSI = r.cd_fips_state
and rlf.County_ANSI = r.cd_fips_county
join ref_lookup_measure rlm
on rlm.cd_measure = r.cd_measure")
return(results)
}
dat <- getLAUCN(states = "WA", counties = "*", measures = "03", startyear = "2000", endyear = "2013")
head(dat)
getLAUCN <- function(states = c("*")
,counties = c("*")
,measures = c("03")
,startyear = NULL
,endyear = NULL
,APIkey = "07d259de03a845fea8736371980a19ad")
{
require(RJSONIO)
require(sqldf)
require(RCurl)
#Load fips codes for states and counties from the US Census
ref_lookup_fips <- read.csv("http://www.census.gov/geo/reference/codes/files/national_county.txt")
#Drop the County Suffix on the text description
ref_lookup_fips$County.Name <- gsub(" County", "", ref_lookup_fips$County.Name)
#Define a table of measure codes (based on http://www.bls.gov/help/def/la.htm#measure)
ref_lookup_measure <- data.frame(cd_measure = c("06", "05", "04", "03")
,tx_measure = c("Labor force",
"Employment",
"Unemployment",
"Unemployment rate")
)
#define a vector of county and state indices based on user specification
county_idx <- which(ref_lookup_fips$County.Name %in% counties)
state_idx <- which(ref_lookup_fips$State %in% states)
#check to see if user wants all counties and states
if (any(states == "*") & any(counties == "*")){
fips_sub <- ref_lookup_fips
#check to see if user wants all states limited by counties
} else if (any(states == "*") & all(counties != "*")){
fips_sub <- ref_lookup_fips[county_idx,]
#check to see if user wants all counties limited by states
} else if (all(states != "*") & any(counties == "*")){
fips_sub <- ref_lookup_fips[state_idx,]
} else{
#otherwise, take the intersection of states and counties
fips_sub <- ref_lookup_fips[intersect(county_idx, state_idx),]
}
#check to see if we need to expand the measurement list to match the
#fips ids
if (nrow(fips_sub) > 1) {
measure_sub <- rep(measures, nrow(fips_sub))
}
#define a seriesid(s) based on user input
seriesids <- paste0("LAUCN"
,fips_sub$State.ANSI
,sprintf("%03d",fips_sub$County.ANSI)
,"00000000"
,measure_sub)
#create a JSON file of all parameters specified
#this registrationKey is limited to 500 queries per day.
jsonbody <- toJSON(list('seriesid'=seriesids,
'startyear'=startyear,
'endyear'=endyear,
'registrationKey' = APIkey))
#the BLS API will fail if a single series is not in brackets.
#add brackets if we only have one seriesid
if (length(seriesids) == 1) {
jsonbody <- gsub(x=jsonbody
,replacement=paste0('[\"', seriesids, '\"]')
,paste0('\"', seriesids, '\"'))
}
#set the header for the http request
httpheader <- c(Accept="application/json; charset=UTF-8",
"Content-Type"="application/json")
#get BLS data with postForm method
bls_content <- postForm("http://api.bls.gov/publicAPI/v2/timeseries/data/"
,.opts=list(httpheader=httpheader
,postfields=jsonbody))
#convert the returned data to a List
bls_content <- fromJSON(bls_content[1])
#check count of series returned
series_cnt <- length(bls_content$Results$series)
#check count of data within the series (should be the same for every series)
data_cnt <- length(bls_content$Results$series[[1]]$data)
#initialize an object for storing data from the JSON object
results <- NULL
#loop through each series and data element
for (i in 1:series_cnt){
for (j in 1:data_cnt) {
nextrow <- c(bls_content$Results$series[[i]]$seriesID
,bls_content$Results$series[[i]]$data[[j]]$year
,bls_content$Results$series[[i]]$data[[j]]$period
,bls_content$Results$series[[i]]$data[[j]]$periodName
,bls_content$Results$series[[i]]$data[[j]]$value)
results <- rbind(results, nextrow)
}
}
#clean data and add some useful fields of information
rownames(results) <- NULL
results <- as.data.frame(results)
colnames(results) <- c("seriesID", "year", "period", "periodName", "value")
results$year <- as.numeric(as.character(results$year))
results$period <- as.numeric(substring(as.character(results$period), 2))
results$value <- as.numeric(as.character(results$value))
results$cd_measure <- substring(as.character(results$seriesID), 19, 20)
results$cd_fips_state <- as.numeric(substring(as.character(results$seriesID), 6, 7))
results$cd_fips_county <- as.numeric(substring(as.character(results$seriesID), 8, 10))
results$month <- as.Date(paste0(results$year, "-", results$period, "-01"))
require(sqldf)
results <- sqldf("select distinct
r.*
,rlf.State tx_state
,rlf.County_Name tx_county
,rlm.tx_measure
from results r
join ref_lookup_fips rlf
on rlf.State_ANSI = r.cd_fips_state
and rlf.County_ANSI = r.cd_fips_county
join ref_lookup_measure rlm
on rlm.cd_measure = r.cd_measure
order by
month, tx_state, tx_county")
return(results)
}
dat <- getLAUCN(states = "WA", counties = "*", measures = "03", startyear = "2000", endyear = "2013")
head(dat)
?bic.surv
?bic.survival
require(knitr)
#source("tex_convert.R")
require(Hmisc)
require(RColorBrewer)
require(Benchmarking)
require(sqldf)
require(BMA)
#require(Amelia) # generate multiple imputations
require(mitools) # for MIextract()
require(mix) # for mi.inference()
require(ggplot2)
#require(extrafont)
require(pocr)
require(gridExtra)
require(plyr)
require(xtable)
require(knitcitations)
biblio <- read.bib("qualpaper_v5.bib")
cols <- brewer.pal(8,"Set1")
opts_chunk$set(echo = FALSE
,warning=FALSE
,results='asis'
,message=FALSE
,error=FALSE
,dpi = 400)
#clear memory
#rm(list=ls(all=TRUE))
setwd("C:/Users/mienkoja/Dropbox/qualpaper")
#setwd("~/Dropbox/qualpaper/")
set.seed(123456)
#load("~/Dropbox/qualpaper/sech_out.RData")
load("C:/Users/mienkoja/Dropbox/repos/qualpaper/sech_out.RData")
# see page 215 from Bogetoft and Otto 2011
sfa_dat <- as.data.frame(na.omit(with(r_dat, cbind(id
,w_ta=w_ta*alpha
,t_ta=t_ta*alpha
,t_tvc
,x_c)
)
)
)
sfa_dat <- subset(sfa_dat, (is.infinite(sfa_dat$w_ta)==FALSE | is.infinite(sfa_dat$t_ta==FALSE)
)
)
sfa_dat <- subset(sfa_dat, !(sfa_dat$w_ta==0 & sfa_dat$t_ta==0))
x <- with(sfa_dat, cbind(w_ta, t_ta))
y1 <- matrix(sfa_dat$t_tvc)
y2 <- matrix(sfa_dat$x_c)
t_tvc_sfa <- sfa(log1p(x), log1p(y1))
summary(t_tvc_sfa)
#percentage of inefficiency variation to total variation
lambda <- lambda.sfa(t_tvc_sfa)
100*lambda^2/(1+lambda^2)
#variance for inefficiency
sigma2u.sfa(t_tvc_sfa)
#variance for random errors
sigma2v.sfa(t_tvc_sfa)
#residuals
e <- residuals(t_tvc_sfa)
#sigma 2
s2 <- sigma2.sfa(t_tvc_sfa)
mustar <- -e*lambda^2/(1+lambda^2)
sstar <- lambda/(1+lambda^2)*sqrt(s2)
tej <- exp(-mustar-sstar*(dnorm(mustar/sstar)/pnorm(mustar/sstar)))
tejt <- data.frame(id=sfa_dat$id, tejt=tej[1:1822])
#tejt <- data.frame(id=sfa_dat$id, tejt=tej[1:927])
r_dat <- sqldf("select
r.*
,tejt
from r_dat r
left join tejt tt
on r.id=tt.id")
#try benchmarking for x_c
x_c_sfa <- sfa(log1p(x), log1p(y2))
summary(x_c_sfa)
#percentage of inefficiency variation to total variation
lambda <- lambda.sfa(x_c_sfa)
100*lambda^2/(1+lambda^2)
#variance for inefficiency
sigma2u.sfa(x_c_sfa)
#variance for random errors
sigma2v.sfa(x_c_sfa)
#residuals
e <- residuals(x_c_sfa)
#sigma 2
s2 <- sigma2.sfa(x_c_sfa)
mustar <- -e*lambda^2/(1+lambda^2)
sstar <- lambda/(1+lambda^2)*sqrt(s2)
tej <- exp(-mustar-sstar*(dnorm(mustar/sstar)/pnorm(mustar/sstar)))
#tejx <- data.frame(id=sfa_dat$id, tejx=tej[1:927])
tejx <- data.frame(id=sfa_dat$id, tejx=tej[1:1822])
r_dat <- sqldf("select
r.*
,tejx
from r_dat r
left join tejx tx
on r.id=tx.id")
X <- with(sfa_dat, cbind(w_ta = w_ta/t_ta, t_tvc, x_c))
Y <- matrix(sfa_dat$t_ta, ncol=1)
dist <- sfa(log1p(X), -log(Y))
tedist <- te.sfa(dist)
sigma2u <- sigma2u.sfa(dist)
sigma2v <- sigma2v.sfa(dist)
te <- data.frame(id=sfa_dat$id, te=tedist)
# commented code is to possibly model te as a hyper parameter
#nsim <- 10
#nobs <- 1823
#te <- data.frame(id=rep(sfa_dat$id, nsim), te=rep(tedist, nsim))
#err <- rep(NA, nsim*nobs)
#err <- rnorm(1, mean=0, sd=sigma2v)-rnorm(1, mean=0, sd=sigma2u), nsim*nobs)
# for (i in 1:(nsim*nobs)){
#     err[i] <- rnorm(1, mean=0, sd=sigma2v)-rnorm(1, mean=0, sd=sigma2u)
# }
# te$err <- err
# te$te_sim <- te$te + te$err
#
# idx <- data.frame(j=rep(seq(1:nobs), nsim), i = rep(seq(1:nsim), nobs), tot = seq(1:(nsim*nobs)))
#
# m <- lm(tot~j+i, dat=idx)
r_dat <- sqldf("select
r.*
,te
from r_dat r
left join te te
on r.id=te.id")
r_dat_sub <- subset(r_dat
,r_dat$c_age > 18
,select = c(w_ta
,neg_count
,neg_count_alt1
,neg_count_alt2
,pos_count
,pos_count_alt1
,pos_count_alt2
,alpha
,te
,c_age
,w_ta
,cnt_ch
,m_white
,m_age
,m_mar
,m_college
,m_hi_frus
,c_health
,c_hi_health
,dev_cnc)
)
#transform some variables
r_dat_sub$te <- as.numeric(r_dat_sub$te)
r_dat_sub$log_w_ta <- log(r_dat_sub$w_ta)
r_dat_sub$log_m_age <- log(r_dat_sub$m_age)
r_dat_sub$log_c_age <- log(r_dat_sub$c_age)
#final calculation for p_all_neg
r_dat_sub$p_all_neg <- (r_dat_sub$neg_count/2)/((r_dat_sub$neg_count/2)+(r_dat_sub$pos_count/3))
r_dat_sub$p_all_neg_alt1 <- (r_dat_sub$neg_count/3)/((r_dat_sub$neg_count/3)+(r_dat_sub$pos_count/2))
r_dat_sub$p_all_neg_alt2 <- (r_dat_sub$neg_count/2)/((r_dat_sub$neg_count/2)+(r_dat_sub$pos_count/2))
r_dat_sub$p_all_neg_d <- ((r_dat_sub$neg_count/2)+(r_dat_sub$pos_count/3))
r_dat_sub$p_all_neg_alt1_d <- ((r_dat_sub$neg_count/3)+(r_dat_sub$pos_count/2))
r_dat_sub$p_all_neg_alt2_d <- ((r_dat_sub$neg_count/2)+(r_dat_sub$pos_count/2))
#calculate some interactions
r_dat_sub$alpha_by_log_w_ta <- r_dat_sub$alpha*log(r_dat_sub$w_ta)
r_dat_sub$te_by_log_w_ta <- r_dat_sub$te*log(r_dat_sub$w_ta)
r_dat_sub$alpha_by_te <- r_dat_sub$te*r_dat_sub$alpha
r_dat_sub$alpha_by_log_w_ta_by_te <- r_dat_sub$te*log(r_dat_sub$w_ta)*r_dat_sub$alpha
x_disp=subset(r_dat_sub, select=c(p_all_neg
,alpha
,te
,w_ta
,cnt_ch
,c_age
,m_white
,m_age
,m_mar
,m_college
,m_hi_frus
,c_hi_health
,dev_cnc))
dtf <- t(sapply(na.omit(x_disp), each(min, max, mean, median)))
colnames(dtf) <- c("Min","Max","Mean", "Median")
rownames(dtf) <- c("Probability of All Type II"
,"Altruism"
,"Efficiency^[Ibid.]"
,"Income"
,"Child Count"
,"Child Age (mos)", "White Mother", "Maternal Age", "Married Mother"
,"Maternal College", "Maternal Frustration", "Child Healthy", "Devolpmental Concerns")
#dtf <- xtable(dtf, label='tabdesc',caption='Descriptive Statistics', hline.after=c(0))
#align(dtf) <- "crrrr"
#print(dtf, sanitize.text.function = function(x){x}, scalebox=1, type='html')
#kable(round(dtf, 2))
# m1 <- glm(p_all_neg ~ alpha +
#                  log_w_ta
#                  ,family=quasibinomial
#                  ,weights=r_dat_sub$p_all_neg_d
#                  ,data=r_dat_sub)
# names(m1$coefficients) <- c("Intercept", "Altruism", "Income")
# #colnames(coefficients(summary(m1)))  <- c("$\\Beta{X}$", "SE", "$t$", "$p$")
# #m1_tbl <- xtable(m1,label='tabmod',caption='Model Results', hline.after=c(0))
# #print(m1_tbl, scalebox=1, type='html')
# kable(coefficients(summary(m1))[,1:3])
#
# m1_alt1 <- glm(p_all_neg_alt1 ~ alpha +
#                  log_w_ta
#                  ,family=quasibinomial
#                  ,weights=r_dat_sub$p_all_neg_alt1_d
#                  ,data=r_dat_sub)
# names(m1_alt1$coefficients) <- c("Intercept", "Altruism", "Income")
# #colnames(coefficients(summary(m1)))  <- c("$\\Beta{X}$", "SE", "$t$", "$p$")
# #m1_tbl <- xtable(m1,label='tabmod',caption='Model Results', hline.after=c(0))
# #print(m1_tbl, scalebox=1, type='html')
# kable(coefficients(summary(m1_alt1))[,1:3])
m1_alt2 <- glm(p_all_neg_alt2 ~ alpha +
log_w_ta
,family=quasibinomial
,weights=r_dat_sub$p_all_neg_alt2_d
,data=r_dat_sub)
names(m1_alt2$coefficients) <- c("Intercept", "Altruism", "Income")
#colnames(coefficients(summary(m1)))  <- c("$\\Beta{X}$", "SE", "$t$", "$p$")
#m1_tbl <- xtable(m1,label='tabmod',caption='Model Results', hline.after=c(0))
#print(m1_tbl, scalebox=1, type='html')
kable(coefficients(summary(m1_alt2))[,1:3])
# m1 <- glm(p_all_neg ~ alpha +
#                  log_w_ta
#                  ,family=quasibinomial
#                  ,weights=r_dat_sub$p_all_neg_d
#                  ,data=r_dat_sub)
# names(m1$coefficients) <- c("Intercept", "Altruism", "Income")
# #colnames(coefficients(summary(m1)))  <- c("$\\Beta{X}$", "SE", "$t$", "$p$")
# #m1_tbl <- xtable(m1,label='tabmod',caption='Model Results', hline.after=c(0))
# #print(m1_tbl, scalebox=1, type='html')
# kable(coefficients(summary(m1))[,1:3])
#
# m1_alt1 <- glm(p_all_neg_alt1 ~ alpha +
#                  log_w_ta
#                  ,family=quasibinomial
#                  ,weights=r_dat_sub$p_all_neg_alt1_d
#                  ,data=r_dat_sub)
# names(m1_alt1$coefficients) <- c("Intercept", "Altruism", "Income")
# #colnames(coefficients(summary(m1)))  <- c("$\\Beta{X}$", "SE", "$t$", "$p$")
# #m1_tbl <- xtable(m1,label='tabmod',caption='Model Results', hline.after=c(0))
# #print(m1_tbl, scalebox=1, type='html')
# kable(coefficients(summary(m1_alt1))[,1:3])
m1_alt2 <- glm(p_all_neg_alt2 ~ alpha +
log_w_ta
,family=quasibinomial
,weights=r_dat_sub$p_all_neg_alt2_d
,data=r_dat_sub)
names(m1_alt2$coefficients) <- c("Intercept", "Altruism", "Income")
#colnames(coefficients(summary(m1)))  <- c("$\\Beta{X}$", "SE", "$t$", "$p$")
#m1_tbl <- xtable(m1,label='tabmod',caption='Model Results', hline.after=c(0))
#print(m1_tbl, scalebox=1, type='html')
xtable(coefficients(summary(m1_alt2))[,1:3])
m1_tbl <- xtable(m1_alt2,label='tabmod',caption='Model Results', hline.after=c(0))
print(m1_tbl, scalebox=1)
library(QuantPsyc)
lm.beta(m1_alt2)
m1_alt2
lm.beta(m1_alt2)
lm.beta(m1_alt2)[1]
std_alpha <- lm.beta(m1_alt2)[1]
lm.beta(m1_alt2)[2]
std_altruism <- lm.beta(m1_alt2)[1]
std_income <- lm.beta(m1_alt2)[2]
std_income
sim_dat_w3
sim_dat <- subset(r_dat_sub, select=c(alpha, log_w_ta))
sim_dat_w1 <- with(sim_dat
,data.frame(alpha = mean(sim_dat$alpha, na.rm=TRUE)
,log_w_ta = rep(seq(from = 4, to = 13, length.out = 1000))
)
)
sim_dat_w2 <- cbind(sim_dat_w1
,predict(m1_alt2, type="response", newdata=sim_dat_w1, se = TRUE))
sim_dat_w3 <- within(sim_dat_w2, {
LL <- fit - (1.96 * se.fit)
UL <- fit + (1.96 * se.fit)
})
sim_dat_a1 <- with(sim_dat
,data.frame(log_w_ta = mean(sim_dat$log_w_ta, na.rm=TRUE)
,alpha = rep(seq(from = 0, to = 1, length.out = 1000))
)
)
sim_dat_a2 <- cbind(sim_dat_a1
,predict(m1_alt2, type="response", newdata=sim_dat_a1, se = TRUE))
sim_dat_a3 <- within(sim_dat_a2, {
LL <- fit - (1.96 * se.fit)
UL <- fit + (1.96 * se.fit)
})
breaks=c(.40, .45, .50, .55, .60)
#
#
#
w_p <- ggplot(sim_dat_w3, aes(x = log_w_ta, y = fit)) +
geom_ribbon(aes(ymin = LL, ymax = UL),alpha = 0.2, fill=cols[2]) +
geom_line(size = 1, colour=cols[2]) +
xlab("Log of Income") +
ylab("Probability of All Type II Discipline") +
scale_y_continuous(labels = percent
#,breaks=breaks
,limits=c(0, 1)) +
theme_bw() +
theme(text=element_text(size=15))
w_p
a_p <- ggplot(sim_dat_a3, aes(x = alpha, y = fit)) +
geom_ribbon(aes(ymin = LL, ymax = UL),alpha = 0.2, fill=cols[2]) +
geom_line(size = 1, colour=cols[2]) +
xlab("Altruism") +
ylab("Probability of All Type II Discipline") +
scale_y_continuous(labels = percent
#,breaks=breaks
,limits=c(0, 1)) +
theme_bw() +
theme(text=element_text(size=15))
a_p
cbind(sim_dat_a3, sim_dat_w3)
head(cbind(sim_dat_a3, sim_dat_w3))
head(rbind(sim_dat_a3, sim_dat_w3))
std_altruism <- round(lm.beta(m1_alt2)[1], 2)
std_income <- round(lm.beta(m1_alt2)[2], 2)
round(lm.beta(m1_alt2)[1], 2)