R/parameter_distr_plots.R
In davewutchiett/nets_demo_spatial_dynamics: Predicting Change and Modeling Relationships In NYC Local Retail and Population.
# coefficient estimates
library(plyr)
library(gplots)
coef_processes = function(set = 'pa'){
if (set == 'pa'){
boot.coefs = readRDS('code/pacoefs_080515.rds')
boot.models = readRDS('code/pamodels_080515.rds')
data.orig = read.csv('code/data_pa.csv')
data.orig = data.orig[,! colnames(data.orig) %in% "yc.pa.2"]
} else if (set == 'hf') {
boot.coefs = readRDS('code/hfcoefs_080515.rds')
boot.models = readRDS('code/hfmodels_080515.rds')
data.orig = read.csv('code/data_hf.csv')
data.orig = data.orig[,! colnames(data.orig) %in% "yc.hf.2"]
} else if (set == 'pop') {
boot.coefs = readRDS('code/popcoefs2015-08-13.rds')
boot.models = readRDS('code/popmodels2015-08-13.rds')
data.orig = read.csv('code/data_pop.csv')
data.orig = data.orig[,! colnames(data.orig) %in% "yc.pop.2"]
}
mod.params = nrow(boot.models[[1]])/2
pen.emp = ldply(boot.coefs, function(x) rbind(x$t[,1:mod.params]))
mle.emp = ldply(boot.coefs, function(x) rbind(x$t[,(mod.params+1):(2*mod.params)]))
ci.pen = lapply(pen.emp, function(x) paste(round(quantile(x,probs=c(.25,.975),na.rm=TRUE),4)))
ci.pen = ldply(ci.pen,rbind)
median.pen = lapply(pen.emp, function(x) round(median(x),4))
mean.pen = lapply(pen.emp, function(x) round(mean(x),4))
pos.dir = as.numeric(paste(ci.pen[,2])) >0 & as.numeric(paste(ci.pen[,3])) >0
neg.dir = as.numeric(paste(ci.pen[,2])) <0 & as.numeric(paste(ci.pen[,3])) <0
pen.sig = pos.dir | neg.dir
pen = cbind(median.pen,as.numeric(paste(ci.pen[,2])),as.numeric(paste(ci.pen[,3])),pen.sig)
ci.mle = lapply(mle.emp, function(x) paste(round(quantile(x,probs=c(.25,.975),na.rm=TRUE),4)))
ci.mle = ldply(ci.mle,rbind)
median.mle = lapply(mle.emp, function(x) round(median(x),4))
mean.mle = lapply(mle.emp, function(x) round(mean(x),4))
pos.dir = as.numeric(paste(ci.mle[,2])) >0 & as.numeric(paste(ci.mle[,3])) >0
neg.dir = as.numeric(paste(ci.mle[,2])) <0 & as.numeric(paste(ci.mle[,3])) <0
mle.sig = pos.dir | neg.dir
mle = cbind(median.mle, as.numeric(paste(ci.mle[,2])),as.numeric(paste(ci.mle[,3])), mle.sig)
names = c(colnames(data.orig), c('nd','rd','aic','rsq','cv_md','cv_ce'))
key_ests = cbind(names,pen,mle)
# pick out variables with non-null estimates
either.sig = key_ests[,'pen.sig'] == TRUE | key_ests[,'mle.sig'] == TRUE
kept = as.data.frame(key_ests[either.sig == TRUE,])
rows = nrow(kept)
plot.row = nrow(kept) - 6
boot.out = as.matrix(kept[c(2:(rows-6)),])
boot.out[,c(2:4,6:8)] = exp(as.numeric(boot.out[,c(2:4,6:8)]))
asn = function(x){as.numeric(x)}
par(oma=c(0,18,0,0))
x.u.lim = max(unlist(boot.out[,4])) + .1*max(unlist(boot.out[,4]))
x.l.lim = min(unlist(boot.out[,4])) - .1*max(unlist(boot.out[,4]))
barplot2((asn(boot.out[,2])), plot.ci= TRUE, ci.l= (asn(boot.out[,3])),
ci.u= (asn(boot.out[,4])),horiz= TRUE,names.arg=boot.out[,1],xlim=c(x.l.lim,x.u.lim) , xlab = "Odds Ratio",
col="white",border="white",las=1,space = 0) #, main="Bootstrap 95% Percentile CI for Model Coefficients (2000 samples)")
abline(v=1)
abline(h=1:plot.row-1,col="light gray")
#abline(h=(2:plot.row-1.5), col="light gray",lty=2)
points((asn(boot.out[,2])),(2:(rows-6))-1.4,pch="|",cex=1)
x.u.lim = max(unlist(boot.out[,8])) + .1*max(unlist(boot.out[,8]))
x.l.lim = min(unlist(boot.out[,8])) - .1*max(unlist(boot.out[,8]))
par(oma=c(0,18,0,0))
barplot2((asn(boot.out[,6])), plot.ci= TRUE, ci.l= (asn(boot.out[,7])),
ci.u= (asn(boot.out[,8])),horiz= TRUE,names.arg=boot.out[,1],xlim=c(x.l.lim,x.u.lim) , xlab = "Odds Ratio",
col="white",border="white",las=1,space=0) #, main="Bootstrap 95% Percentile CI for Model Coefficients (2000 samples)")
abline(v=1)
abline(h=1:plot.row-1,col="light gray")
#abline(h=(2:plot.row-1.5), col="light gray",lty=2)
points((asn(boot.out[,6])),(2:(rows-6))-1.4,pch="|",cex=1)
#pre = c('_','_x',"tX","\\_adj","\\_c",
# ' ',' by',"time by"," Adjacent"," Change",
pre = c("hinc","nFlagged","\\_pa","pasian","phisp","phys_act_b0","pnhblk",
"pop", "ppov", "tr_size", "adj_tracts", "pfb", "phys_act", "healthy_foods", "nBusinesses")
post= c("Median Income","Other Businesses","","% Asian","% Hispanic",
"Previous Physical Activity Facilities",
"% Non-Hispanic Black","Population","% Poverty","Tract Area","Total Adjacent Tracts",
"% Foreign Born", "Physical Activity Facilities", "Healthy Food Outlets", "Businesses (alternate Types)")
boot.out[,1] = sub('tX', 'Time by ', boot.out[,1])
boot.out[,1] = sub('pop_adj', 'adjacent, weighted by population', boot.out[,1])
boot.out[,1] = sub('pop_c_adj', 'adjacent, weighted by population change', boot.out[,1])
boot.out[,1] = sub('_x_', ' by ', boot.out[,1])
boot.out[,1] = sub('_adj', ' adjacent', boot.out[,1])
boot.out[,1] = sub('_c_', ' change ', boot.out[,1])
boot.out[,1] = sub('_c', ' change ', boot.out[,1])
for (i in 1:length(pre)){
boot.out[,1] = sub(pre[i], post[i], boot.out[,1])
}
for (i in 1:length(pre)){
dem.out[,1] = sub(pre[i],post[i],dem.out[,1])
env.out[,1] = sub(pre[i],post[i],env.out[,1])
both.out[,1] = sub(pre[i],post[i],both.out[,1])
}
dem.out = dem.out[nrow(dem.out):1,]
env.out = env.out[nrow(env.out):1,]
both.out = both.out[nrow(both.out):1,]
### plot coefficients
ci_indices = seq(1:rows)
#boot.ci(nets.boot,index=1, type=c("perc"))
#ci_vals_l = unlist(lapply(ci_indices, function(x) boot.ci(nets.boot,index=x, type=c("perc"))$percent[c(4)]))
#ci_vals_u = unlist(lapply(ci_indices, function(x) boot.ci(nets.boot,index=x, type=c("perc"))$percent[c(5)]))
#coef_lab = labels(me_mod_coef)[[1]]
#labs = c(coef_lab) #,"r2" ,"lambda"
full_var_names = c("2000-2010","Adjacent Tract Size (sq. miles)","Tract Size (sq. miles)","% Poverty Change",
"Adjacent % Poverty Change", "% Poverty Adjacent", "% Poverty", "Population Change", "Population Change Adjacent",
"Population Adjacent", "Population","% Non-Hispanic Black Change","% Non-Hispanic Black Adjacent Change",
"% Non-Hispanic Black Adjacent","% Non-Hispanic Black", "% Hispanic Change","% Hispanic Adjacent Change",
"% Hispanic Adjacent", "% Hispanic","% Foreign-born Change", "% Foreign-born Change Adjacent", "% Foreign-born Adjacent","% Foreign-born",
"% Asian Change", "% Asian Adjacent Change", "% Asian Adjacent",
"% Asian", "No Adjacent Tracts", "Other Businesses Change","Other Businesses", "Other Businesses Adjacent Change",
"Other Businesses Adjacent","Median Income Change","Median Income Change Adjacent", "Median Income Adjacent",
"Median Income","Previous Healthy Food Outlets","Total Adjacent Tracts")
full_var_names = rev(full_var_names)
boot.out = cbind(labs,nets.boot$t0,ci_vals_l,ci_vals_u)
boot.out = boot.out[c(2:39),]
boot.out.plot = boot.out[order(boot.out[,1]),]
boot.out.plot[,1] = full_var_names
boot.out.plot[,2:4] = apply(boot.out.plot[,2:4],2, function(x) exp(as.numeric(x)))
asn = function(x){as.numeric(x)}
library(gplots)
par(oma=c(0,18,0,0))
barplot2(asn(boot.out.plot[,2]), plot.ci= TRUE, ci.l= asn(boot.out.plot[,3]),
ci.u= asn(boot.out.plot[,4]),horiz= TRUE,names.arg=boot.out.plot[,1],xlim=c(.4,2.75) , xlab = "Odds Ratio",
col="white",border="white",las=1,space=0) #, main="Bootstrap 95% Percentile CI for Model Coefficients (2000 samples)")
abline(v=1)
abline(h=0:38,col="light gray")
abline(h=(1:38-.5),col="light gray",lty=2)
points(asn(boot.out.plot[,2]),(1:38-.5),pch="|",cex=1.2)
pdf("Y2_coef_plot.pdf", width=8.07/2.54, height=11.07/2.54)
###
full_var_names = c("2000-2010","Tract Area Adjacent","Tract Area","% Poverty Change",
"% Poverty Adjacent Change", "% Poverty Adjacent", "% Poverty", "Population Change", "Population Adjacent Change",
"Population Adjacent", "Population","% Non-Hispanic Black Change","% Non-Hispanic Black Adjacent Change",
"% Non-Hispanic Black Adjacent","% Non-Hispanic Black", "Previous Physical Activity Facilities", "% Hispanic Change","% Hispanic Adjacent Change",
"% Hispanic Adjacent", "% Hispanic", "% Asian Change", "% Asian Adjacent Change", "% Asian Adjacent",
"% Asian", "No Adjacent Tracts", "Other Businesses Change","Other Businesses", "Other Businesses Adjacent Change",
"Other Businesses Adjacent","Median Income Change","Median Income Adjacent Change", "Median Income Adjacent",
"Median Income","Total Adjacent Tracts")
full_var_names = rev(full_var_names)
boot.out = cbind(labs,nets.boot$t0,ci_vals_l,ci_vals_u)
boot.out = boot.out[c(3:36),]
boot.out.plot = boot.out[order(boot.out[,1]),]
boot.out.plot[,1] = full_var_names
boot.out.plot[,2:4] = apply(boot.out.plot[,2:4],2, function(x) exp(as.numeric(x)))
#colnames(boot.out.plot) = c('labs','or','ci_vals_l','ci_vals_u')
#this had removed two outliers, but perhaps split to
#boot.out.plot = boot.out.plot[boot.out.plot[,2]!=1&boot.out.plot[,2]<1.06,]
env_chars = c("Previous Physical Activity Facilities", "Other Businesses", "Other Businesses Change", "Other Businesses Adjacent",
"Other Businesses Adjacent Change", "Tract Area", "Tract Area Adjacent",
"Total Adjacent Tracts","No Adjacent Tracts")
env_output = boot.out.plot[boot.out.plot[,1]%in%env_chars,]
dem_output = boot.out.plot[!boot.out.plot[,1]%in%env_chars,]
boot.out.plot = env_output
asn = function(x){as.numeric(x)}
par(oma=c(0,12,0,0))
barplot2(asn(boot.out.plot[,2]), plot.ci= TRUE, ci.l= asn(boot.out.plot[,3]),
ci.u= asn(boot.out.plot[,4]),horiz= TRUE,names.arg=boot.out.plot[,1],xlim=c(.95,1.135),
col="white",border="white",las=1,space=0) #, main="Bootstrap 95% Percentile CI for Model Coefficients (2000 samples)")
abline(v=1)
abline(h=0:nrow(boot.out.plot),col="light gray")
abline(h=(1:nrow(boot.out.plot)-.5),col="light gray",lty=2)
points(asn(boot.out.plot[,2]),(1:nrow(boot.out.plot)-.5),pch="|",cex=2)
boot.out.plot = dem_output
asn = function(x){as.numeric(x)}
par(oma=c(0,12,0,0))
barplot2(asn(boot.out.plot[,2]), plot.ci= TRUE, ci.l= asn(boot.out.plot[,3]),
ci.u= asn(boot.out.plot[,4]),horiz= TRUE,names.arg=boot.out.plot[,1],xlim=c(.98,1.011),
col="white",border="white",las=1,space=0) #, main="Bootstrap 95% Percentile CI for Model Coefficients (2000 samples)")
abline(v=1)
abline(h=0:nrow(boot.out.plot),col="light gray")
abline(h=(1:nrow(boot.out.plot)-.5),col="light gray",lty=2)
points(asn(boot.out.plot[,2]),(1:nrow(boot.out.plot)-.5),pch="|",cex=2)
####################################################
#names(X) = c("Intercept","Decade Start","Population","Median Household Income","Physical Activity Facilities","Change in Population",
# "Change in Median Household Income","Tract Size")
aa = coef(fit_lambda,lambdaIndex = lamda_indice)[1]$mainEffects$cat
aa.vals = colnames(X)[aa+1]
bb = coef(fit_lambda,lambdaIndex = lamda_indice)[2]$mainEffectsCoef$cat
bb.list = unlist(lapply(bb,function(x) round(exp(x),3)))
cat.main.table = cbind(c("1990-2000","2000-2010"),c(bb.list))
a = coef(fit_lambda,lambdaIndex = lamda_indice)[1]$mainEffects$cont
a.vals = colnames(X)[a+length(aa+1)]
b = coef(fit_lambda,lambdaIndex = lamda_indice)[2]$mainEffectsCoef$cont
b.list = lapply(b,function(x) round(exp(x),3))
#main effects
main.table = as.data.frame(cbind(a.vals,b.list))
colnames(main.table) = c('var','ests')
main.table$var= as.character(main.table$var)
main.table= main.table[order(main.table$var),]
main.table = as.matrix(main.table)
outcome_fill = mat.or.vec((length(coef(fit_lambda,lambdaIndex = lamda_indice)$interactions$contcont[,1])),1)
for (i in 1:length(outcome_fill)){
c = coef(fit_lambda,lambdaIndex = lamda_indice)$interactions$contcont[i,1]
d = coef(fit_lambda,lambdaIndex = lamda_indice)$interactions$contcont[i,2]
outcome_fill[i] = paste(colnames(X)[c+2],colnames(X)[d+2],sep=" X ")
}
inter_values = unlist(coef(fit_lambda,lambdaIndex = lamda_indice)$interactionsCoef$contcont)
inter_values = lapply(inter_values,function(x) round(exp(x),3))
inter.table = cbind(outcome_fill,inter_values)
outcome_fill2a = mat.or.vec((length(coef(fit_lambda,lambdaIndex = lamda_indice)$interactions$catcont[,1])),1)
for (i in 1:length(outcome_fill2a)){
c = coef(fit_lambda,lambdaIndex = lamda_indice)$interactions$catcont[i,1]
d = coef(fit_lambda,lambdaIndex = lamda_indice)$interactions$catcont[i,2]+length(aa)+1
outcome_fill2a[i] = paste('1990-2000','x',colnames(X)[d],sep=" ") #names(X)[c],'
}
outcome_fill2b = mat.or.vec((length(coef(fit_lambda,lambdaIndex = lamda_indice)$interactions$catcont[,1])),1)
for (i in 1:length(outcome_fill2b)){
c = coef(fit_lambda,lambdaIndex = lamda_indice)$interactions$catcont[i,1]
d = coef(fit_lambda,lambdaIndex = lamda_indice)$interactions$catcont[i,2]+length(aa)+1
outcome_fill2b[i] = paste('2000-2010','x',colnames(X)[d],sep=" ")#names(X)[c],
}
outcome_fill2 = rbind(t(t(outcome_fill2a)),t(t(outcome_fill2b)))
cat.inter_values = unlist(coef(fit_lambda,lambdaIndex = lamda_indice)$interactionsCoef$catcont)
cat.inter_values = lapply(cat.inter_values,function(x) round(exp(x),3))
cat.inter_values.transfer = c(cat.inter_values[1],cat.inter_values[2])
cat.inter.table = cbind(outcome_fill2,cat.inter_values.transfer)
cat.inter.table = rbind(cat.inter.table[1:5,],c('',''),cat.inter.table[6:10,])
rbind(cat.main.table,main.table,inter.table,cat.inter.table)
output_table = rbind(cat.main.table,main.table,c('',''),c('Interaction Effects',''),inter.table,c('',''),cat.inter.table)
colnames(output_table) = c("Parameter", "OR")
output_table = rbind(output_table,cbind(c('','Pseudo R^2'),c('',pR2)))
output_table = rbind(output_table,cbind(c('','R^2'),c('',r2)))
output_table[,2] = log(as.numeric(output_table[,2]))
output_table
#write.csv(output_table,"physact_inter_020515.csv") #physact_inter_020515.csv
inter_plot = output_table[c(1:18,21:36,38,39),c(1,2)]
inter_plot[,2] = as.numeric(inter_plot[,2])
par(oma=c(0,12,0,0))
barplot2(asn(inter_plot[,2]), plot.ci= FALSE, horiz= TRUE,names.arg=inter_plot[,1],xlim=c(-.04,.67) ,
col="white",border="white",las=1,space=0, main="Coefficient Estimates for Total Commercial Physical Activity Facilities through Interactions")
abline(v=0)
abline(h=0:36,col="light gray")
abline(h=(1:36-.5),col="light gray",lty=2)
points((inter_plot[,2]),(1:36-.5),pch="|")
##### adjust for dichotmous
aa = coef(fit_lambda2,lambdaIndex = lamda_indice2)[1]$mainEffects$cat
aa.vals = colnames(X)[aa+1]
bb = coef(fit_lambda2,lambdaIndex = lamda_indice2)[2]$mainEffectsCoef$cat
bb.list = unlist(lapply(bb,function(x) round(exp(x),3)))
cat.main.table = cbind(c("1990-2000","2000-2010"),c(bb.list))
a = coef(fit_lambda2,lambdaIndex = lamda_indice2)[1]$mainEffects$cont
a.vals = colnames(X)[a+length(aa+1)]
b = coef(fit_lambda2,lambdaIndex = lamda_indice2)[2]$mainEffectsCoef$cont
b.list = lapply(b,function(x) round(exp(x),3))
#main effects
main.table = as.data.frame(cbind(a.vals,b.list))
colnames(main.table) = c('var','ests')
main.table$var= as.character(main.table$var)
main.table= main.table[order(main.table$var),]
main.table = as.matrix(main.table)
outcome_fill = mat.or.vec((length(coef(fit_lambda2,lambdaIndex = lamda_indice2)$interactions$contcont[,1])),1)
for (i in 1:length(outcome_fill)){
c = coef(fit_lambda2,lambdaIndex = lamda_indice2)$interactions$contcont[i,1]
d = coef(fit_lambda2,lambdaIndex = lamda_indice2)$interactions$contcont[i,2]
outcome_fill[i] = paste(colnames(X)[c+2],colnames(X)[d+2],sep=" X ")
}
inter_values = unlist(coef(fit_lambda2,lambdaIndex = lamda_indice2)$interactionsCoef$contcont)
inter_values = lapply(inter_values,function(x) round(exp(x),3))
inter.table = cbind(outcome_fill,inter_values)
outcome_fill2a = mat.or.vec((length(coef(fit_lambda2,lambdaIndex = lamda_indice2)$interactions$catcont[,1])),1)
for (i in 1:length(outcome_fill2a)){
c = coef(fit_lambda2,lambdaIndex = lamda_indice2)$interactions$catcont[i,1]
d = coef(fit_lambda2,lambdaIndex = lamda_indice2)$interactions$catcont[i,2]+length(aa)+1
outcome_fill2a[i] = paste('1990-2000','x',colnames(X)[d],sep=" ") #names(X)[c],'
}
outcome_fill2b = mat.or.vec((length(coef(fit_lambda2,lambdaIndex = lamda_indice2)$interactions$catcont[,1])),1)
for (i in 1:length(outcome_fill2b)){
c = coef(fit_lambda2,lambdaIndex = lamda_indice2)$interactions$catcont[i,1]
d = coef(fit_lambda2,lambdaIndex = lamda_indice2)$interactions$catcont[i,2]+length(aa)+1
outcome_fill2b[i] = paste('2000-2010','x',colnames(X)[d],sep=" ")#names(X)[c],
}
outcome_fill2 = rbind(t(t(outcome_fill2a)),t(t(outcome_fill2b)))
cat.inter_values = unlist(coef(fit_lambda2,lambdaIndex = lamda_indice2)$interactionsCoef$catcont)
cat.inter_values = lapply(cat.inter_values,function(x) round(exp(x),3))
cat.inter_values.transfer = c(cat.inter_values[1],cat.inter_values[2])
cat.inter.table = cbind(outcome_fill2,cat.inter_values.transfer)
cat.inter.table = rbind(cat.inter.table[1:5,],c('',''),cat.inter.table[6:10,])
rbind(cat.main.table,main.table,inter.table,cat.inter.table)
output_table = rbind(cat.main.table,main.table,c('',''),c('Interaction Effects',''),inter.table,c('',''),cat.inter.table)
colnames(output_table) = c("Parameter", "OR")
output_table = rbind(output_table,cbind(c('','Pseudo R^2'),c('',pR2)))
output_table = rbind(output_table,cbind(c('','R^2'),c('',r2)))
#write.csv(output_table,"physact_inter_021015_binom.csv") #physact_inter_020515.csv
output_table[,2] = round(as.numeric(output_table[,2]),3)
########## trying to quicken/automate subsetting of names & readable changes
output_table = output_table[!output_table[,1]%in%c("","Interaction Effects","R^2"),]
non.int = grep("X|x", output_table[,1],invert=TRUE,value=TRUE)
yes.int = grep("X|x", output_table[,1],value=TRUE)
demv = c("1990-2000|2000-2010|hinc|pasian|phisp|pnhblk|pop|ppov")
envv = c("nFlagged|phys_act_b0|tr_size|adj_tracts|no_adj")
main.dem = grep(demv,non.int,value=TRUE)
main.env =grep(envv,non.int,value=TRUE)
int.env.only = grep(demv,yes.int,invert=TRUE,value=TRUE)
int.dem.only =grep(envv,yes.int,invert=TRUE,value=TRUE)
int.only.either = c(int.env.only,int.dem.only)
int.both.env.dem = grep(paste(int.only.either,collapse="|"),yes.int,invert=TRUE,value=TRUE)
dem.out = output_table[output_table[,1]%in%c(int.dem.only,main.dem),]
env.out = output_table[output_table[,1]%in%c(int.env.only,main.env),]
both.out = output_table[output_table[,1]%in%c(int.both.env.dem),]
pre = c("X","\\.adj","\\.c","hinc","\\.","nFlagged","\\_pa","pasian","phisp","phys_act_b0","pnhblk",
"pop","ppov","tr_size","adj_tracts","no_adj"," adj"," c", "pnhblk")
post= c("x"," Adjacent"," Change","Median Income"," ","Other Businesses","","% Asian","% Hispanic","Previous Physical Activity Facilities",
"% Non-Hispanic Black","Population","% Poverty","Tract Area","Total Adjacent Tracts","No Adjacent Tracts"," Adjacent"," Change","% Non-Hispanic Black")
for (i in 1:length(pre)){
dem.out[,1] = sub(pre[i],post[i],dem.out[,1])
env.out[,1] = sub(pre[i],post[i],env.out[,1])
both.out[,1] = sub(pre[i],post[i],both.out[,1])
}
dem.out = dem.out[nrow(dem.out):1,]
env.out = env.out[nrow(env.out):1,]
both.out = both.out[nrow(both.out):1,]
#inter_plot = output_table[output_table[,2]!=1,]
#inter_plot[,2] = as.numeric(inter_plot[,2])
tables = list(dem.out,env.out,both.out)
i=1
i=2
i=3
#for (i in 1:length(tables)){
inter_plot = as.data.frame(tables[i])
par(oma=c(0,25,0,0))
xmin = min(as.numeric(inter_plot[,2]))-.005
xmax = max(as.numeric(inter_plot[,2]))+.005
barplot2(asn(inter_plot[,2]), plot.ci= FALSE, horiz= TRUE,names.arg=inter_plot[,1],xlim=c(xmin,xmax) ,
col="white",border="white",las=1,space=0)#, main="OR Estimates for Increase in Commercial Physical Activity Facilities through Interactions")
abline(v=1)
abline(h=0:nrow(inter_plot),col="light gray")
abline(h=(1:nrow(inter_plot)-.5),col="light gray",lty=2)
points((inter_plot[,2]),(1:nrow(inter_plot)-.5),pch="|",cex=2)
#}
davewutchiett/nets_demo_spatial_dynamics documentation built on May 25, 2019, 4:22 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
# coefficient estimates
library(plyr)
library(gplots)
coef_processes = function(set = 'pa'){
if (set == 'pa'){
boot.coefs = readRDS('code/pacoefs_080515.rds')
boot.models = readRDS('code/pamodels_080515.rds')
data.orig = read.csv('code/data_pa.csv')
data.orig = data.orig[,! colnames(data.orig) %in% "yc.pa.2"]
} else if (set == 'hf') {
boot.coefs = readRDS('code/hfcoefs_080515.rds')
boot.models = readRDS('code/hfmodels_080515.rds')
data.orig = read.csv('code/data_hf.csv')
data.orig = data.orig[,! colnames(data.orig) %in% "yc.hf.2"]
} else if (set == 'pop') {
boot.coefs = readRDS('code/popcoefs2015-08-13.rds')
boot.models = readRDS('code/popmodels2015-08-13.rds')
data.orig = read.csv('code/data_pop.csv')
data.orig = data.orig[,! colnames(data.orig) %in% "yc.pop.2"]
}
mod.params = nrow(boot.models[[1]])/2
pen.emp = ldply(boot.coefs, function(x) rbind(x$t[,1:mod.params]))
mle.emp = ldply(boot.coefs, function(x) rbind(x$t[,(mod.params+1):(2*mod.params)]))
ci.pen = lapply(pen.emp, function(x) paste(round(quantile(x,probs=c(.25,.975),na.rm=TRUE),4)))
ci.pen = ldply(ci.pen,rbind)
median.pen = lapply(pen.emp, function(x) round(median(x),4))
mean.pen = lapply(pen.emp, function(x) round(mean(x),4))
pos.dir = as.numeric(paste(ci.pen[,2])) >0 & as.numeric(paste(ci.pen[,3])) >0
neg.dir = as.numeric(paste(ci.pen[,2])) <0 & as.numeric(paste(ci.pen[,3])) <0
pen.sig = pos.dir | neg.dir
pen = cbind(median.pen,as.numeric(paste(ci.pen[,2])),as.numeric(paste(ci.pen[,3])),pen.sig)
ci.mle = lapply(mle.emp, function(x) paste(round(quantile(x,probs=c(.25,.975),na.rm=TRUE),4)))
ci.mle = ldply(ci.mle,rbind)
median.mle = lapply(mle.emp, function(x) round(median(x),4))
mean.mle = lapply(mle.emp, function(x) round(mean(x),4))
pos.dir = as.numeric(paste(ci.mle[,2])) >0 & as.numeric(paste(ci.mle[,3])) >0
neg.dir = as.numeric(paste(ci.mle[,2])) <0 & as.numeric(paste(ci.mle[,3])) <0
mle.sig = pos.dir | neg.dir
mle = cbind(median.mle, as.numeric(paste(ci.mle[,2])),as.numeric(paste(ci.mle[,3])), mle.sig)
names = c(colnames(data.orig), c('nd','rd','aic','rsq','cv_md','cv_ce'))
key_ests = cbind(names,pen,mle)
# pick out variables with non-null estimates
either.sig = key_ests[,'pen.sig'] == TRUE | key_ests[,'mle.sig'] == TRUE
kept = as.data.frame(key_ests[either.sig == TRUE,])
rows = nrow(kept)
plot.row = nrow(kept) - 6
boot.out = as.matrix(kept[c(2:(rows-6)),])
boot.out[,c(2:4,6:8)] = exp(as.numeric(boot.out[,c(2:4,6:8)]))
asn = function(x){as.numeric(x)}
par(oma=c(0,18,0,0))
x.u.lim = max(unlist(boot.out[,4])) + .1*max(unlist(boot.out[,4]))
x.l.lim = min(unlist(boot.out[,4])) - .1*max(unlist(boot.out[,4]))
barplot2((asn(boot.out[,2])), plot.ci= TRUE, ci.l= (asn(boot.out[,3])),
ci.u= (asn(boot.out[,4])),horiz= TRUE,names.arg=boot.out[,1],xlim=c(x.l.lim,x.u.lim) , xlab = "Odds Ratio",
col="white",border="white",las=1,space = 0) #, main="Bootstrap 95% Percentile CI for Model Coefficients (2000 samples)")
abline(v=1)
abline(h=1:plot.row-1,col="light gray")
#abline(h=(2:plot.row-1.5), col="light gray",lty=2)
points((asn(boot.out[,2])),(2:(rows-6))-1.4,pch="|",cex=1)
x.u.lim = max(unlist(boot.out[,8])) + .1*max(unlist(boot.out[,8]))
x.l.lim = min(unlist(boot.out[,8])) - .1*max(unlist(boot.out[,8]))
par(oma=c(0,18,0,0))
barplot2((asn(boot.out[,6])), plot.ci= TRUE, ci.l= (asn(boot.out[,7])),
ci.u= (asn(boot.out[,8])),horiz= TRUE,names.arg=boot.out[,1],xlim=c(x.l.lim,x.u.lim) , xlab = "Odds Ratio",
col="white",border="white",las=1,space=0) #, main="Bootstrap 95% Percentile CI for Model Coefficients (2000 samples)")
abline(v=1)
abline(h=1:plot.row-1,col="light gray")
#abline(h=(2:plot.row-1.5), col="light gray",lty=2)
points((asn(boot.out[,6])),(2:(rows-6))-1.4,pch="|",cex=1)
#pre = c('_','_x',"tX","\\_adj","\\_c",
# ' ',' by',"time by"," Adjacent"," Change",
pre = c("hinc","nFlagged","\\_pa","pasian","phisp","phys_act_b0","pnhblk",
"pop", "ppov", "tr_size", "adj_tracts", "pfb", "phys_act", "healthy_foods", "nBusinesses")
post= c("Median Income","Other Businesses","","% Asian","% Hispanic",
"Previous Physical Activity Facilities",
"% Non-Hispanic Black","Population","% Poverty","Tract Area","Total Adjacent Tracts",
"% Foreign Born", "Physical Activity Facilities", "Healthy Food Outlets", "Businesses (alternate Types)")
boot.out[,1] = sub('tX', 'Time by ', boot.out[,1])
boot.out[,1] = sub('pop_adj', 'adjacent, weighted by population', boot.out[,1])
boot.out[,1] = sub('pop_c_adj', 'adjacent, weighted by population change', boot.out[,1])
boot.out[,1] = sub('_x_', ' by ', boot.out[,1])
boot.out[,1] = sub('_adj', ' adjacent', boot.out[,1])
boot.out[,1] = sub('_c_', ' change ', boot.out[,1])
boot.out[,1] = sub('_c', ' change ', boot.out[,1])
for (i in 1:length(pre)){
boot.out[,1] = sub(pre[i], post[i], boot.out[,1])
}
for (i in 1:length(pre)){
dem.out[,1] = sub(pre[i],post[i],dem.out[,1])
env.out[,1] = sub(pre[i],post[i],env.out[,1])
both.out[,1] = sub(pre[i],post[i],both.out[,1])
}
dem.out = dem.out[nrow(dem.out):1,]
env.out = env.out[nrow(env.out):1,]
both.out = both.out[nrow(both.out):1,]
### plot coefficients
ci_indices = seq(1:rows)
#boot.ci(nets.boot,index=1, type=c("perc"))
#ci_vals_l = unlist(lapply(ci_indices, function(x) boot.ci(nets.boot,index=x, type=c("perc"))$percent[c(4)]))
#ci_vals_u = unlist(lapply(ci_indices, function(x) boot.ci(nets.boot,index=x, type=c("perc"))$percent[c(5)]))
#coef_lab = labels(me_mod_coef)[[1]]
#labs = c(coef_lab) #,"r2" ,"lambda"
full_var_names = c("2000-2010","Adjacent Tract Size (sq. miles)","Tract Size (sq. miles)","% Poverty Change",
"Adjacent % Poverty Change", "% Poverty Adjacent", "% Poverty", "Population Change", "Population Change Adjacent",
"Population Adjacent", "Population","% Non-Hispanic Black Change","% Non-Hispanic Black Adjacent Change",
"% Non-Hispanic Black Adjacent","% Non-Hispanic Black", "% Hispanic Change","% Hispanic Adjacent Change",
"% Hispanic Adjacent", "% Hispanic","% Foreign-born Change", "% Foreign-born Change Adjacent", "% Foreign-born Adjacent","% Foreign-born",
"% Asian Change", "% Asian Adjacent Change", "% Asian Adjacent",
"% Asian", "No Adjacent Tracts", "Other Businesses Change","Other Businesses", "Other Businesses Adjacent Change",
"Other Businesses Adjacent","Median Income Change","Median Income Change Adjacent", "Median Income Adjacent",
"Median Income","Previous Healthy Food Outlets","Total Adjacent Tracts")
full_var_names = rev(full_var_names)
boot.out = cbind(labs,nets.boot$t0,ci_vals_l,ci_vals_u)
boot.out = boot.out[c(2:39),]
boot.out.plot = boot.out[order(boot.out[,1]),]
boot.out.plot[,1] = full_var_names
boot.out.plot[,2:4] = apply(boot.out.plot[,2:4],2, function(x) exp(as.numeric(x)))
asn = function(x){as.numeric(x)}
library(gplots)
par(oma=c(0,18,0,0))
barplot2(asn(boot.out.plot[,2]), plot.ci= TRUE, ci.l= asn(boot.out.plot[,3]),
ci.u= asn(boot.out.plot[,4]),horiz= TRUE,names.arg=boot.out.plot[,1],xlim=c(.4,2.75) , xlab = "Odds Ratio",
col="white",border="white",las=1,space=0) #, main="Bootstrap 95% Percentile CI for Model Coefficients (2000 samples)")
abline(v=1)
abline(h=0:38,col="light gray")
abline(h=(1:38-.5),col="light gray",lty=2)
points(asn(boot.out.plot[,2]),(1:38-.5),pch="|",cex=1.2)
pdf("Y2_coef_plot.pdf", width=8.07/2.54, height=11.07/2.54)
###
full_var_names = c("2000-2010","Tract Area Adjacent","Tract Area","% Poverty Change",
"% Poverty Adjacent Change", "% Poverty Adjacent", "% Poverty", "Population Change", "Population Adjacent Change",
"Population Adjacent", "Population","% Non-Hispanic Black Change","% Non-Hispanic Black Adjacent Change",
"% Non-Hispanic Black Adjacent","% Non-Hispanic Black", "Previous Physical Activity Facilities", "% Hispanic Change","% Hispanic Adjacent Change",
"% Hispanic Adjacent", "% Hispanic", "% Asian Change", "% Asian Adjacent Change", "% Asian Adjacent",
"% Asian", "No Adjacent Tracts", "Other Businesses Change","Other Businesses", "Other Businesses Adjacent Change",
"Other Businesses Adjacent","Median Income Change","Median Income Adjacent Change", "Median Income Adjacent",
"Median Income","Total Adjacent Tracts")
full_var_names = rev(full_var_names)
boot.out = cbind(labs,nets.boot$t0,ci_vals_l,ci_vals_u)
boot.out = boot.out[c(3:36),]
boot.out.plot = boot.out[order(boot.out[,1]),]
boot.out.plot[,1] = full_var_names
boot.out.plot[,2:4] = apply(boot.out.plot[,2:4],2, function(x) exp(as.numeric(x)))
#colnames(boot.out.plot) = c('labs','or','ci_vals_l','ci_vals_u')
#this had removed two outliers, but perhaps split to
#boot.out.plot = boot.out.plot[boot.out.plot[,2]!=1&boot.out.plot[,2]<1.06,]
env_chars = c("Previous Physical Activity Facilities", "Other Businesses", "Other Businesses Change", "Other Businesses Adjacent",
"Other Businesses Adjacent Change", "Tract Area", "Tract Area Adjacent",
"Total Adjacent Tracts","No Adjacent Tracts")
env_output = boot.out.plot[boot.out.plot[,1]%in%env_chars,]
dem_output = boot.out.plot[!boot.out.plot[,1]%in%env_chars,]
boot.out.plot = env_output
asn = function(x){as.numeric(x)}
par(oma=c(0,12,0,0))
barplot2(asn(boot.out.plot[,2]), plot.ci= TRUE, ci.l= asn(boot.out.plot[,3]),
ci.u= asn(boot.out.plot[,4]),horiz= TRUE,names.arg=boot.out.plot[,1],xlim=c(.95,1.135),
col="white",border="white",las=1,space=0) #, main="Bootstrap 95% Percentile CI for Model Coefficients (2000 samples)")
abline(v=1)
abline(h=0:nrow(boot.out.plot),col="light gray")
abline(h=(1:nrow(boot.out.plot)-.5),col="light gray",lty=2)
points(asn(boot.out.plot[,2]),(1:nrow(boot.out.plot)-.5),pch="|",cex=2)
boot.out.plot = dem_output
asn = function(x){as.numeric(x)}
par(oma=c(0,12,0,0))
barplot2(asn(boot.out.plot[,2]), plot.ci= TRUE, ci.l= asn(boot.out.plot[,3]),
ci.u= asn(boot.out.plot[,4]),horiz= TRUE,names.arg=boot.out.plot[,1],xlim=c(.98,1.011),
col="white",border="white",las=1,space=0) #, main="Bootstrap 95% Percentile CI for Model Coefficients (2000 samples)")
abline(v=1)
abline(h=0:nrow(boot.out.plot),col="light gray")
abline(h=(1:nrow(boot.out.plot)-.5),col="light gray",lty=2)
points(asn(boot.out.plot[,2]),(1:nrow(boot.out.plot)-.5),pch="|",cex=2)
####################################################
#names(X) = c("Intercept","Decade Start","Population","Median Household Income","Physical Activity Facilities","Change in Population",
# "Change in Median Household Income","Tract Size")
aa = coef(fit_lambda,lambdaIndex = lamda_indice)[1]$mainEffects$cat
aa.vals = colnames(X)[aa+1]
bb = coef(fit_lambda,lambdaIndex = lamda_indice)[2]$mainEffectsCoef$cat
bb.list = unlist(lapply(bb,function(x) round(exp(x),3)))
cat.main.table = cbind(c("1990-2000","2000-2010"),c(bb.list))
a = coef(fit_lambda,lambdaIndex = lamda_indice)[1]$mainEffects$cont
a.vals = colnames(X)[a+length(aa+1)]
b = coef(fit_lambda,lambdaIndex = lamda_indice)[2]$mainEffectsCoef$cont
b.list = lapply(b,function(x) round(exp(x),3))
#main effects
main.table = as.data.frame(cbind(a.vals,b.list))
colnames(main.table) = c('var','ests')
main.table$var= as.character(main.table$var)
main.table= main.table[order(main.table$var),]
main.table = as.matrix(main.table)
outcome_fill = mat.or.vec((length(coef(fit_lambda,lambdaIndex = lamda_indice)$interactions$contcont[,1])),1)
for (i in 1:length(outcome_fill)){
c = coef(fit_lambda,lambdaIndex = lamda_indice)$interactions$contcont[i,1]
d = coef(fit_lambda,lambdaIndex = lamda_indice)$interactions$contcont[i,2]
outcome_fill[i] = paste(colnames(X)[c+2],colnames(X)[d+2],sep=" X ")
}
inter_values = unlist(coef(fit_lambda,lambdaIndex = lamda_indice)$interactionsCoef$contcont)
inter_values = lapply(inter_values,function(x) round(exp(x),3))
inter.table = cbind(outcome_fill,inter_values)
outcome_fill2a = mat.or.vec((length(coef(fit_lambda,lambdaIndex = lamda_indice)$interactions$catcont[,1])),1)
for (i in 1:length(outcome_fill2a)){
c = coef(fit_lambda,lambdaIndex = lamda_indice)$interactions$catcont[i,1]
d = coef(fit_lambda,lambdaIndex = lamda_indice)$interactions$catcont[i,2]+length(aa)+1
outcome_fill2a[i] = paste('1990-2000','x',colnames(X)[d],sep=" ") #names(X)[c],'
}
outcome_fill2b = mat.or.vec((length(coef(fit_lambda,lambdaIndex = lamda_indice)$interactions$catcont[,1])),1)
for (i in 1:length(outcome_fill2b)){
c = coef(fit_lambda,lambdaIndex = lamda_indice)$interactions$catcont[i,1]
d = coef(fit_lambda,lambdaIndex = lamda_indice)$interactions$catcont[i,2]+length(aa)+1
outcome_fill2b[i] = paste('2000-2010','x',colnames(X)[d],sep=" ")#names(X)[c],
}
outcome_fill2 = rbind(t(t(outcome_fill2a)),t(t(outcome_fill2b)))
cat.inter_values = unlist(coef(fit_lambda,lambdaIndex = lamda_indice)$interactionsCoef$catcont)
cat.inter_values = lapply(cat.inter_values,function(x) round(exp(x),3))
cat.inter_values.transfer = c(cat.inter_values[1],cat.inter_values[2])
cat.inter.table = cbind(outcome_fill2,cat.inter_values.transfer)
cat.inter.table = rbind(cat.inter.table[1:5,],c('',''),cat.inter.table[6:10,])
rbind(cat.main.table,main.table,inter.table,cat.inter.table)
output_table = rbind(cat.main.table,main.table,c('',''),c('Interaction Effects',''),inter.table,c('',''),cat.inter.table)
colnames(output_table) = c("Parameter", "OR")
output_table = rbind(output_table,cbind(c('','Pseudo R^2'),c('',pR2)))
output_table = rbind(output_table,cbind(c('','R^2'),c('',r2)))
output_table[,2] = log(as.numeric(output_table[,2]))
output_table
#write.csv(output_table,"physact_inter_020515.csv") #physact_inter_020515.csv
inter_plot = output_table[c(1:18,21:36,38,39),c(1,2)]
inter_plot[,2] = as.numeric(inter_plot[,2])
par(oma=c(0,12,0,0))
barplot2(asn(inter_plot[,2]), plot.ci= FALSE, horiz= TRUE,names.arg=inter_plot[,1],xlim=c(-.04,.67) ,
col="white",border="white",las=1,space=0, main="Coefficient Estimates for Total Commercial Physical Activity Facilities through Interactions")
abline(v=0)
abline(h=0:36,col="light gray")
abline(h=(1:36-.5),col="light gray",lty=2)
points((inter_plot[,2]),(1:36-.5),pch="|")
##### adjust for dichotmous
aa = coef(fit_lambda2,lambdaIndex = lamda_indice2)[1]$mainEffects$cat
aa.vals = colnames(X)[aa+1]
bb = coef(fit_lambda2,lambdaIndex = lamda_indice2)[2]$mainEffectsCoef$cat
bb.list = unlist(lapply(bb,function(x) round(exp(x),3)))
cat.main.table = cbind(c("1990-2000","2000-2010"),c(bb.list))
a = coef(fit_lambda2,lambdaIndex = lamda_indice2)[1]$mainEffects$cont
a.vals = colnames(X)[a+length(aa+1)]
b = coef(fit_lambda2,lambdaIndex = lamda_indice2)[2]$mainEffectsCoef$cont
b.list = lapply(b,function(x) round(exp(x),3))
#main effects
main.table = as.data.frame(cbind(a.vals,b.list))
colnames(main.table) = c('var','ests')
main.table$var= as.character(main.table$var)
main.table= main.table[order(main.table$var),]
main.table = as.matrix(main.table)
outcome_fill = mat.or.vec((length(coef(fit_lambda2,lambdaIndex = lamda_indice2)$interactions$contcont[,1])),1)
for (i in 1:length(outcome_fill)){
c = coef(fit_lambda2,lambdaIndex = lamda_indice2)$interactions$contcont[i,1]
d = coef(fit_lambda2,lambdaIndex = lamda_indice2)$interactions$contcont[i,2]
outcome_fill[i] = paste(colnames(X)[c+2],colnames(X)[d+2],sep=" X ")
}
inter_values = unlist(coef(fit_lambda2,lambdaIndex = lamda_indice2)$interactionsCoef$contcont)
inter_values = lapply(inter_values,function(x) round(exp(x),3))
inter.table = cbind(outcome_fill,inter_values)
outcome_fill2a = mat.or.vec((length(coef(fit_lambda2,lambdaIndex = lamda_indice2)$interactions$catcont[,1])),1)
for (i in 1:length(outcome_fill2a)){
c = coef(fit_lambda2,lambdaIndex = lamda_indice2)$interactions$catcont[i,1]
d = coef(fit_lambda2,lambdaIndex = lamda_indice2)$interactions$catcont[i,2]+length(aa)+1
outcome_fill2a[i] = paste('1990-2000','x',colnames(X)[d],sep=" ") #names(X)[c],'
}
outcome_fill2b = mat.or.vec((length(coef(fit_lambda2,lambdaIndex = lamda_indice2)$interactions$catcont[,1])),1)
for (i in 1:length(outcome_fill2b)){
c = coef(fit_lambda2,lambdaIndex = lamda_indice2)$interactions$catcont[i,1]
d = coef(fit_lambda2,lambdaIndex = lamda_indice2)$interactions$catcont[i,2]+length(aa)+1
outcome_fill2b[i] = paste('2000-2010','x',colnames(X)[d],sep=" ")#names(X)[c],
}
outcome_fill2 = rbind(t(t(outcome_fill2a)),t(t(outcome_fill2b)))
cat.inter_values = unlist(coef(fit_lambda2,lambdaIndex = lamda_indice2)$interactionsCoef$catcont)
cat.inter_values = lapply(cat.inter_values,function(x) round(exp(x),3))
cat.inter_values.transfer = c(cat.inter_values[1],cat.inter_values[2])
cat.inter.table = cbind(outcome_fill2,cat.inter_values.transfer)
cat.inter.table = rbind(cat.inter.table[1:5,],c('',''),cat.inter.table[6:10,])
rbind(cat.main.table,main.table,inter.table,cat.inter.table)
output_table = rbind(cat.main.table,main.table,c('',''),c('Interaction Effects',''),inter.table,c('',''),cat.inter.table)
colnames(output_table) = c("Parameter", "OR")
output_table = rbind(output_table,cbind(c('','Pseudo R^2'),c('',pR2)))
output_table = rbind(output_table,cbind(c('','R^2'),c('',r2)))
#write.csv(output_table,"physact_inter_021015_binom.csv") #physact_inter_020515.csv
output_table[,2] = round(as.numeric(output_table[,2]),3)
########## trying to quicken/automate subsetting of names & readable changes
output_table = output_table[!output_table[,1]%in%c("","Interaction Effects","R^2"),]
non.int = grep("X|x", output_table[,1],invert=TRUE,value=TRUE)
yes.int = grep("X|x", output_table[,1],value=TRUE)
demv = c("1990-2000|2000-2010|hinc|pasian|phisp|pnhblk|pop|ppov")
envv = c("nFlagged|phys_act_b0|tr_size|adj_tracts|no_adj")
main.dem = grep(demv,non.int,value=TRUE)
main.env =grep(envv,non.int,value=TRUE)
int.env.only = grep(demv,yes.int,invert=TRUE,value=TRUE)
int.dem.only =grep(envv,yes.int,invert=TRUE,value=TRUE)
int.only.either = c(int.env.only,int.dem.only)
int.both.env.dem = grep(paste(int.only.either,collapse="|"),yes.int,invert=TRUE,value=TRUE)
dem.out = output_table[output_table[,1]%in%c(int.dem.only,main.dem),]
env.out = output_table[output_table[,1]%in%c(int.env.only,main.env),]
both.out = output_table[output_table[,1]%in%c(int.both.env.dem),]
pre = c("X","\\.adj","\\.c","hinc","\\.","nFlagged","\\_pa","pasian","phisp","phys_act_b0","pnhblk",
"pop","ppov","tr_size","adj_tracts","no_adj"," adj"," c", "pnhblk")
post= c("x"," Adjacent"," Change","Median Income"," ","Other Businesses","","% Asian","% Hispanic","Previous Physical Activity Facilities",
"% Non-Hispanic Black","Population","% Poverty","Tract Area","Total Adjacent Tracts","No Adjacent Tracts"," Adjacent"," Change","% Non-Hispanic Black")
for (i in 1:length(pre)){
dem.out[,1] = sub(pre[i],post[i],dem.out[,1])
env.out[,1] = sub(pre[i],post[i],env.out[,1])
both.out[,1] = sub(pre[i],post[i],both.out[,1])
}
dem.out = dem.out[nrow(dem.out):1,]
env.out = env.out[nrow(env.out):1,]
both.out = both.out[nrow(both.out):1,]
#inter_plot = output_table[output_table[,2]!=1,]
#inter_plot[,2] = as.numeric(inter_plot[,2])
tables = list(dem.out,env.out,both.out)
i=1
i=2
i=3
#for (i in 1:length(tables)){
inter_plot = as.data.frame(tables[i])
par(oma=c(0,25,0,0))
xmin = min(as.numeric(inter_plot[,2]))-.005
xmax = max(as.numeric(inter_plot[,2]))+.005
barplot2(asn(inter_plot[,2]), plot.ci= FALSE, horiz= TRUE,names.arg=inter_plot[,1],xlim=c(xmin,xmax) ,
col="white",border="white",las=1,space=0)#, main="OR Estimates for Increase in Commercial Physical Activity Facilities through Interactions")
abline(v=1)
abline(h=0:nrow(inter_plot),col="light gray")
abline(h=(1:nrow(inter_plot)-.5),col="light gray",lty=2)
points((inter_plot[,2]),(1:nrow(inter_plot)-.5),pch="|",cex=2)
#}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.