从0开始学习R语言--Day43--Wald检验

我们构建模型来研究问题时，我们一般习惯于在研究的末尾计算误差或者其他一些指标来判断模型的效果。但其实这样做我们忽略了一个隐形前提，即模型的构建或研究的公式是有意义的。

比如我们构建一个回归模型来研究学习时间是怎么影响成绩，自变量前的系数代表影响程度，而Wald检验，本质是就是在判断这个系数是不是真的不等于0，也就是学习时间是不是真的影响成绩，这能帮我们提前排查很多问题，如公式编排错误，变量没有标准化等一些准备工作或细节工作。

以下是用于一些模型中的Wald检验的例子：

# 使用car包进行Wald检验
library(car)
# 使用aod包进行Wald检验

library(aod)

# 生成模拟数据集
set.seed(123)
n <- 100
x1 <- rnorm(n)
x2 <- rnorm(n)
y <- 2 + 1.5*x1 + 0.5*x2 + rnorm(n, sd = 0.5)
data <- data.frame(y, x1, x2)

# 拟合线性模型
model <- lm(y ~ x1 + x2, data = data)



# 检验x1的系数是否为0（单个参数检验）
linearHypothesis(model, "x1 = 0")

# 检验x1和x2的系数是否都为0（多个参数检验）
linearHypothesis(model, c("x1 = 0", "x2 = 0"))

# 查看模型摘要（包含Wald检验结果）
summary(model)

# 生成二分类数据
set.seed(123)
n <- 200
x1 <- rnorm(n)
x2 <- rnorm(n)
z <- 0.5 + 1.2*x1 - 0.8*x2
prob <- 1/(1+exp(-z))
y <- rbinom(n, 1, prob)
data_binary <- data.frame(y = as.factor(y), x1, x2)

# 拟合逻辑回归模型
glm_model <- glm(y ~ x1 + x2, family = binomial(), data = data_binary)

# 使用summary函数进行Wald检验
summary(glm_model)


wald.test(b = coef(glm_model), Sigma = vcov(glm_model), Terms = 2:3)


# 生成计数数据
set.seed(123)
n <- 150
x <- rnorm(n)
mu <- exp(1 + 0.5*x)
y <- rpois(n, mu)
data_count <- data.frame(y, x)

# 拟合泊松回归模型
poisson_model <- glm(y ~ x, family = poisson, data = data_count)

# 进行Wald检验
summary(poisson_model)

# 检验x的系数是否为0.5

wald.test(b = coef(poisson_model), Sigma = vcov(poisson_model), Terms = 2, H0 = 0.5)

输出：

Call:
lm(formula = y ~ x1 + x2, data = data)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.93651 -0.33037 -0.06222  0.31068  1.03991 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  2.06753    0.04807   43.01   <2e-16 ***
x1           1.43341    0.05243   27.34   <2e-16 ***
x2           0.51191    0.04950   10.34   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.4756 on 97 degrees of freedom
Multiple R-squared:  0.8952,	Adjusted R-squared:  0.893 
F-statistic: 414.2 on 2 and 97 DF,  p-value: < 2.2e-16

Call:
glm(formula = y ~ x1 + x2, family = binomial(), data = data_binary)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)   0.3473     0.1720   2.019   0.0435 *  
x1            1.1766     0.2188   5.377 7.57e-08 ***
x2           -1.0146     0.2010  -5.049 4.45e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 274.83  on 199  degrees of freedom
Residual deviance: 206.58  on 197  degrees of freedom
AIC: 212.58

Number of Fisher Scoring iterations: 4

Call:
glm(formula = y ~ x, family = poisson, data = data_count)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  0.98722    0.05245   18.82   <2e-16 ***
x            0.55921    0.04883   11.45   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 310.58  on 149  degrees of freedom
Residual deviance: 178.93  on 148  degrees of freedom
AIC: 568.19

Number of Fisher Scoring iterations: 5

输出中的p值都远小于0.001，但重要的是，这是个很好且工作量并不大的检查工作，而且这也可以体现出我们模型结果是否显著。