阻止空数据在R中输入lm()调用



我正试图想出一种机制来防止lm()输出中出现空结果。确切地说,我想先找到它们,然后防止它们被输入到新的lm()调用中。

例如,在下面的示例中,cf.type99:time3&tmp的输出中的cf.type99:time4NAs。因此,在检测到它们之后,我希望防止它们被包括在第二个lm()调用(new_tmp(中。

这在R中可能吗?

p.s.我希望new_tmp不适合cf.type99:time3&cf.type99:time4在其model.matrix()中,而不仅仅是从其输出中物理删除NA。

d5 <- read.csv('https://raw.githubusercontent.com/rnorouzian/m/master/v14.csv')
d5[c("cf.type","time")] <- lapply(d5[c("cf.type","time")], as.factor)
tmp <- lm(dint~cf.type*time, data = d5)
(coef.na <- is.na(coef(tmp))) # detects the `NA` in the output:
# cf.type2:time3  cf.type3:time3  cf.type8:time3 cf.type99:time3  cf.type1:time4 
#          FALSE           FALSE           FALSE            TRUE           FALSE 
# cf.type2:time4  cf.type3:time4  cf.type8:time4 cf.type99:time4 
#          FALSE           FALSE           FALSE            TRUE 
new_tmp <- lm(dint~cf.type*time, data = d5) # NOW PREVENT `cf.type99:time3` & `cf.type99:time4` from entering new_tmp

您可以将列从model.matrix、中删除

X.star <- model.matrix(tmp)[,!is.na(tmp$coefficients)]

然后使用lm.fit计算系数,

y <- d5$dint
fit.star <- lm.fit(X.star, y)
(beta.hat <- fit.star$coe)
# (Intercept)        cf.type1        cf.type2        cf.type3        cf.type8 
# 0.72824752     -0.20633963      0.72317588     -0.08369434      0.17387840 
# cf.type99           time2           time3           time4  cf.type1:time2 
# -0.53932949     -0.49471558     -0.34560481      0.01830757      0.28441094 
# cf.type2:time2  cf.type3:time2  cf.type8:time2 cf.type99:time2  cf.type1:time3 
# -0.44254640      0.77070527      0.71423737      0.51203670      1.15793712 
# cf.type2:time3  cf.type3:time3  cf.type8:time3  cf.type1:time4  cf.type2:time4 
# -0.66306413      0.18818360     -0.12287313     -0.06825029     -1.06941453 
# cf.type3:time4  cf.type8:time4 
# -0.37337606     -0.06063658 

以及可选的手动标准误差。

## Standard errors
y.hat <- X.star %*% beta.hat
e.hat <- y - y.hat
n <- dim(X.star)[2]
m <- dim(X.star)[1]
sig2 <- (t(e.hat) %*% e.hat) / (m - (n + 1))
res.se <- sqrt(sig2)
diag.xx <- diag(solve(t(X.star) %*% X.star))
beta.se <- as.vector(res.se)*sqrt(diag.xx)

结果

cbind(Estimate=beta.hat, `Std. Error`=beta.se, `t value`=beta.hat/beta.se, 
`Pr(>|t|)`=2*pt(-abs(beta.hat/beta.se), df=fit.star$df.residual))
#                    Estimate Std. Error     t value     Pr(>|t|)
# (Intercept)      0.72824752  0.1420687  5.12602302 6.252427e-07
# cf.type1        -0.20633963  0.2460702 -0.83853955 4.025946e-01
# cf.type2         0.72317588  0.4921405  1.46945012 1.430713e-01
# cf.type3        -0.08369434  0.2751149 -0.30421597 7.612372e-01
# cf.type8         0.17387840  0.2231278  0.77927713 4.366141e-01
# cf.type99       -0.53932949  0.6935709 -0.77761265 4.375931e-01
# time2           -0.49471558  0.1813017 -2.72868700 6.847850e-03
# time3           -0.34560481  0.1834099 -1.88432993 6.077612e-02
# time4            0.01830757  0.2391373  0.07655675 9.390424e-01
# cf.type1:time2   0.28441094  0.3152098  0.90229103 3.678421e-01
# cf.type2:time2  -0.44254640  0.5208981 -0.84958339 3.964364e-01
# cf.type3:time2   0.77070527  0.4209305  1.83095614 6.839524e-02
# cf.type8:time2   0.71423737  0.4165119  1.71480662 8.772179e-02
# cf.type99:time2  0.51203670  0.8460178  0.60523161 5.456192e-01
# cf.type1:time3   1.15793712  0.4035018  2.86971958 4.489380e-03
# cf.type2:time3  -0.66306413  0.5522639 -1.20062913 2.311248e-01
# cf.type3:time3   0.18818360  0.4218428  0.44609884 6.559437e-01
# cf.type8:time3  -0.12287313  0.3559152 -0.34523149 7.302345e-01
# cf.type1:time4  -0.06825029  0.3912267 -0.17445200 8.616630e-01
# cf.type2:time4  -1.06941453  0.6066406 -1.76284686 7.924862e-02
# cf.type3:time4  -0.37337606  0.4196728 -0.88968377 3.745613e-01
# cf.type8:time4  -0.06063658  0.3477727 -0.17435692 8.617377e-01

相关内容

  • 没有找到相关文章

最新更新