Nilaya,

You can use lmer and calculate the fitted values, then average over the
predictions as appropriate.

The type of predictions you want are easy with the R version of asreml.
Note, this is non-free software, but worth its cost if you have complex
problems.  See http://www.vsni.co.uk/software/asreml/
See below.

Both versions are below.

Kevin Wright


R> require(lme4)
R> require(SASmixed)
R> data(Animal)
R> anim <- Animal
R> names(anim) <- c('sire','dam','adg')

m2 <- lmer(adg ~ 1 + (1|sire) + (1|dam:sire), data=anim)
anim$fit <- fitted(m2)

# Broad - don't include dam:sire
R> coef(m2)$sire
  (Intercept)
1    2.205746
2    2.271907
3    2.505892
4    2.129903
5    2.536552

# Narrow.  Average these specific dam:sire random effects
R> aggregate(fit~sire, data=anim, mean)
  sire      fit
1    1 2.160931
2    2 2.250955
3    3 2.569332
4    4 2.057733
5    5 2.611050

# dam 1 sire 1
R> anim[anim$sire==1 & anim$dam==1,]
  sire dam  adg      fit
1    1   1 2.24 2.100192
2    1   1 1.85 2.100192



R> require(asreml)
R> m1 <- asreml(adg~ 1, data=anim, random= ~ sire + dam:sire)
# Broad space predictions
R> predict(m1, classify='sire')$predictions$pvals

Notes:
- The predictions are obtained by averaging across the hypertable
  calculated from model terms constructed solely from factors in
  the averaging and classify sets.
- Use "average" to move ignored factors into the averaging set.

- The ignored set:  dam

  sire predicted.value standard.error est.status
1    1        2.205747      0.1420762  Estimable
2    2        2.271908      0.1420762  Estimable
3    3        2.505890      0.1420762  Estimable
4    4        2.129905      0.1420762  Estimable
5    5        2.536550      0.1420762  Estimable

# Narrow space
R> predict(m1, classify='sire', average=list(dam=NULL))$predictions$pvals

Notes:
- The predictions are obtained by averaging across the hypertable
  calculated from model terms constructed solely from factors in
  the averaging and classify sets.
- Use "average" to move ignored factors into the averaging set.

- The SIMPLE averaging set:  dam

  sire predicted.value standard.error est.status
1    1        2.160933     0.09344944  Estimable
2    2        2.250956     0.09344944  Estimable
3    3        2.569328     0.09344944  Estimable
4    4        2.057737     0.09344944  Estimable
5    5        2.611046     0.09344944  Estimable

# Dam 1 Sire 1
R> predict(m1, classify='sire:dam')$predictions$pvals


Notes:
- The predictions are obtained by averaging across the hypertable
  calculated from model terms constructed solely from factors in
  the averaging and classify sets.
- Use "average" to move ignored factors into the averaging set.


   dam sire predicted.value standard.error est.status
1    1    1        2.100194      0.1228243  Estimable
2    1    2        2.067096      0.1228243  Estimable
3    1    3        2.505306      0.1228243  Estimable
4    1    4        2.185782      0.1228243  Estimable
5    1    5        2.663577      0.1228243  Estimable
6    2    1        2.221673      0.1228243  Estimable
7    2    2        2.434815      0.1228243  Estimable
8    2    3        2.633351      0.1228243  Estimable
9    2    4        1.929692      0.1228243  Estimable
10   2    5        2.558515      0.1228243  Estimable




> --
> | David Duffy (MBBS PhD)                                         ,-_|\
> | email: davidD@qimr.edu.au  ph: INT+61+7+3362-0217 fax: -0101  /     *
> | Epidemiology Unit, Queensland Institute of Medical Research   \_,-._/
> | 300 Herston Rd, Brisbane, Queensland 4029, Australia  GPG 4D0B994A v
>




On Thu, Apr 21, 2011 at 2:52 AM, David Duffy <davidD@qimr.edu.au> wrote:

> On Mon, 18 Apr 2011, Nilaya Sharma wrote:
>
>  Dear Mixed Model Experts
>>
>> I was trying to fit a mixed model in animal experiment but stuck at simple
>> point. The following similar example is from SAS mixed model pp 212.
>>
>> require(lme4); mod <- lmer(adg ~ 1 + (1|sire) + (1|dam:sire),
>> data=genetic_evaluation)
>> ranef (mod)
>>
>> provides random effects associated with dam:sire and sire,  but do not
>> provide the comparisons I  need. It looks like I need to define and apply
>> own contrasts or there is any other functions to do so?
>>
>
>  ****************Expected outputs according to SAS
>>                     Estimate
>> sire 1 BLUP "broad    2.2037
>> sire 1 BLUP "narrow"  2.1609
>> sire 1 BLUP with dam1 2.1002
>>
>
> If you look at mod and the output from ranef(mod), you will see that you
> have all the necessary coefficients to form the BLUPs as defined in Littell
> et al., recalling that they are centred.  It's easy to write your own
> routines to automate this if you have a large application.
>
> Cheers, David Duffy.
>
> _______________________________________________
> R-sig-mixed-models@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>

	[[alternative HTML version deleted]]


