||
数据:
dat <- read.table("wang9.1.csv",head=T,sep=",")
str(dat)
# 可以看出cul为数字型,需要转化为factor
# > str(dat)
# 'data.frame': 360 obs. of 4 variables:
# $ loc : Factor w/ 10 levels "L1","L10","L2",..: 1 1 1 1 3 3 3 3 4 4 ...
# $ cul : int 1 1 1 1 1 1 1 1 1 1 ...
# $ block: int 1 2 3 4 1 2 3 4 1 2 ...
# $ yield: num 2.21 1.92 2.44 2.19 2.4 ...
for( i in 1:3) dat[,i] <- as.factor(dat[,i]) # for 循环批量转化
str(dat)
# > str(dat)
# 'data.frame': 360 obs. of 4 variables:
# $ loc : Factor w/ 10 levels "L1","L10","L2",..: 1 1 1 1 3 3 3 3 4 4 ...
# $ cul : Factor w/ 9 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
# $ block: Factor w/ 4 levels "1","2","3","4": 1 2 3 4 1 2 3 4 1 2 ...
# $ yield: num 2.21 1.92 2.44 2.19 2.4 ...
dat <- dat[order(dat$loc),]
library(asreml)
summary(dat)
# 模型一,假定地点间方差其次
model1 <- asreml(yield ~1 + cul*loc+loc/block,
rcov =~ units,data=dat)
# 模型二,假定地点间方差不齐次,然后检测两个模型的显著性LRT检验,如果显著,说明方差不齐次,不能# # 进行联合方差分析
model2 <- asreml(yield ~1 + cul*loc+loc/block,
rcov =~ at(loc):units,data=dat)
# > model1 <- asreml(yield ~1 + cul*loc+loc/block,
# + rcov =~ units,data=dat)
# ASReml: Sat Nov 05 08:55:14 2016
#
# LogLik S2 DF wall cpu
# 267.8811 0.0189 240 08:55:14 0.0
# 267.8811 0.0189 240 08:55:14 0.0
#
# Finished on: Sat Nov 05 08:55:14 2016
#
# LogLikelihood Converged
# > model2 <- asreml(yield ~1 + cul*loc+loc/block,
# + rcov =~ at(loc):units,data=dat)
# ASReml: Sat Nov 05 08:55:21 2016
#
# LogLik S2 DF wall cpu
# 258.7141 1.0000 240 08:55:21 0.0
# 270.2500 1.0000 240 08:55:21 0.0
# 274.2532 1.0000 240 08:55:21 0.0
# 275.5956 1.0000 240 08:55:21 0.0
# 275.6803 1.0000 240 08:55:21 0.0
# 275.6811 1.0000 240 08:55:21 0.0
# 275.6811 1.0000 240 08:55:21 0.0
#
# Finished on: Sat Nov 05 08:55:21 2016
#
# LogLikelihood Converged
# 手动的卡方检验,和AAfun中的LRT检验是一样的
1- pchisq(2*(model2$loglik-model1$loglik),9)
library(AAfun)
model.comp(model1,model2,LRT = T)
# 手动计算,概率为0.075,不显著,说明方差齐次
# > 1- pchisq(2*(model2$loglik-model1$loglik),9)
# [1] 0.07571481
# 利用AAfun中的model.comp检验,概率为0.076,结果和手动计算的概率一致。
# > library(AAfun)
# > model.comp(model1,model2,LRT = T)
# Attension:
# Fixed factors should be the same!
#
#
# Model LogL Npm AIC BIC AIC.State BIC.State
# 1 model1 267.8811 1 -533.7621 -530.2815 better better
# 2 model2 275.6811 10 -531.3623 -496.5559
# -----------------------------
# Lower AIC and BIC is better model.
#
#
#
# Attension: Please check every asreml results' length is 43;
# if the length < 43, put the object at the end of Nml.
# In the present, just allow one object's length < 43.
# =====================================
# Likelihood ratio test (LRT) results:
#
#
# Model compared between model1 -- model2 :
# Model LogL Npm AIC Pr(>F) Sig.level
# 1 model1 267.8811 1 -533.7621
# 2 model2 275.6811 10 -531.3623 0.076 Not signif
# ---------------
# Sig.level: 0'***' 0.001 '**' 0.01 '*' 0.05 'Not signif' 1
#
# =====================================
# Attension: Ddf did not minus 0.5.
# When for corr model, against 0.
说明地点间方差异质,可以进行联合方差分析。
如果方差不齐次,需要手动定义地点间的误差,需要用模型二:
model2 <- asreml(yield ~1 + cul*loc+loc/block,
rcov =~ at(loc):units,data=dat)
Archiver|手机版|科学网 ( 京ICP备07017567号-12 )
GMT+8, 2024-9-20 15:17
Powered by ScienceNet.cn
Copyright © 2007- 中国科学报社