r语言 - AOV在使用日期和因子时的功能不同

Question

我正在运行一个简单的方差分析，但当Year是一个因素时，p值没有出现。如果我将Year编码为数字，p值就会显示出来。我真的希望Year是一个因子，而不是像日期那样的连续变量。

dat <- structure(list(Year = structure(1:26, levels = c("1994", "1995", 
"1996", "1997", "1998", "1999", "2001", "2002", "2003", "2004", 
"2005", "2006", "2007", "2008", "2009", "2010", "2011", "2012", 
"2013", "2014", "2015", "2016", "2017", "2018", "2019", "2021"
), class = "factor"), no = c(1, 1, 0, 1, 2, 3, 4, 14, 28, 0, 
2, 2, 6, 1, 0, 1, 0, 0, 0, 2, 5, 0, 4, 3, 0, 0), yes = c(3L, 
5L, 2L, 1L, 8L, 16L, 30L, 19L, 23L, 2L, 40L, 16L, 23L, 11L, 2L, 
5L, 9L, 1L, 2L, 12L, 7L, 5L, 8L, 10L, 6L, 5L), percentage = c(75, 
83.3333333333333, 100, 50, 80, 84.2105263157895, 88.2352941176471, 
57.5757575757576, 45.0980392156863, 100, 95.2380952380952, 88.8888888888889, 
79.3103448275862, 91.6666666666667, 100, 83.3333333333333, 100, 
100, 100, 85.7142857142857, 58.3333333333333, 100, 66.6666666666667, 
76.9230769230769, 100, 100), total = c(4, 6, 2, 2, 10, 19, 34, 
33, 51, 2, 42, 18, 29, 12, 2, 6, 9, 1, 2, 14, 12, 5, 12, 13, 
6, 5)), class = "data.frame", row.names = c(NA, -26L))

以上数据集以Year为因子。下面是我的输出:

summary(aov(percentage ~ Year, data = dat)) # not significant
#            Df Sum Sq Mean Sq
#Year        25   7030   281.2

有什么想法会有帮助的!

Answer 1

由于这是生物数据，因此将每年视为个性化数据很重要，因此我们认为一个因素将是最好的方法。

head(dat)
#  Year no yes percentage total
#1 1994  1   3   75.00000     4
#2 1995  1   5   83.33333     6
#3 1996  0   2  100.00000     2
#4 1997  1   1   50.00000     2
#5 1998  2   8   80.00000    10
#6 1999  3  16   84.21053    19

这里有很多问题。我将讨论:

• 线性模型与lmpercentage ~ Year/aov;

• logistic回归cbind(yes, no) ~ Yearwithglm.

使用lm/aov有什么问题

从统计学上讲，使用线性回归对手动计算的百分比("是"的百分比)进行建模肯定不是一个好主意。但是这里有一个额外的问题:你每年只有一次观测。将Year作为一个因素当然是不合理的。这样，你就有了和数据一样多的回归系数，这样你就能得到一个所有残差都为0的完美拟合。因此，所有的检验统计量和p值将是NA或NaN。

lmfit <- aov(percentage ~ Year, data = dat)
## use `summary.lm()` for aov() fit to show coefficient table
summary.lm(fit)
#Call:
#aov(formula = percentage ~ Year, data = dat)
#
#Residuals:
#ALL 26 residuals are 0: no residual degrees of freedom!  ## dang!!!
#
#Coefficients:
#            Estimate Std. Error t value Pr(>|t|)
#(Intercept)   75.000        NaN     NaN      NaN
#Year1995       8.333        NaN     NaN      NaN
#Year1996      25.000        NaN     NaN      NaN
#Year1997     -25.000        NaN     NaN      NaN
#Year1998       5.000        NaN     NaN      NaN
#Year1999       9.211        NaN     NaN      NaN
#... trimmed; all NaN ...
#
#Residual standard error: NaN on 0 degrees of freedom
#Multiple R-squared:      1,    Adjusted R-squared:    NaN 
#F-statistic:   NaN on 25 and 0 DF,  p-value: NA
anova(lmfit)
#Analysis of Variance Table
#
#Response: percentage
#          Df Sum Sq Mean Sq F value Pr(>F)
#Year      25 7029.5  281.18     NaN    NaN
#Residuals  0    0.0     NaN               
#Warning message:
#In anova.lm(lmfit) :
#  ANOVA F-tests on an essentially perfect fit are unreliable  ## dang!!!

切换到逻辑回归与glm

原则上，我们想要下面的逻辑回归。

glmfit1 <- glm(cbind(yes, no) ~ Year, family = binomial(), data = dat)

然而，每个Year仍然只有一个观察值，所以你仍然得到一个完美的拟合。在这种情况下，偏差残差均为0。

anova(glmfit1, test = "Chisq")
#Analysis of Deviance Table
#
#Model: binomial, link: logit
#
#Response: cbind(yes, no)
#
#Terms added sequentially (first to last)
#
#
#     Df Deviance Resid. Df Resid. Dev Pr(>Chi)    
#NULL                    25     70.107             
#Year 25   70.107         0      0.000 3.71e-06 ***
#---
#Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

我不确定这里的p值是否有效，因为我们有0剩余自由度(即完美拟合)。

一种解决方法(或者可能是作弊)是将数据转换为二进制格式。

Year_no <- with(dat, rep(Year, no))
Year_yes <- with(dat, rep(Year, yes))
fctr <- rep(c("no", "yes"), c(length(Year_no), length(Year_yes)))
fctr <- factor(fctr, levels = c("no", "yes"))
Year <- c(Year_no, Year_yes)
cheat <- data.frame(fctr = fctr, Year = Year)
rm(Year_no, Year_yes, fctr, Year)
head(cheat)
#  fctr Year
#1   no 1994
#2   no 1995
#3   no 1997
#4   no 1998
#5   no 1998
#6   no 1999

在此格式中，样本量上升到351，允许您将Year作为一个因素，而不会最终得到完美的拟合。

glmfit2 <- glm(fctr ~ Year, family = binomial(), data = cheat)

我认为glm的这两个规格应该是一样的，但实际上，它们不是。

all.equal(glmfit1$coef, glmfit2$coef)
#[1] "Mean relative difference: 0.3250444"

注意两个模型拟合都是收敛的。

glmfit1$converged
#[1] TRUE
glmfit2$converged
#[1] TRUE

所以我认为这真的很奇怪(可以在https://stats.stackexchange.com/上提出一个好问题)。但无论如何，这是方差分析表。

anova(glmfit2, test = "Chisq")
#Analysis of Deviance Table
#
#Model: binomial, link: logit
#
#Response: fctr
#
#Terms added sequentially (first to last)
#
#
#     Df Deviance Resid. Df Resid. Dev Pr(>Chi)    
#NULL                   350     376.80             
#Year 25   70.107       325     306.69 3.71e-06 ***
#---
#Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

在这种情况下，剩余自由度为325。我倾向于相信这个结果。

关闭评论这是我可以帮助堆栈溢出。下一步应该在https://stats.stackexchange.com/上寻找这两个逻辑回归的解释。请在问题链接中与我分享，在你发布问题之后。