Visualobesity3.R

#Stepwise regression was used to predict obesity rate. Starting from saturated model
#variables not found significant were removed using stepwise regression.
data = fulldb[,-c(1:3)]
complete.data = data[complete.cases(data),]
model.saturated = lm(PCT_OBESE_ADULTS10 ~ ., data = complete.data)
model.empty = lm(PCT_OBESE_ADULTS10 ~ 1, data = complete.data)
scope = list(lower = formula(model.empty), upper = formula(model.saturated))
backwardAIC = step(model.saturated, scope, direction = "backward", k = 2)

#Used to predict obesity rate with multiple linear regression
predfunc <- function(GROC, Conv, FF, LACCESS, MEDHHIN, PCT18, FOODINS, FARMRT, 
                     VEGFARM, DIABETE, HSACT,POVRT, PCT65, RECFAC, Full) {
    newdata = data.frame(PCT_LACCESS_POP10 = LACCESS, GROCPTH12 = GROC, CONVSPTH12 = Conv,
                         FFRPTH12 = FF, FOODINSEC_10_12 = FOODINS, VEG_FARMS07 = VEGFARM,
                         FMRKTPTH13 = FARMRT,  PCT_DIABETES_ADULTS10 = DIABETE, FSRPTH12 = Full,
                         PCT_HSPA09 = HSACT, MEDHHINC10 = MEDHHIN, POVRATE10 = POVRT,
                         PCT_65OLDER10 = PCT65, PCT_18YOUNGER10 = PCT18, RECFACPTH12 = RECFAC)
    predict(backwardAIC, newdata, interval = "prediction")
}

#Used for slider bar for prediction, low bounds is mean - 3*standard dev 
rlow <- function(xinput) {
    if (mean(xinput, na.rm = T) > 3*sd(xinput, na.rm = T)) {
        round(mean(xinput, na.rm = T)-3*sd(xinput, na.rm = T),2)
    } else {
        print(0)
    }
}   
#Used for slider bar for prediction, high bounds is mean + 3*standard dev 
rhi <- function(xinput) {
    round(mean(xinput, na.rm = T)+3*sd(xinput, na.rm = T),2)
}   
#Used for slider bar for prediction, starting value is mean 
valme <- function(xinput) {
    round(mean(xinput, na.rm = T),2)
}   

coef = summary(backwardAIC)$coefficients
coef = as.data.frame(coef)
coefnames = c("Intercept", "Percent low access to store", 
                   "Convenience stores/1,000 pop, 2012", "Fast-food restaurants/1,000 pop, 2012",
                   "Household food insecurity (%)", "Full-service restaurants/1,000 pop, 2012",
                   "Farmers' markets/1,000 pop, 2013", "Vegetable farms, 2007", "Adult diabetes 
                   rate, 2010", "High schoolers physically active (%), 2009", "Rec & fitness 
                   facilities/1,000 pop, 2012", "Median household income, 2010","Poverty rate, 2010", 
                   "% Population 65 years or older, 2010", "% Population under age 18, 2010")
coef = mutate(coef, signif = c("***", "**", ".", "***","***","***", ".", "***","***","***",
                               "*", "***","***","***","***"))
coef = cbind(coefnames, coef)
coef[,2:4] = round(coef[,2:4], 5)
sigcodes = "Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1"
sumreg = "Residual standard error: 2.549 on 2499 degrees of freedom
Multiple R-squared:  0.6671,	Adjusted R-squared:  0.6652 
F-statistic: 357.6 on 14 and 2499 DF,  p-value: < 2.2e-16"

vifs = car::vif(backwardAIC)