R语言一般线性模型（涉及因变量是虚拟变量（哑变量））

Posted 2022-11-20 影月风格厂

R语言的一般线性模型

R语言的一般线性模型用函数：lm()，即可轻松实现。

例子

建立一般线性模型

NC.glm1 = lm(fmri.SFG_R_CerebellumGM_L ~ age + gender + race,data = NC)
summary(NC.glm1)

这里解释一下变量（我直接copy我项目里面的两行代码），因变量y就是fmri.SFG_R_CerebellumGM_L，自变量x有三个age，gender（因子变量：female、male）和race（因子变量：AA、Asian、Cauc、Hispanic、Other），这些都是数据集NC的某一列。建立y和x的一般线性模型，得到的结果如下所示：

Call:
lm(formula = fmri.SFG_R_CerebellumGM_L ~ age + gender + race, 
    data = NC)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.54058 -0.20018 -0.01524  0.15821  1.00124 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)  
(Intercept)   0.442044   0.168781   2.619   0.0103 *
age           0.003510   0.006861   0.512   0.6102  
genderMale    0.023009   0.060557   0.380   0.7049  
raceAsian     0.008780   0.205792   0.043   0.9661  
raceCauc     -0.073530   0.067110  -1.096   0.2761  
raceHispanic  0.035361   0.148205   0.239   0.8120  
raceOther    -0.107146   0.168367  -0.636   0.5261  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.2834 on 90 degrees of freedom
Multiple R-squared:  0.02018,	Adjusted R-squared:  -0.04514 
F-statistic: 0.309 on 6 and 90 DF,  p-value: 0.9308

对模型结果进行分析

一般线性模型（针对普通数值变量类型）

我们知道回归分析就是建立一个一般的线性函数来拟合：y=a1x1+a2x2+a3*x3+b
从上面summary的结果中，我们要得到几个重要的系数（coefficients）：
intercept -> b;
age -> a1；
写到这里，可能部分同学会不懂了，为什么后面的系数有这么多呢？理论上应该就b，a1，a2，a3四个系数啊！

一般线性模型（针对因子变量类型）

在R中，如果自变量是因子变量类型，lm()函数会自动处理成虚拟变量，也叫哑变量（dummy variable）。所以本例中的自变量gender和race在lm()函数中自动处理成了虚拟变量。
其实虚拟变量就是将因子变量的几个不同的level处理成对应的01变量集合。比如：gender

gender	genderFemale	genderMale
female	1	0
male	0	1

从表中我们可以看见本来是因子变量的gender变量，如今映射成为两个01变量genderFemale和genderMale。
同理：race

race	raceAA	raceAsian	raceCauc	raceHispanic	raceOther
AA	1	0	0	0	0
Asian	0	1	0	0	0
Cauc	0	0	1	0	0
Hispanic	0	0	0	1	0
Other	0	1	0	0	1

从表中可见因子变量映射成为虚拟变量之后，一般因子变量中有几个level就会生成几个虚拟变量。
那么问题来了，为什么从之前的summary的结果我们没有看到genderFemale和raceAA呢?
答案是：为了避免多重共线性，对于level=n的分类变量只需选取其任意n-1个虚拟变量。因此需要选取一个reference level，gender中选取的是female，race中选的是AA，因此summary的结果就可以完全解释了。

利用模型进行回归分析

有时候我们会通过predict()函数对该模型在新的测试集中进行回归。但我们需要做一个矫正的时候，需要自己写一个线性模型进行回归分析。了解了因子变量和虚拟变量（dummy 变量）的关系之后，这并不难。代码如下：

NCgenderD = nnet::class.ind(NC$gender)    #generate dummy variable
NCgenderD = NCgenderD[,-1    ]#to avoid collinear,we drop the fisrt variate
NCraceD = nnet::class.ind(NC$race)    #generate dummy variable
NCraceD = NCraceD[,-1]     #to avoid collinear,we drop the fisrt variate
NC_D = cbind(NC$age,NCgenderD,NCraceD)       # merge independent variables

#Actually,I have lots of dependent variables to do correction(or regression)
coef_lm = sapply(37:3232,function(i)
  NC.glm = lm(NC[,i] ~ age + gender + race, data = NC)
  return(NC.glm$coefficients)
)
coef_lm1 = coef_lm[-1,]#drop the intercept ,so as to do matrix multiplication
RegressNC = sapply(37:3232,function(i)
  return(coef_lm[1,i]+(NC_D %*% coef_lm1[,i-36]))
)

总结

在R语言实现一般线性模型是很方便的，本文通过手写带虚拟变量的线性回归函数更好地理解R语言中对因变量是因子变量的处理方式，即转化为虚拟变量（哑变量dummy variable）进行处理。我们在利用模型的coefficients进行回归的时候也需要注意一下。
PS：为了理解这个虚拟变量，我找了很多资料，说起容易还是挺花时间的。差不多有大半天专门去找虚拟变量的相关知识，为了使得别人不再像我一样，所以又花了小半天时间写了这个markdown。