数据：

wang9.1.csv

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)