Re: Calculating and using uniterrors in ASREML-R

Thanks Brian....

Don't worry - I am not using weights from BLUPs analyses....only from BLUEs.... shouldnt have scared people with the code below...

But thanks for the code. Shame that "samm.Ginverse" hasn't been reincarnated to save us from some list code....

Am still tempted to use weights as:

weights <- (pea05.pvs$predictions$'EXPT:NAME'$pvals$standard.error^2)

s

________________________________
From: ASReml users discussion group [mailto:ASREML-L_at_DPI.NSW.GOV.AU] On Behalf Of brian.cullis_at_DPI.NSW.GOV.AU
Sent: Friday, 31 July 2009 10:08 AM
To: ASREML-L_at_DPI.NSW.GOV.AU
Subject: 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 - 16:26:19 EST