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.
原因是響應(yīng)變量不是二項(xiàng)分布，而是有多個(gè)分類。這是要用將響應(yīng)變量轉(zhuǎn)因子后，用asr_multinomial()。
林元震老師教材的第二版中10.4.2.5（p485）中提到泊松分布，但所用的響應(yīng)變量str應(yīng)該是干直度的分級(jí)，屬于多元分類變量而非泊松分布，實(shí)際的運(yùn)行效果如下：
- 泊松分布
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
泊松分布?xì)埐顖D
- 多項(xiàng)分類
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)
# 出這個(gè)圖時(shí)，會(huì)報(bào)錯(cuò)Error in log(y[nz]) : non-numeric argument to mathematical function，
# 在原函數(shù)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
# 遺傳力估計(jì)值和上面模型是近似的，但標(biāo)準(zhǔn)誤差得比較大。
多項(xiàng)分類殘差圖
這個(gè)圖怎么解讀呢？？？
最新的版本殘差圖是這樣的：
image.png
也需是前面修改plot的代碼導(dǎo)致殘差值提取錯(cuò)誤導(dǎo)致的。
