#asreml #error asr_binomial() Er
2019-07-11 本文已影响67人
董八七
Invalid Binomial/mltinomial proportion ,2.000000/1.000000 in record 559Error in asreml(str ~ 1, random = ~Mum, family = asr_binomial(), subset = Spacing == : Errors in GLM.
原因是响应变量不是二项分布，而是有多个分类。这是要用将响应变量转因子后，用asr_multinomial()。
林元震老师教材的第二版中10.4.2.5（p485）中提到泊松分布，但所用的响应变量str应该是干直度的分级，属于多远分类变量而非泊松分布，实际的运行效果如下：
- 泊松分布
bm_asr <- asreml(str~1,
random = ~Mum,
family = asr_poisson(),
subset = Spacing==3,
data = dfm2)
# Model fitted using the sigma parameterization.
# ASReml4 Beta-release 4.0.0.9005 Tue Jul 11 09:16:32 2019
# Poisson; Log Mu=exp(XB) V=Mu;
# Note: The LogLik value is unsuitable for comparing GLM models
#
# LogLik Sigma2 DF wall cpu
# 1 -89.0416 1.0 558 09:16:32 0.0 (1 restrained)
# 2 -66.0837 1.0 558 09:16:32 0.0
# 3 -71.2797 1.0 558 09:16:32 0.0
# 4 -71.4648 1.0 558 09:16:32 0.0
# 5 -71.5269 1.0 558 09:16:32 0.0
# 6 -71.5594 1.0 558 09:16:32 0.0
# 7 -71.5616 1.0 558 09:16:32 0.0
# 8 -71.5617 1.0 558 09:16:32 0.0
# Deviance from GLM fit: 802.73
# Variance heterogenity factor (Deviance/df): 1.44
# (assuming 558 degrees of freedom)
summary(bm_asr)$varcomp
# component std.error z.ratio bound %ch
# Mum 0.01318554 0.01032939 1.276507 P 0
# units(R) 1.00000000 NA NA F 0
plot(bm_asr)
h2 <- vpredict(bm_asr, h2~4*V1/(V1+V2))
# h2 h2.se
# 0.052 0.040
泊松分布残差图
- 多项分类
dfm2$str <- dfm2$str %>% factor
bm_asr <- asreml(str~trait, #注意这里是trait，不是1
random = ~Mum,
family = asr_multinomial(),
subset = Spacing==3,
data = dfm2)
summary(bm_asr)$varcomp
# Model fitted using the sigma parameterization.
# ASReml4 Beta-release 4.0.0.9005 Tue Jul 11 09:14:16 2019
# Cum. Multinomial; Logit P=1/(1+exp(-XB));
# Note: The LogLik value is unsuitable for comparing GLM models
#
# LogLik Sigma2 DF wall cpu
# 1 -1388.772 1.0 2790 09:14:16 0.0
# 2 -1391.806 1.0 2790 09:14:16 0.0
# 3 -1393.413 1.0 2790 09:14:16 0.0
# 4 -1394.429 1.0 2790 09:14:16 0.0
# 5 -1395.058 1.0 2790 09:14:16 0.0
# 6 -1395.444 1.0 2790 09:14:16 0.0
# 7 -1395.913 1.0 2790 09:14:16 0.0
# 8 -1396.005 1.0 2790 09:14:16 0.0
# 9 -1396.034 1.0 2790 09:14:16 0.0
# 10 -1396.041 1.0 2790 09:14:16 0.0
# 11 -1396.043 1.0 2790 09:14:16 0.0
# Deviance from GLM fit: 1986.22
# Variance heterogenity factor (Deviance/df): 0.71
# (assuming 2790 degrees of freedom)
summary(bm_asr)$varcomp
# component std.error z.ratio bound %ch
# Mum 0.04835465 0.06831176 0.7078524 P 0
# units:trait(R) 1.00000000 NA NA F 0
# units:trait!trait_0:0 1.00000000 NA NA F 0
# units:trait!trait_1:0 0.00000000 NA NA F NA
# units:trait!trait_1:1 1.00000000 NA NA F 0
# units:trait!trait_2:0 0.00000000 NA NA F NA
# units:trait!trait_2:1 0.00000000 NA NA F NA
# units:trait!trait_2:2 1.00000000 NA NA F 0
# units:trait!trait_3:0 0.00000000 NA NA F NA
# units:trait!trait_3:1 0.00000000 NA NA F NA
# units:trait!trait_3:2 0.00000000 NA NA F NA
# units:trait!trait_3:3 1.00000000 NA NA F 0
# units:trait!trait_4:0 0.00000000 NA NA F NA
# units:trait!trait_4:1 0.00000000 NA NA F NA
# units:trait!trait_4:2 0.00000000 NA NA F NA
# units:trait!trait_4:3 0.00000000 NA NA F NA
# units:trait!trait_4:4 1.00000000 NA NA F 0
plot(bm_asr)
# 出这个图时，会报错Error in log(y[nz]) : non-numeric argument to mathematical function，
# 在原函数plot.asreml中，将rr <- resid(x, type = res, spatial = spatial)改成rr <- resid(x, spatial = spatial)
h2 <- vpredict(bm_asr, h2~4*V1/(V1+V2*3.29))
# h2 h2.se
# 0.058 0.081
# 遗传力估计值和上面模型是近似的，但标准误差得比较大。
多项分类残差图
这个图怎么解读呢？？？
