|||
作者:林元震( https://yzhlinscau2018.netlify.com)
说明:由于科学网对代码的支持不好,感兴趣的读者建议浏览我的netlify官方博客。
y=xb+zb+e
Var(ξ)=σ2ξ[Σc(ρc)Σr(ρr)]
1ρcρ2c…ρnc1cρc1ρc…ρ2cρc1…ρnc1c…1
1ρrρ2r…ρnr1rρr1ρr…ρ2rρr1…ρnr1r…1
Var(e)=σ2ξ[Σc(ρc)Σr(ρr)]+σ2ηI
library(breedR)
library(breedRPlus)
library(ggplot2)
data(stroup.nin,package='agridat')
dat <- stroup.nin
names(dat)
## [1] "gen" "rep" "yield" "col" "row"str(dat)
## 'data.frame': 242 obs. of 5 variables:
## $ gen : Factor w/ 56 levels "Arapahoe","Brule",..: 12 12 12 12 12 12 12 12 12 12 ...
## $ rep : Factor w/ 4 levels "R1","R2","R3",..: NA NA NA NA NA NA NA NA NA NA ...
## $ yield: num NA NA NA NA NA NA NA NA NA NA ...
## $ col : int 1 2 3 4 5 6 7 8 9 10 ...
## $ row : int 1 1 1 1 1 1 1 1 1 1 ...
ggplot(dat, aes(yield)) + geom_histogram(col='blue',fill='white')
ggplot(dat, aes(row, col, fill = yield)) +
geom_tile(show.legend = TRUE) +
coord_fixed() + coord_flip()+
scale_fill_gradientn(colours = topo.colors(50))
RCB.bdR <- remlf90(fixed = yield ~ 1+gen, random = ~ rep, data = dat)
breedRPlus::var(RCB.bdR) ## gamma component std.error z.ratio constraint ## rep 0.1993244 9.8829 8.7928 1.123976 Positive ## Residual 1.0000000 49.5820 5.4588 9.082949 Positive summary(RCB.bdR)$model.fit[c(1,3)] ## AIC logLik ## 1253.938 -624.9691
RCB1.bdR <- remlf90(fixed = yield ~ 1+gen, spatial = list(model = 'blocks', coord = dat[, c('row','col')], id = "rep"), data = dat)
## gamma component std.error z.ratio constraint ## spatial 0.1993244 9.8829 8.7928 1.123976 Positive ## Residual 1.0000000 49.5820 5.4588 9.082949 Positive ## AIC logLik ## 1253.938 -624.9691
SUM.bdR <- remlf90(fixed = yield ~ 1+gen, spatial = list(model = 'AR', coord = dat[,c('row','col')], rho =c(.66,.44)), data = dat)
## gamma component std.error z.ratio constraint ## spatial 6.459346 35.7570 6.8376 5.229466 Positive ## Residual 1.000000 5.5357 2.5391 2.180182 Positive ## AIC logLik ## 1147.881 -571.9406
plot(SUM.bdR, type = 'spatial')+ coord_flip()+
scale_fill_gradientn(colours = topo.colors(50))
ColAR.bdR <- remlf90(fixed = yield ~ 1+gen, spatial = list(model = 'AR', coord = dat[,c('row','col')], rho = c(0,.44)), data = dat)
## gamma component std.error z.ratio constraint ## spatial 26138.25 43.7110000 4.7691000 9.1654610 Positive ## Residual 1.00 0.0016723 0.0082369 0.2030254 Positive ## AIC logLik ## 1198.539 -597.2693
#plot(ColAR.bdR, type = 'spatial')
RowAR.bdR <- remlf90(fixed = yield ~ 1+gen, spatial = list(model = 'AR', coord = dat[,c('row','col')], rho = c(0.66,0)), data = dat)
## gamma component std.error z.ratio constraint ## spatial 20272.6 57.0390000 6.22540 9.16230282 Positive ## Residual 1.0 0.0028136 0.06103 0.04610192 Positive ## AIC logLik ## 1198.989 -597.4946
#plot(RowAR.bdR, type = 'spatial')
注意:RCB1模型是把区组放到空间分析方法里,因此所计算的空间相关误差其实就是区组方差。
compare.plots(
list(RCB = plot(RCB1.bdR, type = 'fullspatial'),
SUM = plot(SUM.bdR, type = 'fullspatial'),
ColAR = plot(ColAR.bdR, type = 'fullspatial'),
RowAR = plot(RowAR.bdR, type = 'fullspatial')
)) + coord_flip() +
scale_fill_gradientn(colours = topo.colors(50))
sres<-function(object){
df<-summary(object)$model.fit[c(1,3)]
df[,c('Spatial','Residual')]=var(object)[,2] return(df)
}
M1<-sres(RCB1.bdR)
M2<-sres(SUM.bdR)
M3<-sres(ColAR.bdR)
M4<-sres(RowAR.bdR)
res<-rbind (M1,M2,M3,M4)
row.names(res)<-c('RCB','SUM','ColAR','RowAR')
res
Model | AIC | logLik | Spatial | Residual |
---|---|---|---|---|
RCB | 1253.938 | -624.969 | 9.883 | 49.582 |
SUM | 1147.881 | -571.941 | 35.757 | 5.536 |
ColAR | 1198.539 | -597.269 | 43.711 | 0.002 |
RowAR | 1198.989 | -597.495 | 57.039 | 0.003 |
mvres<-function(object){
df<-data.frame(mv=fixef(object)$gen,mv.se=attr(fixef(object)$gen,'se'))
df$ID<-row.names(df)
df<-dplyr::arrange(df,desc(mv))
df$mv.rank<-1:nrow(df) return(df)
}
fM1<-mvres(RCB.bdR)
names(fM1)<-c('m1','m1.se','ID','m1.rank')
fM2<-mvres(SUM.bdR)
names(fM2)<-c('m2','m2.se','ID','m2.rank')
fM3<-mvres(ColAR.bdR)
names(fM3)<-c('m3','m3.se','ID','m3.rank')
fM4<-mvres(RowAR.bdR)
names(fM4)<-c('m4','m4.se','ID','m4.rank')
fM.res11<-merge(fM1,fM2,by='ID')
fM.res12<-merge(fM3,fM4,by='ID')
fM.res<-merge(fM.res11,fM.res12,by='ID')
fM.res1<-dplyr::arrange(fM.res,desc(m1))
fM.res2<-dplyr::arrange(fM.res,desc(m2))
head(fM.res1,8)
ID | m1 | m1.se | m1.rank | m2 | m2.se | m2.rank | m3 | m3.se | m3.rank | m4 | m4.se | m4.rank |
---|---|---|---|---|---|---|---|---|---|---|---|---|
NE86503 | 32.65 | 3.856 | 1 | 26.49 | 2.473 | 14 | 29.33 | 2.846 | 5 | 27.56 | 2.632 | 11 |
NE87619 | 31.26 | 3.856 | 2 | 29.43 | 2.552 | 2 | 29.84 | 2.849 | 2 | 30.00 | 2.720 | 2 |
NE86501 | 30.94 | 3.856 | 3 | 25.13 | 2.554 | 27 | 28.42 | 2.904 | 11 | 25.56 | 2.752 | 22 |
Redland | 30.50 | 3.856 | 4 | 28.41 | 2.574 | 5 | 29.54 | 2.891 | 3 | 28.45 | 2.820 | 6 |
Centurk78 | 30.30 | 3.856 | 5 | 26.39 | 2.584 | 16 | 28.07 | 2.897 | 13 | 26.60 | 2.724 | 12 |
NE83498 | 30.12 | 3.856 | 6 | 29.34 | 2.609 | 3 | 29.41 | 2.882 | 4 | 29.85 | 2.850 | 3 |
Siouxland | 30.11 | 3.856 | 7 | 23.89 | 2.526 | 37 | 28.01 | 2.895 | 14 | 24.11 | 2.647 | 36 |
NE86606 | 29.76 | 3.856 | 8 | 27.43 | 2.503 | 9 | 28.42 | 2.850 | 10 | 28.12 | 2.658 | 7 |
head(fM.res2,8)
ID | m1 | m1.se | m1.rank | m2 | m2.se | m2.rank | m3 | m3.se | m3.rank | m4 | m4.se | m4.rank |
---|---|---|---|---|---|---|---|---|---|---|---|---|
Buckskin | 25.56 | 3.856 | 28 | 33.78 | 2.531 | 1 | 31.61 | 2.941 | 1 | 32.11 | 2.710 | 1 |
NE87619 | 31.26 | 3.856 | 2 | 29.43 | 2.552 | 2 | 29.84 | 2.849 | 2 | 30.00 | 2.720 | 2 |
NE83498 | 30.12 | 3.856 | 6 | 29.34 | 2.609 | 3 | 29.41 | 2.882 | 4 | 29.85 | 2.850 | 3 |
NE85556 | 26.39 | 3.856 | 23 | 28.98 | 2.559 | 4 | 28.82 | 2.844 | 8 | 28.76 | 2.732 | 5 |
Redland | 30.50 | 3.856 | 4 | 28.41 | 2.574 | 5 | 29.54 | 2.891 | 3 | 28.45 | 2.820 | 6 |
Brule | 26.07 | 3.856 | 25 | 28.08 | 2.553 | 6 | 26.28 | 2.888 | 22 | 29.20 | 2.734 | 4 |
NE86507 | 23.79 | 3.856 | 39 | 27.67 | 2.473 | 7 | 26.36 | 2.844 | 20 | 27.77 | 2.654 | 9 |
Scout66 | 27.52 | 3.856 | 17 | 27.53 | 2.554 | 8 | 28.82 | 2.843 | 7 | 25.91 | 2.777 | 17 |
Archiver|手机版|科学网 ( 京ICP备07017567号-12 )
GMT+8, 2024-10-19 23:03
Powered by ScienceNet.cn
Copyright © 2007- 中国科学报社