心理健康正日益成为一个传统上被掩盖的话题。 我们已经开始理解心理健康对生产力，总体健康，人际关系和身体健康的影响，并将更多的注意力放在心理健康上。 甚至雇主也开始更加重视提供工作环境和条件，以保持员工尽可能快乐和健康，而先进的公司则提供津贴，例如每月按摩，餐饮，免费和补贴的健身房以及无限制的假期。 鉴于此，作为一个对发展我的编程和数据科学技能感兴趣的人，我正在研究的Coursera Statistics专业的数据可视化和分析项目也部分激发了我的动力。 我决定使用疾病控制与预防中心（ CDC ）的数据，提取2013年行为风险因素监视系统数据 （BRFSS）。 该数据是通过每月一次电话采访收集的，是最大的数据集之一，致力于了解预防健康实践和风险行为的趋势，这些趋势与慢性病伤害和影响美国成年人的可预防传染病有关。 我分析的目的是了解收入水平和所报告的心理健康之间的关系，以及每年收入多少，身体活动程度以及所吃水果和蔬菜的量对心理健康的影响。
提取数据并了解结构
install.packages("ggplot2")
install.packages("mice")
install.packages("effects")
library(ggplot2)
library(MASS)
library(mice)
library(data.table)
library(dplyr)
library(effects)
load("~/Documents/R Stats Assignments/brfss2013.RData")
上传必要的库并加载数据后，我想查看数据集的一般形状和结构，以了解维数，每列中的数据类型以及为获取数据而需要进行的清理以一致的格式进行分析。
> dim(brfss2013)
[1] 491775 330
#check structure of columns
str(brfss2013$columname)
#check for NA values
sum(is.na(brfss2013$columname))
在查看了提供的代码本/索引之后 ，我对需要与之交互的列有了很好的了解，并决定专门关注这些列。 其中包括具有以下数据的列：每周进行身体活动的分钟，每周食用的水果，每周食用的蔬菜，BMI，年收入范围以及当然的心理健康状况，这些数据由一个人据称身体状况不佳的天数来衡量在一定时期内的心理健康。 由于NA值会混淆分析的准确性并可能导致无响应偏差，因此我在列中寻找NA值，这些值介于数据的4.8％至34％之间。 我认为这些是不可忽略的缺失值，不能简单地从我的分析中删除。 对于在缺少34％数据的情况下的数值，我假设仅填写标准均值或中位数会歪曲我的分析并降低其准确性。 我选择了一种易于实现且相对准确的插补方法来代替我的NA值。
贝叶斯线性回归的多重插补
准确估算值的最快方法之一是通过线性回归算法，该算法尝试根据预定义的预测变量列表预测潜在值。 使用频繁方法（常规线性回归），我们假设自变量之间存在线性关系，另外还假设存在足够的值来进行准确的预测。 但是，对于这种估算，我选择了贝叶斯模型，通过这种方法，我们假设响应和预测变量都是随机的（将估算的不确定性考虑在内），尝试查找特定事件的概率，并根据先验概率更新概率。数据并根据这些概率进行预测。 合并先前信息以更新概率的能力使贝叶斯模型在做出预测时更灵活，并且可能更准确，尤其是对于较小的样本。
多重插补让我遍历插补模型n次，根据所使用的预测模型提供不同的可能插补。 通常，这是一个三阶段过程，涉及实际估算，用可能的预测替换每个NA，然后分析并汇总所有结果。
brfss2 <- brfss2013[,c('X_bmi5','X_drnkmo4','X_frutsum','X_vegesum','fc60_','pa1min_')]
imp_model <- mice(brfss2,method="norm",m=5)
#fill na function
fillna_fun <- function(data,columns){
df <- setNames(data.frame(rowMeans(squeeze(imp_model$imp[[columns]], bounds = c(0,max(brfss2013[[columns]]))))),"col2")
brf <- setNames(data.frame(data[[columns]]),"col2")
brf$col1 <- rownames(brf)
df$col1 <- rownames(df)
setDT(brf)[df,col2 :=i.col2,on=.(col1)]
brf$col2
}
brfss2013$pa1min_ <- fillna_fun(brfss2,"pa1min_")
brfss2013$X_bmi5 <- fillna_fun(brfss2,"X_bmi5")
brfss2013$X_drnkmo4 <- fillna_fun(brfss2,"X_drnkmo4")
brfss2013$X_frutsum <- fillna_fun(brfss2,"X_frutsum")
brfss2013$X_vegesum <- fillna_fun(brfss2,"X_vegesum")
通过链式方程式（小鼠）进行多元插补可使用贝叶斯线性回归轻松进行插补。 我选择缩小适用于具有分析所需数值的列的插补，并将插补运行默认次数-5。
我注意到贝叶斯线性回归返回了一些估算值的负数，考虑到这些列所测量的数据，这是不可能的。 例如，一个人不可能每周花费-300分钟进行锻炼。 为了解决这个问题，我使用了挤压功能来添加约束，以将预测限制在合理的自定义范围内。
与小鼠文档相反，我选择获取所有5个估算值的平均值，并用其填充NA值。 我完全理解这不是应该使用机器学习的方式，我错过了合并合并归因的基本步骤，并且有人认为我最好选择一个归因于平均值的归因。 尽管如此，我还是采用了该选项以简化应用程序，并希望收到读者对该选项的进一步反馈。
我通过使用setDT函数将函数中的向量强制到表中并按索引将向量连接起来来完成函数 ，因为均值插补向量仅具有对应于我原始数据帧中NA值的索引值。
映射
我分析的重要步骤涉及将某些列的数值映射到组和范围中。 例如，我没有看到BMI编号，而是想查看特定BMI属于哪个组/范围（例如25 =正常）。
brfss2013$income2 <- as.character(brfss2013$income2)
brfss2013$X_bmi5 <- brfss2013$X_bmi5/100
brfss2013$healtheat <- (brfss2013$X_frutsum+brfss2013$X_vegesum)/100
labels <- c('Excellent','Good','Ok','Bad','Very Bad')
breaks <- c(0,5,10,15,25,10000)
bmiLabs <- c('10','20','30','40','50','60','>60')
bmiBreaks <-c(0,10,20,30,40,50,60,10000)
activLabs <-c('0-200','200-500','500-1000','1000-2000','2000-4000','4000-10000','>10000')
activBreaks <-c(0,200,500,1000,2000,4000,10000,100000)
brfss2013 <- brfss2013 %>%
mutate(mentalHealth = cut(menthlth,breaks=breaks,labels=labels,include.lowest=TRUE)) %>%
mutate(bmiLev = cut(X_bmi5, breaks=bmiBreaks,labels=bmiLabs,include.lowest = TRUE)) %>%
mutate(physLev = cut(pa1min_, breaks=activBreaks,labels=activLabs,include.lowest = TRUE)) %>%
mutate(incomeLev = case_when(grepl("15|20",income2)~"0-$20k",
grepl("25|35",income2)~"25-$35k",
grepl("50",income2)~"35-$50k",
income2 %in% "Less than $75,000" ~ "50-$75k",
grepl("more",income2)~">$75k"
))
由于cut旨在采用数值矢量并将其根据自定义断点集将其拆分为bin，因此我决定使用cut函数 ，定义中断和标签以根据标签定义映射数据。 对于收入水平列，我将值转换为字符以使我能够使用grepl 。 这是一种模式匹配功能，我通过查找与用于标识每一行收入范围的自定义单词匹配的关键字来创建收入范围。
需要注意的重要一点是，我对心理健康作了一些假设。 我假设某个人在给定时期内仅报告过0至5次心理健康问题，则被归类为处于良好的心理健康状态，具有5至10的健康状况，10至15的正常，15至25的不良状况以及任何超出那非常糟糕。 这是基于观察数据并找到理想的断裂点进行分类的基础。
心理健康与年收入的关系
我认为叠加条形图是可视化心理健康与年收入范围之间关系的最有效和最有吸引力的方法。
#replace NA with missing
brfss2013$mentalHealth <- forcats::fct_explicit_na(brfss2013$mentalHealth, na_level = "Missing")
#convert income back to factor
brfss2013$incomeLev <- as.factor(brfss2013$incomeLev)
brfss2013 <- subset(brfss2013, !is.na(incomeLev))
brfss2013 %>%
add_count(incomeLev) %>%
rename(count_inc = n) %>%
count(incomeLev, mentalHealth, count_inc) %>%
rename(count_mentalHealth = n) %>%
mutate(percent= count_mentalHealth / count_inc) %>%
mutate(incomeLev = factor(incomeLev,
levels=c('0-$20k','25-$35k','35-$50k','50-$75k','>$75k')))%>%
ggplot(aes(x= incomeLev,
y= count_mentalHealth,
group= mentalHealth)) +
xlab('Annual Income')+ylab('Number of People')+
geom_bar(aes(fill=mentalHealth),
stat="identity",na.rm=TRUE)+
# Using the scales package does the percent formatting for me
geom_text(aes(label = scales::percent(percent)),position = position_stack(vjust = 0.5))+
theme_minimal()
使用dplyr软件包 ，我使用％>％运算符编写了多个运算，这些运算将按组统计收入水平和心理健康，按收入水平查找报告的心理健康数字的相对百分比，根据我的自定义对条形图进行分组，并可视化堆叠的条形图，并在每个图上添加百分比标签。
使用我对构成良好心理健康的假设，数据似乎可以证实人们的预期，您赚的钱越多，您处于更好的心理状态的可能性就越大。
收入水平的NA值该怎么办？
从年收入范围可视化显示的心理健康分布中，很明显在“收入水平”列中有很多NA值。 尽管我可以简单地降低NA值，但假设它们不会改变我的分析准确性并继续前进，但我相信很多人可能不愿意谈论其收入水平。 这意味着很大一部分资产净值来自不愿意通过电话分享其年收入的人。 可能的情况是，那些不愿意分享其收入信息的人处于较低的收入范围（0至20k）或收入很高的人群（超过70k），而忽略了这些行，这些行占收入数据的15％以上，容易引入无响应偏差。 我想通过机器学习来估算这些NA的潜在价值。 这次我选择通过多元回归模型来估算我的有序数据。 这种方法使用比例赔率逻辑回归模型，其机制将在以下段落中详细讨论。
#filling na values of income level column
brfss <- brfss2013[,c('incomeLev','healtheat','X_age_g','employ1','renthom1','sex','physLev')]
ordered_brfss <-mice(brfss, m=1, method='polr', maxit=1)
fillna_inc <- function(data,columns){
df <- setNames(data.frame(ordered_brfss$imp[[columns]]),"col2")
brf <- setNames(data.frame(data[[columns]]),"col2")
brf$col1 <- rownames(brf)
df$col1 <- rownames(df)
setDT(brf)[df,col2 :=i.col2,on=.(col1)]
brf$col2
}
brfss2013$incomeLev_ <- fillna_inc(brfss2013,"incomeLev")
为了进行这种估算，我重点介绍了一些我认为可以最大程度预测潜在收入的专栏。 我特别关注了个人报告是否租房或拥有房屋，是否雇用他/她，性别，年龄，体育活动水平以及受访者报告的水果和蔬菜数量。 我选择运行单个插补，纯粹是为了易于使用。
按收入水平可视化心理分布
我想看看按收入水平划分的心理健康分布如何，并建立了另一个可视化视图。
比例赔率
将心理健康划分为好，好，优秀等组的行为，将本专栏转化为一个具有5级水平的因子。 这些因素根据级别排序。 考虑到此列的结构，我选择使用比例赔率逻辑回归。 为了进行分析，我需要选择一个模型，该模型可以帮助我理解在给定的独立变量的情况下某个人的心理健康处于特定类别的可能性，并了解这些变量对概率的影响。 尽管线性回归在某种程度上可以解决二进制分类问题（例如，如果响应变量是电子邮件还是垃圾邮件），但是当您为响应变量有多个类别时，线性回归就会崩溃。 序数逻辑回归使用logit函数将线性模型转换为满足有序响应类别，从而确保返回的概率在0到1的范围内。 我本来可以选择更灵活的多项式逻辑回归模型 ，但这是基于这样的假设，即响应变量中的类别不能以任何有意义的方式进行排序-我认为该假设不适用于我的心理健康响应变量。
我的模型依赖于计算比例优势比。 简而言之，比例优势比是一种用于有序逻辑回归的工具，可通过预测在特定值类别下给定响应变量属于特定类别的条件概率来帮助对自变量与有序响应因子/类别之间的关系做出假设。观察到的独立变量。
在这种情况下，J表示我们要预测的类别/因素。 上述公式吸收了n个因子，其中n取决于我在因变量中拥有的因子数。 数字指的是我用来构建模型的自变量。
我使用此公式来了解与我们的结果（心理健康）关系最大的自变量，并试图了解自变量对心理健康的影响程度。
brfss2_model = polr(mentalHealth ~ incomeLev+bmiLev+X_drnkmo4+healtheat+physLev,data=brfss2013,Hess=TRUE) #Hessian used to get standard errors
我使用polr（）函数将比例优势物流回归拟合到响应变量（心理健康）和预测变量。 在这种情况下，我决定将BMI，一个月内消耗的水果和蔬菜总量，体力活动水平，年收入范围和一个月内饮用的酒精饮料包括在内。 我通过将TRUE参数添加到Hessian中来包括观察到的信息矩阵，以便获得标准误差并尝试评估该模型对数据的适用性。
(ctable<-coef(summary(brfss2_model)))
#calculating p_value by comparing the t-value against the stnd norm distr similar to a z-test
p<-pnorm(abs(ctable[, "t value"]),lower.tail = FALSE)*2
#combining p-value
(ctable<-cbind(ctable,"p value"=p))
我着手找到回归系数并计算p_value以确定这些结果的重要性。 回归系数的范围从每月含酒精饮料的0.001到一组BMI水平的2.55略高，p值表明具有统计学意义。 根据我的理解，回归系数越接近零，就越表明预测变量的分布对于我们的响应变量的每个级别完全相同。 对于接近零的系数，这宽松地表示在响应变量上向自变量添加更多单位的效果接近于零。
可视化自变量的影响
为了清楚显示自变量对心理健康的影响，我决定使用“ 影响”软件包 。 由于polr函数返回估计的回归系数，因此我使用此程序包来帮助我解释polr结果 。
plot(Effect(focal.predictors = c("incomeLev","bmiLev"),brfss2_model))
在尝试了几种组合之后，我选择显示年收入和BMI效果图。 我将这两个变量用作焦点预测变量，使其他预测变量具有典型值（固定和平均）。
根据结果，无论您的收入水平如何，随着BMI的增加，处于良好心理健康状态的可能性都会降低。 有趣的是，当您接近年收入$ 75,000时，下降的趋势逐渐减弱。 我们可以看到，在非常糟糕的心理健康状态下，类似的影响较小，您越富有，BMI升高对您的心理健康的影响就越小。 如果您要增加体重，那么在增加收入的同时，您的精神状态可能会更好。
除了我所知甚少之外，我正在学习的最重要的课程之一是数据分析和数据科学是一个迭代过程。 您可以尝试各种变量和模型，根据不断变化的假设清理方法和可视化效果，并测试这对结果及其准确性有何影响。
随时联系Twitter或在Twitter @Emmoemm上发送任何反馈
1094
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
