今天处理数据，用
predict()获得SED，但预测值、标准误及SED都缺失。估计状态是aliased，即不可估计，上论坛查找原因。
View topic - Predicting over three seasons for many individuals | Forum | VSN International
kiwi Posted: Wed Mar 24, 2010 2:02 am
数据示例：
# Season Vine Block Flower DM
# 2007-52-04-14g 2007 52-04-14g 52 first 18.87
# 2008-52-04-14g 2008 52-04-14g 52 previously 18.64
# 2009-52-04-14g 2009 52-04-14g 52 previously 17.70
# 2007-52-04-15a 2007 52-04-15a 52 first 18.08
# 2008-52-04-15a 2008 52-04-15a 52 previously 17.98
# 2009-52-04-15a 2009 52-04-15a 52 previously 17.20
# 2007-52-04-15b 2007 52-04-15b 52 first 17.27
# 2008-52-04-15b 2008 52-04-15b 52 previously 18.52
# 2009-52-04-15b 2009 52-04-15b 52 previously 18.50
# 2007-52-04-15d 2007 52-04-15d 52 first NA
# 2008-52-04-15d 2008 52-04-15d 52 previously 18.20
#
# with(model.df, levels(Block))
# [1] "3" "39" "45" "52" "54"
# with(model.df, levels(Season))
# [1] "2007" "2008" "2009"
# with(model.df, levels(Flower))
# [1] "first" "previously"
# length(levels(model.df$Vine))
# [1] 2345
模型：
dm.asr <- try(asreml(fixed = DM ~ Season*Flower + Block,
random = ~ped(Vine):corv(Block) +
ped(Vine):idv(Season),
rcov = ~id(Vine):idv(Season),
na.method.X = "include",
ginverse = list(Vine = G.inv$ginv),
data = model.df,
workspace = 50e6, maxiter = 50))
预测
dict <- predict(dm.asr, classify = "Vine:Season", maxiter = 1)
use.preds <- dict$predictions$pvals[,1:4]
出现2个问题：
- Season和Vine的顺序和原始数据中不一致
tail(use.preds)
# Season Vine predicted.value standard.error
# 7030 2009 54-13-01b 19.17404 1.100054
# 7031 2009 54-13-04e 20.88258 1.100054
# 7032 2009 54-13-12b 19.23295 1.100054
# 7033 2009 54-14-04b 18.05464 1.100054
# 7034 2009 54-14-04e 20.76475 1.100054
# 7035 2009 54-14-04g 19.35078 1.100054
tail(model.df)
# Season Vine Block Flower DM
# 2007-54-14-04e 2007 54-14-04e <NA> <NA> NA
# 2008-54-14-04e 2008 54-14-04e <NA> <NA> NA
# 2009-54-14-04e 2009 54-14-04e 54 first 23.1
# 2007-54-14-04g 2007 54-14-04g <NA> <NA> NA
# 2008-54-14-04g 2008 54-14-04g <NA> <NA> NA
# 2009-54-14-04g 2009 54-14-04g 54 first 20.7
有点费解，但不是大问题。但是，2007年的每一个Vine的预测都是NA，而其他Season的每个Vine都预测，无论它们是否有读数（根据需要）。 如果我去掉Season 2007，那么2008年的每个预测就都会变成NA。 难道这两种现象是相关的，得到的预测值是可信的吗？ 我也想知道如何在模型的随机部分中进行三向交互。 我知道可以在空间模型中完成，但我能够收敛的唯一模型具有“singular”约束，而且z ratios极高。根据我的经验，这个模型所得到的预测值并不可信（即使残差图看起来很好）。
david.butler Posted: Wed Mar 31, 2010 4:15 am
你的模型在固定效应中有Flower x Season，然而，就你的预测而言，“问题”是这是一个不平衡的组合，Season 1 (2007)没有Flower的水平2。 这些项在平均集中并将导致aliasing，在这种情况下，所有2007年预测都将aliased【这几句是解决问题的关键】。 有几个方法可以处理，但一个例子是使用“present”平均，即对数据中存在的组合进行“平均”预测。 在asreml-R中，可以使用'present'参数指定它以进行预测。 如果你添加present=c('Season','Flower')，然后，预测将由Flower x Season列表中的成员组成。
不清楚你是如何去掉2007的，但利用levels=list(Season=c("2008","2009"))你可以只获得2008和2009的预测值。
当前版本的asreml-R不按照分类集中指定项的顺序，目前这由R遇到项的顺序决定。预测值应该是正确的。
具体到我的问题
查看预测结果，可以得到如下信息：
asr_pv$predictions$pvals
Notes:
- The predictions are obtained by averaging across the hypertable
calculated from model terms constructed solely from factors in
the averaging and classify sets.
- Use "average" to move ignored factors into the averaging set.
- The SIMPLE averaging set: Rep
- Aliased cells: 79.000000
- The overall SED includes non-estimable predictions
Clone predicted.value standard.error est.status
1 delta NA NA Aliased
2 orchard2aCK NA NA Aliased
3 RC302xCHANG238-12 NA NA Aliased
4 RC302xCHANGhun10-18 NA NA Aliased
知道，平均集是Rep，套到上面那个例子中，就是预测因子和Rep的组合不平衡，即有些水平中有，有些没有，所以加一个present=c()即可解决，尽管我的Rep是固定因子。
