2017-02-20 72 views
1

我在理解最小二乘均值的行为时遇到问题。以下是一个玩具示例,它使用随机数据集来演示我的问题。该方案是这样的:存在被1999年至2015年最小二乘方表示

# Number of observations in data set 
n.obs <- 1000 

# Create dummy data set 
df.tst <- data.frame(density = runif(n.obs, 0, 1), 
       year = as.factor(round(runif(n.obs, 1999, 2015))), 
       season = sample(c("Fall", "Spring"), n.obs, replace= TRUE), 
       station = sample(paste("Stat", 1:10), n.obs, replace= TRUE)) 

这里之间取样用于在任一春季或秋季的度量称为密度10个站,我已经随机分配采样特意表示片状数据集。由于数据集是不完整的,因此我使用最小二乘平均值的方法,以便可以在避免采样偏差影响的同时估计一个季节和一年的密度。

# Fit linear model to data 
fitted.model <- lm(density ~ year + season + station, data=df.tst) 

# Calculate least square means 
seasonal.model <- summary(lsmeans(fitted.model, c("year", "season"))) 

比较的结果,我创建了我称之为“异常”,也就是某年/季减均值为所有年(以中心值)和归一化的标准偏差值的度量。这个想法是,这提供了一些标准化的度量标准,比如2005年的春天与其他春季的观测值有所不同。所以...

# Anomaly for spring 
seasonal.model %>% 
    filter(season == "Spring") %>% 
    mutate(anom = (lsmean - mean(lsmean))/sd(lsmean)) 

这是伟大的。我不明白的是,为什么每当我这样做时,我都会得到同样的答案。例如...

# Anomaly for fall 
seasonal.model %>% 
    filter(season == "Fall") %>% 
    mutate(anom = (lsmean - mean(lsmean))/sd(lsmean)) 

这些给,

year season lsmean   SE df lower.CL upper.CL   anom 
1 1999 Spring 0.5966423 0.05242108 973 0.4937709 0.6995137 1.879970679 
2 2000 Spring 0.4510249 0.03717688 973 0.3780688 0.5239810 -1.617190892 
3 2001 Spring 0.4683691 0.03713725 973 0.3954908 0.5412474 -1.200649943 
4 2002 Spring 0.4451014 0.03730148 973 0.3719008 0.5183020 -1.759450098 
5 2003 Spring 0.5035975 0.03881258 973 0.4274315 0.5797635 -0.354600778 
6 2004 Spring 0.5263538 0.03505567 973 0.4575604 0.5951472 0.191916022 
7 2005 Spring 0.5553849 0.03703036 973 0.4827163 0.6280535 0.889129265 
8 2006 Spring 0.5182714 0.03626301 973 0.4471087 0.5894341 -0.002191364 
9 2007 Spring 0.4765210 0.04097960 973 0.3961024 0.5569395 -1.004874623 
10 2008 Spring 0.5465877 0.04483499 973 0.4586033 0.6345721 0.677854169 
11 2009 Spring 0.4959347 0.03185768 973 0.4334170 0.5584523 -0.538633207 
12 2010 Spring 0.5359122 0.03934735 973 0.4586968 0.6131276 0.421471255 
13 2011 Spring 0.5533493 0.03461044 973 0.4854296 0.6212690 0.840241309 
14 2012 Spring 0.5465454 0.03697864 973 0.4739783 0.6191124 0.676838641 
15 2013 Spring 0.5594985 0.03436268 973 0.4920650 0.6269320 0.987923047 
16 2014 Spring 0.5320895 0.03825361 973 0.4570205 0.6071586 0.329666086 
17 2015 Spring 0.5009818 0.04771979 973 0.4073363 0.5946274 -0.417419568 

year season lsmean   SE df lower.CL upper.CL   anom 
1 1999 Fall 0.5683228 0.05226823 973 0.4657514 0.6708943 1.879970679 
2 2000 Fall 0.4227054 0.03638423 973 0.3513048 0.4941060 -1.617190892 
3 2001 Fall 0.4400496 0.03752816 973 0.3664042 0.5136951 -1.200649943 
4 2002 Fall 0.4167819 0.03741172 973 0.3433649 0.4901989 -1.759450098 
5 2003 Fall 0.4752781 0.04039258 973 0.3960115 0.5545447 -0.354600778 
6 2004 Fall 0.4980343 0.03492573 973 0.4294959 0.5665728 0.191916022 
7 2005 Fall 0.5270654 0.03591814 973 0.4565795 0.5975514 0.889129265 
8 2006 Fall 0.4899519 0.03618696 973 0.4189385 0.5609654 -0.002191364 
9 2007 Fall 0.4482015 0.04026627 973 0.3691828 0.5272202 -1.004874623 
10 2008 Fall 0.5182682 0.04545050 973 0.4290759 0.6074605 0.677854169 
11 2009 Fall 0.4676152 0.064 973 0.4046207 0.5306096 -0.538633207 
12 2010 Fall 0.5075927 0.03880628 973 0.4314391 0.5837464 0.421471255 
13 2011 Fall 0.5250298 0.03440547 973 0.4575123 0.5925473 0.840241309 
14 2012 Fall 0.5182259 0.03755452 973 0.4445287 0.5919231 0.676838641 
15 2013 Fall 0.5311791 0.03463023 973 0.4632205 0.5991376 0.987923047 
16 2014 Fall 0.5037701 0.03826525 973 0.4286782 0.5788620 0.329666086 
17 2015 Fall 0.4726624 0.04793952 973 0.3785856 0.5667391 -0.417419568 

我希望我误解的东西很简单,但为什么春季和秋季异常是完全一样的即使实际的最小二乘方法实际上是不同的。任何见解将不胜感激...

编辑: 在响应bethanyP。这是来自另一次运行,因为数据是随机的,所以可能会有所不同。

Classes ‘summary.ref.grid’ and 'data.frame': 34 obs. of 7 variables: 
$ year : Factor w/ 17 levels "1999","2000",..: 1 2 3 4 5 6 7 8 9 10 ... 
$ season : Factor w/ 2 levels "Fall","Spring": 1 1 1 1 1 1 1 1 1 1 ... 
$ lsmean : num 0.484 0.461 0.441 0.495 0.575 ... 
$ SE  : num 0.046 0.0361 0.0355 0.0408 0.0347 ... 
$ df  : num 973 973 973 973 973 973 973 973 973 973 ... 
$ lower.CL: num 0.394 0.39 0.371 0.415 0.507 ... 
$ upper.CL: num 0.574 0.532 0.511 0.575 0.643 ... 
- attr(*, "estName")= chr "lsmean" 
- attr(*, "clNames")= chr "lower.CL" "upper.CL" 
- attr(*, "pri.vars")= chr "year" "season" 
- attr(*, "mesg")= chr "Results are averaged over the levels of: station" "Confidence level used: 0.95" 
+0

你可以运行str(seasonal.model)并在此显示结果吗?我想我可以帮你找到需要的地方。 – sconfluentus

+0

@bethanyP我已经将该信息添加到原始问题中了。 – Lyngbakr

回答

0

试试这个,我想可能不会被过滤,因为你给了它一个字符串,但赛季是一个因素:

seasonal.model %>% 
    filter(season == as.factor("Spring")) %>% 
    mutate(anom = (lsmean - mean(lsmean))/sd(lsmean)) 

如果是这样的话,这应该工作,你可以设置: spring = as.factor("Spring")然后只需将spring送入您的管道。

试试这个:seasonal.model %>% group_by(year, season) %>% summarize(anom = (lsmean - mean(lsmean))/sd(lsmean))

+0

感谢您的回复,bethanyP。当我运行您的代码时,出现以下错误: filter_impl(。数据,点):水平集的因素是不同的 但是,当我运行seasonal.model%>%filter(season ==“Spring”)时,它确实过滤正确(即只给我弹簧值)。 – Lyngbakr

+0

是的,因为秋天仍然是一个字符串......有道理 – sconfluentus

+0

由于某种原因,第二块代码会产生NaN。 – Lyngbakr

1

的原因是因为你安装一个添加剂模型,这意味着你承担了年度的影响是在每个赛季一样。换句话说,预测之间的关系每个季节都完全相同。

最小二乘法意味着完全基于模型预测 - 所以模型很重要!如果你适合一个有交互作用的模型,那么你每个季节都会得到不同的异常情况。

+0

@lyngbakr你甚至看过这个答案吗? – rvl