Re: Calculating and using uniterrors in ASREML-R

Dear scott
ASREML in GenStat allows use of 'uniterror' as part of 2nd stage type
analyses. From what I understand these uniterror ('weights') equal:
 
1/v
 
where
 
v = the diagonal from the inverse of the Sigma (G) matrix.
 
My shortcut to do this has usually been to take the standard errors from
the Genotype predictions (and could therefore use BLUEs or BLUPs to do
it).
But should be able to do from the asreml object.
****I have real trouble with two stage with BLUPS!@ This does not make any
sense at all to me.
 ****I can cope with a short cut using ses on BLUEs - read a recent
offering by Hans Peter for all sorts of mulacky

1. I did have some code in the 'samm' days that could calculate the
ginverse.
 
    How can I get the uniterrors (for Genotypes) from some of the
following types of models?
 
    b.asreml <- asreml(msg ~ 1, random=~G + sbloc + sbloc:bloc, data=b)
 
    b2.asreml <- asreml(msg ~ G, random=~ sbloc + sbloc:bloc, data=b)
    b3.asreml.corh <- asreml(msg~E, random=~corh(E):G + at(E):sbloc +
at(E):sbloc:bloc, data=all, predict=predict.asreml(classify=list("E:G")))
****example code for a MET two stage scam
****note some of the predict syntax may be slightly different with ASReml
3.0

as this is 2.0 code

##################
# 18 sites, all over the shot
# analyse as a MET
#################
require(asreml)
pea05 <- pea05[order(pea05$EXPT,pea05$RANGE,pea05$ROW),]
pea05.asr <- asreml(yield~EXPT+EXPT:NAME,random=~at(EXPT):BLOC,
 rcov=~at(EXPT):ar1(RANGE):ar1(ROW),data=pea05,workspace=12e7)
plot(pea05.asr,forceMET="together")
pea05[pea05$EXPT=='PS3WW05' & pea05$RANGE==1 & pea05$ROW==45,'yield']<- NA
pea05[pea05$EXPT=='PS3MC05' & pea05$RANGE==10 & pea05$ROW==17,'yield']<-
NA
pea05.asr <-
  asreml(yield~EXPT+EXPT:NAME+at(EXPT,c(1)):lin(RANGE) +
         at(EXPT,c(7,8)):lin(ROW),
         random=~at(EXPT):BLOC+at(EXPT,c(1,2,9)):ROW
         + at(EXPT,c(5,9,12,13,14,15,17)):RANGE,
 rcov=~at(EXPT):ar1(RANGE):ar1(ROW),data=pea05,workspace=50e6)
plot(pea05.asr,forceMET="together")
#################
# now for the predictions
#
###############
pea05.pvs <-
  predict(pea05.asr,classify='EXPT:NAME',vcov=list("EXPT:NAME"=T),
          control=asreml.control(workspace=50e6,pworkspace=50e6))

pea05.pvs <-
  predict(pea05.asr,classify='EXPT:NAME',vcov=list("EXPT:NAME"=T),
          levels=list('EXPT:NAME'=list(EXPT="PS3BA05")))
# now we have the data frame of pvals and we add to this two extra
# columns
# ie weights and true reps
#
blues.pea05 <- pea05.pvs$predictions$'EXPT:NAME'$pvals
true.reps <-
  tapply(pea05$yield,list(pea05$NAME,pea05$EXPT),function(x)
         {length(x[!is.na(x)])})
wts <- pea05.pvs$predictions$'EXPT:NAME'$vpv
ns <- length(levels(pea05$EXPT))
nv <- length(levels(pea05$NAME))
wts.lst <- list()
for(i in 1:ns){
  st.col <- 1 + (i-1)*nv
  end.col <- i*nv
  wts.lst[[i]] <- wts[st.col:end.col,st.col:end.col]
  wts.lst[[i]] <-
wts.lst[[i]][!is.na(true.reps[,i]),!is.na(true.reps[,i])]
  wts.lst[[i]] <- diag(solve(wts.lst[[i]]))
}
blues.pea05$truereps <- as.vector(true.reps)
blues.pea05 <- blues.pea05[!is.na(blues.pea05$truereps),]
blues.pea05$weights <- unlist(wts.lst)

########################################
# end of code
#
########################################

#########################
# two stage MET
# with weights
#
#########################
fieldpeas.asr <- asreml(Yield~EXPT,random=~ us(link(~smy)):Genotype +
      Genotype:Year +
      Genotype:State + Genotype:Region + Genotype:Location +
      Genotype:Year:State +
      Genotype:Year:Region -
      Genotype:Year:Location +
      units,
      family=asreml.gaussian(dispersion=1),weights=wt,
                     G.param='fieldpeas.sv1.csv',
                     R.param='fieldpeas.sv1.csv',
 data=fieldpeas.met,control=asreml.control(workspace=30e6,maxiter=15))

2. How do I use these errors in 2nd stage analyses? (i.e. if data in 3rd
model was BLUEs and weights from a series of analyses with the 2nd
model....)

warm regards

Brian Cullis
Research Leader, Biometrics &
Senior Principal Research Scientist
NSW Department of Primary Industries
Wagga Wagga Agricultural Institute

Visiting Professorial Fellow
School of Mathematics and Applied Statistics
Faculty of Informatics
University of Wollongong

Professor,
Faculty of Agriculture, Food & Natural Resources
The University of Sydney

Adjunct Professor
School of Computing and Mathematics
Charles Sturt University

Phone: 61 2 6938 1855
Fax: 61 2 6938 1809
Mobile: 0439 448 591
 
Received on Sat Jul 31 2009 - 10:08:07 EST