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

Question

不确定这个问题在交叉验证中是否会更好，但我认为它既是一个纯粹的统计问题，也是一个编程问题。

我有一个 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。

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

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，但回归似乎并不关心。

我有一种感觉，这可以通过对我的脚本进行小幅调整来完成。我怎么做？

Answer 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 的年份没有条形图。为零，但这只是因为绘图准确地反射(reflect)了值——当条形的高度为零时，您将看不到条形。