##Publication Table Formulas - China

##Load Packages and Data
library(Deducer)
library(psych)
library(catspec)
library(gdata)
library(Hmisc)


library(xtable)
library(psych)
x <- read.csv(file = "/Users/erinschneider/Desktop/IDFDataFiles/China/CHINA_PERSON_DATA.csv", sep = ",")



splice3 <- function (x, y, w,
                    colnames.splice = NULL,
                    pc = TRUE) {
  z <- data.frame(matrix(rbind(x, y, w), nrow(x)))
  rownames(z) <- rownames(x)
  if (ncol(z) / 2 == length(colnames.splice)) colnames.splice <- rep(colnames.splice, 2)
  colnames(z) <- colnames.splice
  z
}

splice <- function (x, y,
                    colnames.splice=NULL,
                    pc=TRUE) {
  z<-data.frame(matrix(rbind(x,y),nrow(x)))
  rownames(z)<-rownames(x)
  if (pc) colnames.splice <- matrix(rbind(colnames(x), rep("%", ncol(x))), nrow=1)
  if (is.null(colnames.splice)) {
    colnames.splice <- matrix(rbind(paste(deparse(substitute(x)), ".", colnames(x), sep=""),
                                    paste(deparse(substitute(y)), ".", colnames(y), sep="")), nrow=1)
  }
  colnames(z) <- colnames.splice
  z
}




ci <- function(x){
    mean <- mean(x, na.rm=TRUE)
    s <- sd(x, na.rm=TRUE)
    n <- length(na.omit(x))
    error <- qnorm(0.975)*s/sqrt(n)
    left <- mean-error
    right <- mean+error
    conf <- c(left, right)
    conf
}
yesno <- function (x) {
    if (is.na(x)) {
        y <- NA
    } else {
        if (x == "Yes") {
            y <- 1
        } else {
            y <- 0
        }
        y
    }
}


y <- x[x$case == "Diabetes" | x$case == "No Diabetes",]
y$case <- factor(y$case, levels = c("Diabetes", "No Diabetes"))
dm <- y[y$case == "Diabetes",]
nd <- y[y$case == "No Diabetes",]


##China Sites
bj <- y[y$site == "Beijing",]
sy <- y[y$site == "liaoning",]
sn <- y[y$site == "Shandong",]
sx <- y[y$site == "Shanxi",]
sg <- y[y$site == "Shanghai",]
hu <- y[y$site == "Hunan",]
fu <- y[y$site == "Fujian",]
si <- y[y$site == "Sichuan",]
sz <- y[y$site == "Shaanxi",]
xi <- y[y$site == "Xinjiang",]


##Location Ns
fm <- table(y$site, y$location)
fm

fm <- table(y$location)
fm

#Mean Age & s.d.
describe.by(bj$age, bj$location)
describe.by(bj$age)
describe.by(sy$age, sy$location)
describe.by(sy$age)
describe.by(sn$age, sn$location)
describe.by(sn$age)
describe.by(sx$age, sx$location)
describe.by(sx$age)
describe.by(sg$age, sg$location)
describe.by(sg$age)
describe.by(hu$age, hu$location)
describe.by(hu$age)
describe.by(fu$age, fu$location)
describe.by(fu$age)
describe.by(si$age, si$location)
describe.by(si$age)
describe.by(sz$age, sz$location)
describe.by(sz$age)
describe.by(xi$age, xi$location)
describe.by(xi$age)
describe.by(y$age, y$location)
describe.by(y$age)

#Percent Female
fm <- table(bj$location, bj$sex)
p.fm <- prop.table(fm, 1)*100
p.fm <- p.fm[, 1]
p.fm

fm <- table(sy$location, sy$sex)
p.fm <- prop.table(fm, 1)*100
p.fm <- p.fm[, 1]
p.fm

fm <- table(sn$location, sn$sex)
p.fm <- prop.table(fm, 1)*100
p.fm <- p.fm[, 1]
p.fm

fm <- table(sx$location, sx$sex)
p.fm <- prop.table(fm, 1)*100
p.fm <- p.fm[, 1]
p.fm

fm <- table(sg$location, sg$sex)
p.fm <- prop.table(fm, 1)*100
p.fm <- p.fm[, 1]
p.fm

fm <- table(hu$location, hu$sex)
p.fm <- prop.table(fm, 1)*100
p.fm <- p.fm[, 1]
p.fm

fm <- table(fu$location, fu$sex)
p.fm <- prop.table(fm, 1)*100
p.fm <- p.fm[, 1]
p.fm

fm <- table(si$location, si$sex)
p.fm <- prop.table(fm, 1)*100
p.fm <- p.fm[, 1]
p.fm

fm <- table(sz$location, sz$sex)
p.fm <- prop.table(fm, 1)*100
p.fm <- p.fm[, 1]
p.fm

fm <- table(xi$location, xi$sex)
p.fm <- prop.table(fm, 1)*100
p.fm <- p.fm[, 1]
p.fm

fm <- table(y$location, y$sex)
p.fm <- prop.table(fm, 1)*100
p.fm <- p.fm[, 1]
p.fm

#Percent Female TOTAL 
fm <- table(y$sc, y$sex)
p.fm <- prop.table(fm, 1)*100
p.fm <- p.fm[, 1]
p.fm


#Percent Urban
rr <- table(bj$site, bj$location)
p.rr <- prop.table(rr, 1)*100
p.rr <- p.rr[, 2]
p.rr

rr <- table(sy$site, sy$location)
p.rr <- prop.table(rr, 1)*100
p.rr <- p.rr[, 2]
p.rr

rr <- table(sn$site, sn$location)
p.rr <- prop.table(rr, 1)*100
p.rr <- p.rr[, 2]
p.rr

rr <- table(sx$site, sx$location)
p.rr <- prop.table(rr, 1)*100
p.rr <- p.rr[, 2]
p.rr

rr <- table(sg$site, sg$location)
p.rr <- prop.table(rr, 1)*100
p.rr <- p.rr[, 2]
p.rr

rr <- table(hu$site, hu$location)
p.rr <- prop.table(rr, 1)*100
p.rr <- p.rr[, 2]
p.rr

rr <- table(fu$site, fu$location)
p.rr <- prop.table(rr, 1)*100
p.rr <- p.rr[, 2]
p.rr

rr <- table(si$site, si$location)
p.rr <- prop.table(rr, 1)*100
p.rr <- p.rr[, 2]
p.rr

rr <- table(sz$site, sz$location)
p.rr <- prop.table(rr, 1)*100
p.rr <- p.rr[, 2]
p.rr

rr <- table(xi$site, xi$location)
p.rr <- prop.table(rr, 1)*100
p.rr <- p.rr[, 2]
p.rr

rr <- table(y$site, y$location)
p.rr <- prop.table(rr, 1)*100
p.rr <- p.rr[, 2]
p.rr

rr <- table(y$location)
p.rr <- prop.table(rr)*100
p.rr

#Percent Urban TOTALS by Site
rr <- table(y$site, y$location)
p.rr <- prop.table(rr, 1)*100
p.rr <- p.rr[, 2]
p.rr


#Percent Complication by Location - NGTS
com <- table(nd$ncd4, nd$location)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(nd$ncd1, nd$location)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(nd$ncd17, nd$location)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(nd$macro, nd$location)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(nd$micro, nd$location)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(nd$acute, nd$location)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(nd$anyncd, nd$location)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

#Percent Complications by age - NGTS
com <- table(nd$ncd4, nd$age.b)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(nd$ncd1, nd$age.b)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(nd$ncd17, nd$age.b)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(nd$macro, nd$age.b)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(nd$micro, nd$age.b)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(nd$acute, nd$age.b)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(nd$anyncd, nd$age.b)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com


##TOTALS - NGT complications
com <- table(nd$ncd4, nd$sc)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(nd$ncd1, nd$sc)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(nd$ncd17, nd$sc)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(nd$macro, nd$sc)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(nd$micro, nd$sc)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(nd$acute, nd$sc)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(nd$anyncd, nd$sc)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com


#Percent Complication by Location - DMs
com <- table(dm$ncd4, dm$location)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$ncd1, dm$location)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$ncd2, dm$location)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$ncd5, dm$location)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$kidney, dm$location)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$macro, dm$location)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$micro, dm$location)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$acute, dm$location)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$anyncd, dm$location)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

#Percent Complications by age - DMs
com <- table(dm$ncd4, dm$age.b)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$ncd1, dm$age.b)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$ncd2, dm$age.b)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$ncd5, dm$age.b)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$kidney, dm$age.b)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$macro, dm$age.b)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$micro, dm$age.b)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$acute, dm$age.b)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$anyncd, dm$age.b)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

#Percent Complications by DM Duration
com <- table(dm$ncd4, dm$dur.gp)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$ncd1, dm$dur.gp)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$ncd2, dm$dur.gp)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$ncd5, dm$dur.gp)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$kidney, dm$dur.gp)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$macro, dm$dur.gp)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$micro, dm$dur.gp)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$acute, dm$dur.gp)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$anyncd, dm$dur.gp)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

##TOTALS - DM Complications
com <- table(dm$ncd4, dm$sc)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$ncd1, dm$sc)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$ncd17, dm$sc)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$macro, dm$sc)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$micro, dm$sc)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$acute, dm$sc)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$anyncd, dm$sc)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

##Relevel case variable
y$case <- relevel(y$case, ref="No Diabetes")

##Mean Inpatient Admissions Last 90 Days
#Urban/Rural
hosp.ad <- y$hosp.ad
aggregate.table(hosp.ad, y$case, y$location, FUN=mean, na.rm=TRUE)
#Age Bracket
aggregate.table(hosp.ad, y$case, y$age.b, FUN=mean, na.rm=TRUE)
#Duration
aggregate.table(hosp.ad, y$case, y$dur.gp, FUN=mean, na.rm=TRUE)
##P-Values

summary(lm(hosp.ad ~ age + sex + location + dm.yrs, data = dm))
summary(lm(hosp.ad ~ age + sex + location, data = nd))
summary(lm(hosp.ad ~ case + sex + age + location, data=y))

#Mean Annual Inpatient Admissions Per Person
y$niwa[is.na(y$niwa)] <- 0
dur.hosp.yr <- y$hosp.ad * y$niwa * 4
#Urban/Rural
aggregate.table(dur.hosp.yr, y$case, y$location, FUN=mean, na.rm=TRUE)
#Age Bracket
aggregate.table(dur.hosp.yr, y$case, y$age.b, FUN=mean, na.rm=TRUE)
#Duration
aggregate.table(dur.hosp.yr, y$case, y$dur.gp, FUN=mean, na.rm=TRUE)

##Mean Inpatient Length of Stay if Admitted
#Urban/Rural
y$nights.hosp[y$nights.hosp == 0] <- NA
hosp.nights <- y$nights.hosp
#Urban/Rural
aggregate.table(hosp.nights, y$case, y$location, FUN=mean, na.rm=TRUE)
#Age Bracket
aggregate.table(hosp.nights, y$case, y$age.b, FUN=mean, na.rm=TRUE)
#Duration
aggregate.table(hosp.nights, y$case, y$dur.gp, FUN=mean, na.rm=TRUE)
#p-values

lm1 <- summary(lm(nights.hosp ~ age + sex + location + dm.yrs data, = dm))
lm2 <- summary(lm(nights.hosp ~ age + sex + location, data = nd))
lm3 <- summary(lm(nights.hosp ~ case + age + sex + location, data = y))

##Mean Inpatient Expenditures per admission
#Urban/Rural
inpat.exp <- y$tpyad
inpat.exp1 <- dm$tpyad
inpat.exp2 <- nd$tpyad
#Urban/Rural
aggregate.table(inpat.exp, y$case, y$location, FUN=mean, na.rm=TRUE)
#Age Bracket
aggregate.table(inpat.exp, y$case, y$age.b, FUN=mean, na.rm=TRUE)
#Duration
aggregate.table(inpat.exp, y$case, y$dur.gp, FUN=mean, na.rm=TRUE)

summary(lm(inpat.exp1 ~ age + sex + location + dm.yrs, data = dm))
summary(lm(inpat.exp2 ~ age + sex + location, data = nd))
summary(lm(inpat.exp~ case + age + sex + location, data = y))

##Mean Inpatient Expenditures per person
#Urban/Rural

y$tpyad[is.na(y$tpyad)] <- 0
dm$tpyad[is.na(dm$tpyad)] <- 0
nd$tpyad[is.na(nd$tpyad)] <- 0
inpat.exp <- y$tpyad
inpat.exp1 <- dm$tpyad
inpat.exp2 <- nd$tpyad
#Urban/Rural
aggregate.table(inpat.exp, y$case, y$location, FUN=mean, na.rm=TRUE)
#Age Bracket
aggregate.table(inpat.exp, y$case, y$age.b, FUN=mean, na.rm=TRUE)
#Duration
aggregate.table(inpat.exp, y$case, y$dur.gp, FUN=mean, na.rm=TRUE)

summary(lm(inpat.exp1 ~ age + sex + location + dm.yrs, data = dm))
summary(lm(inpat.exp2 ~ age + sex + location, data = nd))
summary(lm(inpat.exp~ case + age + sex + location, data = y))
##Mean Outpatient Admissions Last 90 Days per person
#Urban/Rural
hosp.op <- y$hosp.op
#Urban/Rural
aggregate.table(hosp.op, y$case, y$location, FUN=mean, na.rm=TRUE)
#Age Bracket
aggregate.table(hosp.op, y$case, y$age.b, FUN=mean, na.rm=TRUE)
#Duration
aggregate.table(hosp.op, y$case, y$dur.gp, FUN=mean, na.rm=TRUE)

y$hosp.op1[y$hosp.op > 0] <- 1
y$hosp.op1[y$hosp.op <= 0] <- 0
dm$hosp.op1[dm$hosp.op > 0] <- 1
dm$hosp.op1[dm$hosp.op <= 0] <- 0
nd$hosp.op1[nd$hosp.op > 0] <- 1
nd$hosp.op1[nd$hosp.op <= 0] <- 0

summary(glm(hosp.op1 ~ age + sex + location + dm.yrs, data=dm, family=binomial(link= "logit")))
summary(glm(hosp.op1 ~ age + sex + location, data=nd, family=binomial(link= "logit")))
summary(glm(hosp.op1 ~ case + sex + age + location, data=y, family=binomial(link= "logit")))

##Mean Western Med Visits Last 90 Days
#Urban/Rural
wm.op <- y$wm90
wm.op1 <- dm$wm90
wm.op2 <- nd$wm90
aggregate.table(wm.op, y$case, y$location, FUN=mean, na.rm=TRUE)
#Age Bracket
aggregate.table(wm.op, y$case, y$age.b, FUN=mean, na.rm=TRUE)
#Duration
aggregate.table(wm.op, y$case, y$dur.gp, FUN=mean, na.rm=TRUE)
#Total
aggregate.table(wm.op, y$case, y$sc, FUN=mean, na.rm=TRUE)


summary(lm(wm.op1 ~ age + sex + location + dm.yrs, data = dm))
summary(lm(wm.op2 ~ age + sex + location, data = nd))
summary(lm(wm.op ~ case + sex + age + location, data=y))

##Mean TCM Visits Last 90 Days
#Urban/Rural
th.op <- y$th90
th.op1 <- dm$th90
th.op2 <- nd$th90
aggregate.table(th.op, y$case, y$location, FUN=mean, na.rm=TRUE)
#Age Bracket
aggregate.table(th.op, y$case, y$age.b, FUN=mean, na.rm=TRUE)
#Duration
aggregate.table(th.op, y$case, y$sc, FUN=mean, na.rm=TRUE)
aggregate.table(th.op, y$case, y$dur.gp, FUN=mean, na.rm=TRUE)

summary(lm(th.op1 ~ age + sex + location + dm.yrs, data = dm))
summary(lm(th.op2 ~ age + sex + location, data = nd))
summary(lm(th.op ~ case + sex + age + location, data=y))

##Mean Com. Health Worker Visists Last 90 Days
#Urban/Rural
ch.op <- y$ch90
ch.op1 <- dm$ch90
ch.op2 <- nd$ch90
aggregate.table(ch.op, y$case, y$location, FUN=mean, na.rm=TRUE)
#Age Bracket
aggregate.table(ch.op, y$case, y$age.b, FUN=mean, na.rm=TRUE)
#Duration
aggregate.table(ch.op, y$case, y$dur.gp, FUN=mean, na.rm=TRUE)
aggregate.table(ch.op, y$case, y$sc, FUN=mean, na.rm=TRUE)

summary(lm(ch.op1 ~ age + sex + location + dm.yrs, data = dm))
summary(lm(ch.op2 ~ age + sex + location, data = nd))
summary(lm(ch.op ~ case + sex + age + location, data=y))

##Mean Total Outpatient Medical Visists Last 90 Days
#Urban/Rural
tot.op <- y$hosp.op + y$wm90 + y$th90 + y$ch90
tot.op1 <- dm$hosp.op + dm$wm90 + dm$th90 + dm$ch90
tot.op2 <- nd$hosp.op + nd$wm90 + nd$th90 + nd$ch90
aggregate.table(tot.op, y$case, y$location, FUN=mean, na.rm=TRUE)
#Age Bracket
aggregate.table(tot.op, y$case, y$age.b, FUN=mean, na.rm=TRUE)
#Duration
aggregate.table(tot.op, y$case, y$dur.gp, FUN=mean, na.rm=TRUE)
aggregate.table(tot.op, y$case, y$sc, FUN=mean, na.rm=TRUE)


summary(lm(tot.op1 ~ age + sex + location + dm.yrs, data = dm))
summary(lm(tot.op2 ~ age + sex + location, data = nd))
summary(lm(tot.op ~ case + sex + age + location, data=y))

##Mean Outpatient Expenditures Last 90 Days if Admitted
#Urban/Rural
y$tyop1[is.na(y$tyop1)] <- 0
dm$tyop1[is.na(dm$tyop1)] <- 0
nd$tyop1[is.na(nd$tyop1)] <- 0

y$wmco[is.na(y$wmco)] <- 0
dm$wmco[is.na(dm$wmco)] <- 0
nd$wmco[is.na(nd$wmco)] <- 0

y$thco[is.na(y$thco)] <- 0
dm$thco[is.na(dm$thco)] <- 0
nd$thco[is.na(nd$thco)] <- 0

y$chco[is.na(y$chco)] <- 0
dm$chco[is.na(dm$chco)] <- 0
nd$chco[is.na(nd$chco)] <- 0

op.exp <- y$tyop1 
aggregate.table(op.exp, y$case, y$location, FUN=mean, na.rm=TRUE)
#Age Bracket
aggregate.table(op.exp, y$case, y$age.b, FUN=mean, na.rm=TRUE)
#Duration
aggregate.table(op.exp, y$case, y$dur.gp, FUN=mean, na.rm=TRUE)

summary(glm(hosp.op1 ~ age + sex + location + dm.yrs, data=dm, family=binomial(link= "logit")))
summary(glm(hosp.op1 ~ age + sex + location, data=nd, family=binomial(link= "logit")))
summary(glm(op.exp ~ case + sex + age + location, data=y, family=binomial(link= "logit")))

##Mean TOTAL Outpatient Expenditures Last 90 Days Per Visit
#Urban/Rural

y$tyop1[is.na(y$tyop1)] <- 0
op.exp <- y$tyop1
aggregate.table(op.exp, y$case, y$location, FUN=median, na.rm=TRUE)
#Age Bracket
aggregate.table(op.exp, y$case, y$age.b, FUN=median, na.rm=TRUE)
#Duration
aggregate.table(op.exp, y$case, y$dur.gp, FUN=median, na.rm=TRUE)

y$wmco[y$wmco > 1999] <- NA
wm.op <- y$wmco
wm.op1 <- dm$wmco
wm.op2 <- nd$wmco
aggregate.table(wm.op, y$case, y$location, FUN=mean, na.rm=TRUE)
#Age Bracket
aggregate.table(wm.op, y$case, y$age.b, FUN=mean, na.rm=TRUE)
#Duration
aggregate.table(wm.op, y$case, y$dur.gp, FUN=mean, na.rm=TRUE)
#Total
aggregate.table(wm.op, y$case, y$sc, FUN=mean, na.rm=TRUE)


summary(lm(wm.op1 ~ age + sex + location + dm.yrs, data = dm))
summary(lm(wm.op2 ~ age + sex + location, data = nd))
summary(lm(wm.op ~ case + sex + age + location, data=y))

##Mean PH Expenditures Last 90 Days
#Urban/Rural
ph.op <- y$phco
ph.op1 <- dm$phco
ph.op2 <- nd$phco
aggregate.table(ph.op, y$case, y$location, FUN=mean, na.rm=TRUE)
#Age Bracket
aggregate.table(ph.op, y$case, y$age.b, FUN=mean, na.rm=TRUE)
#Duration
aggregate.table(ph.op, y$case, y$sc, FUN=mean, na.rm=TRUE)
aggregate.table(ph.op, y$case, y$dur.gp, FUN=mean, na.rm=TRUE)

summary(lm(ph.op1 ~ age + sex + location + dm.yrs, data = dm))
summary(lm(ph.op2 ~ age + sex + location, data = nd))
summary(lm(ph.op ~ case + sex + age + location, data=y))


##Mean Traditional Healer Expenditures Last 90 Days
#Urban/Rural
th.op <- y$thco
th.op1 <- dm$thco
th.op2 <- nd$thco
aggregate.table(th.op, y$case, y$location, FUN=mean, na.rm=TRUE)
#Age Bracket
aggregate.table(th.op, y$case, y$age.b, FUN=mean, na.rm=TRUE)
#Duration
aggregate.table(th.op, y$case, y$sc, FUN=mean, na.rm=TRUE)
aggregate.table(th.op, y$case, y$dur.gp, FUN=mean, na.rm=TRUE)

summary(lm(th.op1 ~ age + sex + location + dm.yrs, data = dm))
summary(lm(th.op2 ~ age + sex + location, data = nd))
summary(lm(th.op ~ case + sex + age + location, data=y))


##Mean Com. Health Worker Visists Last 90 Days
#Urban/Rural
ch.op <- y$chc
ch.op1 <- dm$chc
ch.op2 <- nd$chco
aggregate.table(ch.op, y$case, y$location, FUN=mean, na.rm=TRUE)
#Age Bracket
aggregate.table(ch.op, y$case, y$age.b, FUN=mean, na.rm=TRUE)
#Duration
aggregate.table(ch.op, y$case, y$dur.gp, FUN=mean, na.rm=TRUE)
aggregate.table(ch.op, y$case, y$sc, FUN=mean, na.rm=TRUE)

summary(lm(ch.op1 ~ age + sex + location + dm.yrs, data = dm))
summary(lm(ch.op2 ~ age + sex + location, data = nd))
summary(lm(ch.op ~ case + sex + age + location, data=y))



##TOTALS
describe.by(y$hosp.ad, y$case)
dur.hosp.yr <- y$dur.hosp.yr
describe.by(y$dur.hosp.yr, y$case)
describe.by(y$tpyad.ad.yr, y$case)
describe.by(y$nights.hosp, y$case)
describe.by(y$tpyad, y$case)
describe.by(y$hosp.op, y$case)
describe.by(y$tyop1, y$case)

##Mean Drugs Taken
#Urban/Rural
any.drug <- y$any.drug
aggregate.table(any.drug, y$case, y$location, FUN=mean, na.rm=TRUE)
#Age Bracket
aggregate.table(any.drug, y$case, y$age.b, FUN=mean, na.rm=TRUE)
#Duration
aggregate.table(any.drug, y$case, y$dur.gp, FUN=mean, na.rm=TRUE)

###MEDICATIONS TABLE
#Percent Medications by Location
com <- table(dm$met, dm$location)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$su, dm$location)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$ins, dm$location)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$acar, dm$location)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$other.gllow, dm$location)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$statin, dm$location)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$other.lipid, dm$location)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$diur, dm$location)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$ace, dm$location)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$arb, dm$location)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$bb, dm$location)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$ccb, dm$location)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$other.bp, dm$location)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$asp, dm$location)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$no.asp.anticoag, dm$location)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$other.nondm, dm$location)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$any.herb, dm$location)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com



#Percent Complications by age
com <- table(dm$met, dm$age.b)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$su, dm$age.b)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$ins, dm$age.b)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$acar, dm$age,b)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$other.gllow, dm$age.b)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$statin, dm$age.b)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$other.lipid, dm$age.b)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$diur, dm$age.b)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$ace, dm$age.b)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$arb, dm$age.b)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$bb, dm$age.b)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$ccb, dm$age.b)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$other.bp, dm$age.b)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$asp, dm$age.b)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$no.asp.anticoag, dm$age.b)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$other.nondm, dm$age.b)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$any.herb, dm$age.b)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

#Percent Complications by DM Duration
com <- table(dm$met, dm$dur.gp)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$su, dm$dur.gp)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$ins, dm$dur.gp)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$acar, dm$dur.gp)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$other.gllow, dm$dur.gp)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$statin, dm$dur.gp)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$other.lipid, dm$dur.gp)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$diur, dm$dur.gp)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$ace, dm$dur.gp)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$arb, dm$dur.gp)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$bb, dm$dur.gp)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$ccb, dm$dur.gp)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$other.bp, dm$dur.gp)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$asp, dm$dur.gp)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$no.asp.anticoag, dm$dur.gp)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$other.nondm, dm$dur.gp)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$any.herb, dm$dur.gp)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

##Meds by history of HTN
com <- table(dm$met, dm$ncd5)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$su, dm$ncd5)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$ins, dm$ncd5)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$acar, dm$ncd5)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$other.gllow, dm$ncd5)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$statin, dm$ncd5)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$other.lipid, dm$ncd5)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$diur, dm$ncd5)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$ace, dm$ncd5)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$arb, dm$ncd5)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$bb, dm$ncd5)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$ccb, dm$ncd5)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$other.bp, dm$ncd5)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$asp, dm$ncd5)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$no.asp.anticoag, dm$ncd5)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$other.nondm, dm$ncd5)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$any.herb, dm$ncd5)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

##TOTALS for drugs taken
com <- table(dm$met, dm$sc)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$su, dm$sc)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$ins, dm$sc)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$other.gllow, dm$sc)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$statin, dm$sc)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$other.lipid, dm$sc)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$diur, dm$sc)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$ace, dm$sc)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$arb, dm$sc)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$bb, dm$sc)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$ccb, dm$sc)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$other.bp, dm$sc)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$asp, dm$sc)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$no.asp.anticoag, dm$sc)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$other.nondm, dm$sc)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$any.herb, dm$sc)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com



##P-Values for drugs taken
dm$any.gllow1[dm$any.gllow == "Yes"] <- 1
dm$any.gllow1[dm$any.gllow == "No"] <- 0
summary(glm(any.gllow2 ~ age + sex + location + dm.yrs, data=dm, family=binomial(link= "logit")))

dm$any.met1[dm$met== "Yes"] <- 1
dm$any.met1[dm$met == "No"] <- 0
summary(glm(any.met1 ~ age + sex + location + dm.yrs, data=dm, family=binomial(link= "logit")))


dm$su1[dm$su == "Yes"] <- 1
dm$su1[dm$su == "No"] <- 0
summary(glm(su1 ~ age + sex + location + dm.yrs, data=dm, family=binomial(link= "logit")))


dm$ins1[dm$ins == "Yes"] <- 1
dm$ins1[dm$ins == "No"] <- 0
summary(glm(ins1 ~ age + sex + location + dm.yrs, data=dm, family=binomial(link= "logit")))


dm$acar1[dm$acar == "Yes"] <- 1
dm$acar1[dm$acar == "No"] <- 0
summary(glm(acar1 ~ age + sex + location + dm.yrs, data=dm, family=binomial(link= "logit")))

dm$other.gllow1[dm$other.gllow == "Yes"] <- 1
dm$other.gllow1[dm$other.gllow == "No"] <- 0
summary(glm(other.gllow1 ~ age + sex + location + dm.yrs, data=dm, family=binomial(link= "logit")))


dm$any.lipid1[dm$any.lipid == "Yes"] <- 1
dm$any.lipid1[dm$any.lipid == "No"] <- 0
summary(glm(any.lipid1 ~ age + sex + location + dm.yrs, data=dm, family=binomial(link= "logit")))

dm$statin1[dm$statin == "Yes"] <- 1
dm$statin1[dm$statin== "No"] <- 0
summary(glm(statin1 ~ age + sex + location + dm.yrs, data=dm, family=binomial(link= "logit")))

dm$other.lipid1[dm$other.lipid== "Yes"] <- 1
dm$other.lipid1[dm$other.lipid == "No"] <- 0
summary(glm(other.lipid1 ~ age + sex + location + dm.yrs, data=dm, family=binomial(link= "logit")))


dm$diur1[dm$diur == "Yes"] <- 1
dm$diur1[dm$diur == "No"] <- 0
summary(glm(diur1 ~ age + sex + location + dm.yrs, data=dm, family=binomial(link= "logit")))


dm$any.bp1[dm$any.bp == "Yes"] <- 1
dm$any.bp1[dm$any.bp == "No"] <- 0
summary(glm(any.bp1 ~ age + sex + location + dm.yrs, data=dm, family=binomial(link= "logit")))


dm$ace1[dm$ace == "Yes"] <- 1
dm$ace1[dm$ace == "No"] <- 0
summary(glm(ace1 ~ age + sex + location + dm.yrs, data=dm, family=binomial(link= "logit")))

dm$arb1[dm$arb == "Yes"] <- 1
dm$arb1[dm$arb == "No"] <- 0
summary(glm(arb1 ~ age + sex + location + dm.yrs, data=dm, family=binomial(link= "logit")))


dm$bb1[dm$bb== "Yes"] <- 1
dm$bb1[dm$bb == "No"] <- 0
summary(glm(bb1 ~ age + sex + location + dm.yrs, data=dm, family=binomial(link= "logit")))

dm$ccb1[dm$ccb== "Yes"] <- 1
dm$ccb1[dm$ccb == "No"] <- 0
summary(glm(ccb1 ~ age + sex + location + dm.yrs, data=dm, family=binomial(link= "logit")))

dm$other.bp1[dm$other.bp == "Yes"] <- 1
dm$other.bp1[dm$other.bp == "No"] <- 0
summary(glm(other.bp1 ~ age + sex + location + dm.yrs, data=dm, family=binomial(link= "logit")))

dm$asp1[dm$asp == "Yes"] <- 1
dm$asp1[dm$asp == "No"] <- 0
summary(glm(asp1 ~ age + sex + location + dm.yrs, data=dm, family=binomial(link= "logit")))

dm$no.asp.anticoag1[dm$no.asp.anticoag == "Yes"] <- 1
dm$no.asp.anticoag1[dm$no.asp.anticoag == "No"] <- 0
summary(glm(no.asp.anticoag1 ~ age + sex + location + dm.yrs, data=dm, family=binomial(link= "logit")))

dm$any.anticoag1[dm$any.anticoag == "Yes"] <- 1
dm$any.anticoag1[dm$any.anticoag == "No"] <- 0
summary(glm(any.anticoag1 ~ age + sex + location + dm.yrs, data=dm, family=binomial(link= "logit")))

dm$anal1[dm$anal == "Yes"] <- 1
dm$anal1[dm$anal == "No"] <- 0
summary(glm(anal1 ~ age + sex + location + dm.yrs, data=dm, family=binomial(link= "logit")))

dm$other.nondm1[dm$other.nondm== "Yes"] <- 1
dm$other.nondm1[dm$other.nondm == "No"] <- 0
summary(glm(other.nondm1 ~ age + sex + location + dm.yrs, data=dm, family=binomial(link= "logit")))

dm$any.herb1[dm$any.herb== "Yes"] <- 1
dm$any.herb1[dm$any.herb == "No"] <- 0
summary(glm(any.herb1 ~ age + sex + location + dm.yrs, data=dm, family=binomial(link= "logit")))


##BG monitoring by Location
dm$any.hosp1[dm$any.hosp > 0] <- "Yes"
dm$any.hosp1[dm$any.hosp <= 0 | is.na(dm$any.hosp)] <- "No"
com <- table(dm$any.hosp1, dm$location)
p.com <- prop.table(com, 2)*100
any.d <- splice(com, p.com)
any.d

dm$any.hosp2[dm$any.hosp1== "Yes"] <- 1
dm$any.hosp2[dm$any.hosp1 == "No"] <- 0
summary(glm(any.hosp2 ~ age + sex + location + dm.yrs, data=dm, family=binomial(link= "logit")))


com <- table(dm$hp.bg.any, dm$location)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

dm$hp.bg.any1[dm$hp.bg.any == "Yes"] <- 1
dm$hp.bg.any1[dm$hp.bg.any == "No"] <- 0
summary(glm(hp.bg.any1 ~ age + sex + location + dm.yrs, data=dm, family=binomial(link= "logit")))

com <- table(dm$ed.eye.exam, dm$location)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

dm$ed.eye.exam1[dm$ed.eye.exam == "Yes"] <- 1
dm$ed.eye.exam1[dm$ed.eye.exam == "No"] <- 0
summary(glm(ed.eye.exam1 ~ age + sex + location + dm.yrs, data=dm, family=binomial(link= "logit")))


com <- table(dm$ed.foot, dm$location)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

dm$ed.foot1[dm$ed.foot == "Yes"] <- 1
dm$ed.foot1[dm$ed.foot == "No"] <- 0
summary(glm(ed.foot1 ~ age + sex + location + dm.yrs, data=dm, family=binomial(link= "logit")))

com <- table(dm$hm.less.wk, dm$location)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

dm$hm.less.wk1[dm$hm.less.wk == "Yes"] <- 1
dm$hm.less.wk1[dm$hm.less.wk == "No"] <- 0
summary(glm(hm.less.wk1 ~ age + sex + location + dm.yrs, data=dm, family=binomial(link= "logit")))

com <- table(dm$hm.less.day, dm$location)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

dm$hm.less.day1[dm$hm.less.day == "Yes"] <- 1
dm$hm.less.day1[dm$hm.less.day == "No"] <- 0
summary(glm(hm.less.day1 ~ age + sex + location + dm.yrs, data=dm, family=binomial(link= "logit")))

com <- table(dm$hm.day, dm$location)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

dm$hm.day1[dm$hm.day == "Yes"] <- 1
dm$hm.day1[dm$hm.day == "No"] <- 0
summary(glm(hm.day1 ~ age + sex + location + dm.yrs, data=dm, family=binomial(link= "logit")))




##BG monitoring by Age
com <- table(dm$any.hosp1, dm$age.b)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$hp.bg.any, dm$age.b)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$bp.test, dm$age.b)
p.com <- prop.table(com, 2)*100
p.com <- p.com[3, ]
p.com

com <- table(dm$ed.eye.exam, dm$age.b)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$ed.foot, dm$age.b)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$hm.less.wk, dm$age.b)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$hm.less.day, dm$age.b)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$hm.day, dm$age.b)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$hm.any, dm$age.b)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com


##BG monitoring by Duration
com <- table(dm$any.hosp1, dm$dur.gp)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$hp.bg.any, dm$dur.gp)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$ed.eye.exam, dm$dur.gp)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$ed.foot, dm$dur.gp)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$hm.less.wk, dm$dur.gp)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$hm.less.day, dm$dur.gp)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$hm.day, dm$dur.gp)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$hm.any, dm$dur.gp)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

##BG monitoring Totals
com <- table(dm$any.hosp, dm$sc)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$hp.bg.any, dm$sc)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$ed.eye.exam, dm$sc)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$ed.foot, dm$sc)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$hm.less.wk, dm$sc)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$hm.less.day, dm$sc)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$hm.day, dm$sc)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$hm.any, dm$sc)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com


##Length of hospital stay
dm$nights.hosp[dm$nights.hosp == 0] <- NA
hosp.nights <- dm$nights.hosp
aggregate.table(hosp.nights, dm$case, dm$reason.ad, FUN=mean, na.rm=TRUE)

##Payments
##Mean Inpatient Expenditures per admission
#Urban/Rural
inpat.exp <- y$tpyad
aggregate.table(inpat.exp, y$case, y$reason.ad, FUN=mean, na.rm=TRUE)

opat.exp <- y$tyop1
aggregate.table(opat.exp, y$case, y$reason.ad, FUN=mean, na.rm=TRUE)
##Mean number of meds currently taking

##Table 1a
##Ns
table(y$case)
table(y$Outcome)
## % Female
fm <- table(y$case, y$sex)
p.fm <- prop.table(fm, 1)*100
p.fm <- p.fm[, 1]
p.fm

fm <- table(y$Outcome, y$gender)
p.fm <- prop.table(fm, 1)*100
p.fm <- p.fm[, 2]
p.fm
##Complications
com <- table(y$ncd4, y$case)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(y$ncd1, y$case)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(y$ncd5, y$case)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(y$ncd2, y$case)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(y$ncd17, y$case)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(y$ncd13, y$case)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(y$e251, y$Outcome)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(y$e211, y$Outcome)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(y$e41, y$Outcome)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(y$e241, y$Outcome)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(y$e311, y$Outcome)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(y$e321, y$Outcome)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

##Recent DM
com <- table(y$dur.gp, y$case)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com


##Mean Age
describe.by(y$age, y$case)
describe.by(y$age, y$Outcome)

##Adjustments for Age and Sex
ctab(nd$age.b, nd$sex, nd$acute, type = "r", percentages=TRUE)

wm.90 <- nd$wm90
aggregate.table(wm.90, nd$age.b, nd$sex, FUN=mean, na.rm=TRUE)

com <- table(dm$age.b, dm$sex)
com <- prop.table(com)*100
com
##Medicines
op.hos.pharm <- y[y$m1g == "hospital pharmacy" & y$hosp.op > 0,]

##Mean Number of DM Meds per Person
many.drug.dm <- y$many.drug.dm
many.drug.dm1 <- dm$many.drug.dm
many.drug.dm2 <- nd$many.drug.dm
#Urban/Rural
aggregate.table(many.drug.dm, y$case, y$location, FUN=mean, na.rm=TRUE)
#Age Bracket
aggregate.table(many.drug.dm, y$case, y$age.b, FUN=mean, na.rm=TRUE)
#Duration
aggregate.table(many.drug.dm, y$case, y$dur.gp, FUN=mean, na.rm=TRUE)
#Total
aggregate.table(many.drug.dm, y$case, y$sc, FUN=mean, na.rm=TRUE)

summary(lm(ch.op1 ~ age + sex + location + dm.yrs, data = dm))
summary(lm(ch.op2 ~ age + sex + location, data = nd))
summary(lm(ch.op ~ case + sex + age + location, data=y))

##Mean Number of NonDM Meds per Person
many.drug.nondm <- y$many.drug.nondm
many.drug.nondm1 <- dm$many.drug.nondm
many.drug.nondm2 <- nd$many.drug.nondm
#Urban/Rural
aggregate.table(many.drug.nondm, y$case, y$location, FUN=mean, na.rm=TRUE)
#Age Bracket
aggregate.table(many.drug.nondm, y$case, y$age.b, FUN=mean, na.rm=TRUE)
#Duration
aggregate.table(many.drug.nondm, y$case, y$dur.gp, FUN=mean, na.rm=TRUE)
#Total
aggregate.table(many.drug.nondm, y$case, y$sc, FUN=mean, na.rm=TRUE)

summary(lm(ch.op1 ~ age + sex + location + dm.yrs, data = dm))
summary(lm(ch.op2 ~ age + sex + location, data = nd))
summary(lm(ch.op ~ case + sex + age + location, data=y))


##Mean Incremental Annual Meds Cost
y$mean.med.tot.pr[is.na(y$mean.med.tot.pr)] <- 0
mean.med <- y$mean.med.tot.pu
#Urban/Rural
aggregate.table(mean.med, y$case, y$location, FUN=mean, na.rm=TRUE)
#Age Bracket
aggregate.table(mean.med, y$case, y$age.b, FUN=mean, na.rm=TRUE)
#Duration
aggregate.table(mean.med, y$case, y$dur.gp, FUN=mean, na.rm=TRUE)

##Mean Incremental Annual strips Cost
dm$oop.strips.yr[is.na(dm$oop.strips.yr)] <- 0
mean.strips <- dm$oop.strips.yr
#Urban/Rural
aggregate.table(mean.strips, dm$case, dm$location, FUN=mean, na.rm=TRUE)
#Age Bracket
aggregate.table(mean.strips, dm$case, dm$age.b, FUN=mean, na.rm=TRUE)
#Duration
aggregate.table(mean.strips, dm$case, dm$dur.gp, FUN=mean, na.rm=TRUE)

##Proportion of people taking at least 1 drug
y$total.drug1 <- y$total.drug
y$total.drug1[y$total.drug1 > 0] <- 1
dm$total.drug1 <- dm$total.drug
dm$total.drug1[dm$total.drug1 > 0] <- 1
nd$total.drug1 <- nd$total.drug
nd$total.drug1[nd$total.drug1 > 0] <- 1

total.drug <- y$total.drug1
total.drug1 <- dm$total.drug1
total.drug2 <- nd$total.drug1

com <- table(dm$total.drug1, dm$location)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(nd$total.drug1, nd$location)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$total.drug1, dm$age.b)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(nd$total.drug1, nd$age.b)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$total.drug1, dm$dur.gp)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(dm$total.drug1, dm$sc)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com

com <- table(nd$total.drug1, nd$sc)
p.com <- prop.table(com, 2)*100
p.com <- p.com[2, ]
p.com


##Mean annual incremental exp for meds if taken
y$mean.med.tot.pr[y$total.drug == 0] <- NA
dm$mean.med.tot.pr[dm$total.drug == 0] <- NA
nd$mean.med.tot.pr[nd$total.drug == 0] <- NA
y$mean.med <- y$mean.med.tot.pr
dm$mean.med <- dm$mean.med.tot.pr
nd$mean.med <- nd$mean.med.tot.pr

mean.medy <- y$mean.med.tot.pr
mean.med1 <- dm$mean.med.tot.pr
mean.med2 <- nd$mean.med
#Urban/Rural
aggregate.table(mean.med, y$case, y$location, FUN=mean, na.rm=TRUE)
#Age Bracket
aggregate.table(mean.med, y$case, y$age.b, FUN=mean, na.rm=TRUE)
#Duration
aggregate.table(mean.med, y$case, y$dur.gp, FUN=mean, na.rm=TRUE)
#Total
aggregate.table(mean.med, y$case, y$sc, FUN=mean, na.rm=TRUE)

summary(lm(mean.med1 ~ age + sex + location + dm.yrs, data = dm))
summary(lm(mean.med2 ~ age + sex + location, data = nd))
summary(lm(mean.medy~ case + sex + age + location, data=y))


## Mean number western meds per person
y$total.drug[is.na(y$total.drug)] <- 0
dm$total.drug[is.na(dm$total.drug)] <- 0
nd$total.drug[is.na(nd$total.drug)] <- 0
total.drug <- y$total.drug
total.drug1 <- dm$total.drug
total.drug2 <- nd$total.drug
#Urban/Rural
aggregate.table(total.drug, y$case, y$location, FUN=mean, na.rm=TRUE)
#Age Bracket
aggregate.table(total.drug, y$case, y$age.b, FUN=mean, na.rm=TRUE)
#Duration
aggregate.table(total.drug, y$case, y$dur.gp, FUN=mean, na.rm=TRUE)
#Total
aggregate.table(total.drug, y$case, y$sc, FUN=mean, na.rm=TRUE)

tot.drug1 <- summary(hurdle(total.drug1 ~ age + sex + location + dm.yrs | age + sex + location + dm.yrs , data = dm))
tot.drug2 <- summary(hurdle(total.drug2 ~ age + sex + location | age + sex + location, data = dm))
tot.drug <- summary(hurdle(total.drug ~ case +age + sex + location + dm.yrs | case + age + sex + location + dm.yrs , data = dm))


##MEOP Vistis
y$meop1a <- y$meop1
y$meop1a[y$hosp.op == 0] <- NA
com <- table(y$meop1a, y$sc)
p.com <- prop.table(com, 2)*100
p.com

y$meop2a <- y$meop2
y$meop2a[y$hosp.op == 0] <- NA
com <- table(y$meop2a, y$sc)
p.com <- prop.table(com, 2)*100
p.com

y$meop3a <- y$meop3
y$meop3a[y$hosp.op == 0] <- NA
com <- table(y$meop3a, y$sc)
p.com <- prop.table(com, 2)*100
p.com

y$meop4a <- y$meop4
y$meop4a[y$hosp.op == 0] <- NA
com <- table(y$meop4a, y$sc)
p.com <- prop.table(com, 2)*100
p.com