强制 R 在计数与年份的回归中包含 0 作为值答案

【问题标题】：Force R to include 0 as a value in a regression of counts vs year强制 R 在计数与年份的回归中包含 0 作为值
【发布时间】：2017-11-02 17:54:07
【问题描述】：

不确定这个问题在 Cross Validated 中是否会更好，但我认为它既是一个编程问题，也是一个纯统计问题。

我有一个 102 x 1147 的数据框，其中有年份（1960 年到 2016 年之间），每条记录都是一篇科学论文。我计算每年在某些主题内发表的论文数量（以特定列中的值为指导），我想计算一年的线性斜率和论文数量的年度计数。

这是我的脚本，首先是线性模型，然后是情节：

# THEME 1 (POPABU)
sub2=subset(as.data.frame(table(sysrev60[,c("YR","POPABU")])),
        POPABU==1,select=c(1,3))
sub2$YR<-as.numeric(paste(sub2$YR))

lm_eqn <- function(df){
  m <- lm(Freq ~ YR, sub2);
  eq <- substitute(italic(y) == a + b %.% italic(x)*","~~italic(r)^2~"="~r2, 
               list(a = format(coef(m)[1], digits = 2), 
                    b = format(coef(m)[2], digits = 2),
                    r2 = format(summary(m)$r.squared, digits = 3)))
  as.character(as.expression(eq));                 
}

ggplot(sub2, aes(x=YR,y=Freq)) + 
  scale_y_continuous(limit=c(0,20),expand=c(0, 0)) +
  scale_x_continuous(breaks=c(1960,1965,1970,1975,1980,1985,1990,1995,2000,
                          2005,2010,2015),labels=c(1960,1965,1970,1975,1980,1985,
                                                   1990,1995,2000,2005,2010,2015)) +
  geom_bar(stat='identity') + 
  geom_text(x = 1960, y = 16, label = lm_eqn(df), size=5,hjust=0, parse = TRUE) +
  stat_smooth(method="lm",col="red") +
  xlab(" ") + ylab("No of papers") +
  annotate("text",x=1960,y=18,label="THEME 1",
       family="serif",size=7,hjust=0,color="darkred")

我的问题是这个程序只计算年份和计数之间的线性关系 > 0。有很多年的论文数等于 0，我需要回归来覆盖同一时期（1960- 2016）对于我正在研究的所有 25 个不同主题，即我需要强制回归包含 0，因为每年论文数为 0。

我已经制作了与我想研究其发表率的每个主题相对应的大型数据框的子集。这是我的“sub2”数据框的DPUT：

dput(sub2)
structure(list(YR = c(1960, 1961, 1962, 1963, 1964, 1965, 1966, 
1967, 1968, 1969, 1970, 1971, 1972, 1973, 1974, 1975, 1976, 1977, 
1978, 1979, 1980, 1981, 1982, 1983, 1984, 1985, 1986, 1987, 1988, 
1989, 1990, 1991, 1992, 1993, 1994, 1995, 1996, 1997, 1998, 1999, 
2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 
2011, 2012, 2013, 2014, 2015, 2016), Freq = c(0L, 0L, 0L, 0L, 
0L, 1L, 0L, 0L, 0L, 0L, 1L, 1L, 1L, 0L, 0L, 0L, 2L, 1L, 0L, 1L, 
3L, 0L, 1L, 0L, 2L, 0L, 3L, 0L, 1L, 0L, 1L, 0L, 0L, 1L, 1L, 2L, 
0L, 2L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 2L, 0L, 1L, 
1L, 1L, 2L, 3L, 5L)), .Names = c("YR", "Freq"), row.names = 58:114, class = "data.frame")

如您所见，我的数据框中似乎有明确的 0，但回归似乎并不在意。

我感觉这可以通过对我的脚本稍作调整来完成。我该怎么做？

【问题讨论】：

听起来你需要用显式零值填充隐式零值。看tidyr::complete。
同意；除非您拥有庞大的稀疏数据集（其中“巨大”= 数百或数千列、数百万行），否则完成数据要比诱使回归隐含地包含这些点要容易 1000 倍
这可能只是一个错字吗？ lm_eqn 没有使用它的参数 df。
你为什么用lm建模？您应该使用 GLM，例如 Poisson 回归。

标签： r linear-regression

【解决方案1】：

到目前为止，您所做的考虑了零，我们可以通过手动计算系数来仔细检查，以防您认为 lm() 出于某种原因做了一些奇怪的事情：

# Make sure zeros are there:
sub2$Freq
[1] 0 0 0 0 0 1 0 0 0 0 1 1 1 0 0 0 2 1 0 1 3 0 1 0 2 0 3 0 1 0 1 0 0 1 1 2 0 2
[39] 0 0 0 1 0 0 0 0 0 1 0 2 0 1 1 1 2 3 5
# Yep
X <- cbind(rep(1, nrow(sub2)), sub2$YR) # add a column of 1s for intercept
solve(t(X) %*% X) %*% t(X) %*% sub2$Freq # (X'X)^-1 X'Y -- OLS formula

            [,1]
[1,] -38.1778584
[2,]   0.0195748

考虑到四舍五入，这与您发布的代码产生的图上显示的内容相匹配：

当我们使用包括零在内的所有值时，截距约为-38，年份系数约为0.02。所以，那里绝对没有错。可能导致您认为它忽略零的原因是Freq 为零的年份没有条形图，但这只是因为该图准确地反映了这些值——当条形图的高度为零时，您将无法看到栏。

【讨论】：