Re: [breedR] Multi-trait model not working

Richard Whittet

unread,

Jun 17, 2020, 5:43:55 AM6/17/20

to H G D, breedR

Hello,

Could just be syntax..? Try replacing the "+" in the cbind with a ","

So cbind(var1, var2)

Best wishes

Richard

On Wed, 17 Jun 2020, 10:38 H G D, <hg.dal...@gmail.com> wrote:

Dear Facundo and others,

I am trying to run a Multi-trait model, according to the proposed mehtod using cbind for the traits. Unfortunately, the function does not seem to recognize one of the traits?
Here is an example, where i stripped down the number of arguments to the bare minumum:

Best Regards, Hermann

> head(dat[,c("ROW","COL","GEN","FZ","AUSWFL")])
    ROW COL   GEN FZ AUSWFL
742   1   1 19939 220      1
743   3   1 23003 62      3
744   4   1 44398 62      6
745   5   1 44398 62      5
[ reached 'max' / getOption("max.print") -- omitted 2 rows ]
> res <- remlf90(dat = dat, fixed = cbind(FZ + AUSWFL) ~ GEN) #fixed = FZ ~ GEN,
Using default initial variances given by default_initial_variance()
See ?breedR.getOption.

> summary(res)
Formula: cbind(FZ + AUSWFL) ~ 0 + GEN
   Data: dat
AIC BIC logLik
2742 2747 -1370

Variance components:
         Estimated variances S.E.
Residual                2226 170.5

Fixed effects:
            value    s.e.
GEN.19171   0.000 0.0000
GEN.19187 152.000 47.1788
GEN.19188   0.000 0.0000
GEN.19191   0.000 0.0000
GEN.19456 67.000 33.3604
GEN.19493 99.750 23.5894
GEN.19495   0.000 0.0000
GEN.19496   0.000 0.0000
GEN.19497   0.000 0.0000
GEN.19502 217.800 21.0990
[ reached getOption("max.print") -- omitted 1021 rows ]

--
Report issues at https://github.com/famuvie/breedR/issues
---
You received this message because you are subscribed to the Google Groups "breedR" group.
To unsubscribe from this group and stop receiving emails from it, send an email to breedr+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/breedr/8e6c2329-f9ce-4b43-928a-25c3ff90eb84o%40googlegroups.com.

Facundo Muñoz

unread,

Jun 17, 2020, 5:52:36 AM6/17/20

to bre...@googlegroups.com

Indeed! Thanks Richard for intervening so quickly.

Hermann, you are probably fitting a univariate model for a trait which is the sum of FZ and AUSWFL.

The syntax requires:

fixed = cbind(FZ, AUSWFL) ~ GEN

ƒacu.-

To view this discussion on the web visit https://groups.google.com/d/msgid/breedr/CAM6Eiuo%3DAnZ4Zfnn99TENVraCTKoErK1%2BJ6DgMJBA0NdK8b32A%40mail.gmail.com.

--

Facundo Muñoz
Biostatistician
Département de Systèmes Biologiques
UMR ASTRE: Animal, Santé, Territoires, Risques et Écosystèmes
Campus international de Baillarguet, TA A-117 / E (E-203)
34398 Montpellier Cedex 5 (France)
+33(0) 467 593 868
@famuvie

jesus valdes hernandez

unread,

Jun 18, 2020, 4:54:45 AM6/18/20

to breedR

Hello,

Based on the previous suggestions I simulate an example of your case:

# Load sample data:

dat <- read.table(header = TRUE,

stringsAsFactors = FALSE,

text="ROW COL GEN FZ AUSWFL

742 1 1 19939 220 1

743 3 1 23003 62 3

744 4 1 44398 62 6

745 5 1 44398 62 5")

head(dat[,c("ROW","COL","GEN","FZ","AUSWFL")])

str(dat)# Here your traits are int type

library(breedR)

res <- remlf90(fixed = cbind(FZ, AUSWFL) ~ GEN,

dat = dat)# Here, Note that AI-REML has been used by default: i.e, method = 'ai'. Then, if you need, the last can be changed by EM algorithm: i.e, method = 'em'

#Exploring the results:

summary(res)

#Formula: cbind(FZ, AUSWFL) ~ 0 + Intercept + GEN

"Variance components:

Estimated variances S.E.

Residual.FZ 5386.7000 5386.7000

Residual.FZ_Residual.AUSWFL -55.9800 61.7020

Residual.AUSWFL 0.8318 0.8317

Fixed effects:

value s.e.

Intercept.FZ 2.2901e+02 111.1928

Intercept.AUSWFL -1.4235e+00 1.3817

GEN.FZ -3.8716e-03 0.0032

GEN.AUSWFL 1.5708e-04 0.0000"

#Variance components:

res[["var"]][[1]]

res[["var"]][[2]]

#Fixed effects:

res[["fixed"]]

#Apart, also you can specify: fixed and random effects like "as.factor ()", where the random is considered mere labels

#For example:

attach(dat)

dat = cbind.data.frame(as.numeric (FZ), as.numeric(AUSWFL), as.factor(GEN))

setNames(dat,c("FZ","AUSWFL","GEN"))

detach(dat)

str(dat)

#Run again the estimate breedR ...

Best regards,

Jesús

El miércoles, 17 de junio de 2020, 5:43:55 (UTC-4), Richard Whittet escribió:

Hello,

Could just be syntax..? Try replacing the "+" in the cbind with a ","

So cbind(var1, var2)

Best wishes

Richard

To unsubscribe from this group and stop receiving emails from it, send an email to bre...@googlegroups.com.

H G D

unread,

Jun 18, 2020, 5:20:06 AM6/18/20

to breedR

Dear Richard,

I realized my mistake shortly after posting. It was indeed a mixup, I should have use a "comma" instead of "+": cbind( v1, v2) rather than cbind(v1 + v2).

Thank you all very much for the help!

Best Regards,

Hermann

Am Mittwoch, 17. Juni 2020 11:43:55 UTC+2 schrieb Richard Whittet:

Hello,

Could just be syntax..? Try replacing the "+" in the cbind with a ","

So cbind(var1, var2)

Best wishes

Richard

To unsubscribe from this group and stop receiving emails from it, send an email to bre...@googlegroups.com.

Message has been deleted

H G D

unread,

Jun 20, 2020, 5:05:41 PM6/20/20

to breedR

Dear Richard and Facundo,

I have a followup question on multivariate models. I have tried fitting a multivariate model including AR1 spatial component and a genetic random effect to predict four traits.

I have tried using both, "em" and "ai" methods, both seem to converge, but "ai" is not providing model fit (AIC, LL) and only providing an intercept for the first response, "em" seems to work better here. Is there any problem with my code or onterpretation?

Best Regards, Hermann

> res <- remlf90(dat = dat, fixed = cbind(FZ, AUSWFL, AUWP, AUWPS) ~ 1,
+                random = ~ GEN, method = 'ai', # ai is faster, and em does sometimes not converge
+                spatial = list(model = 'AR', coord = dat[, c('ROW','COL')],
+                               rho = c(.8,.8)))

Using default initial variances given by default_initial_variance()
See ?breedR.getOption.

> summary(res)

Formula: cbind(FZ, AUSWFL, AUWP, AUWPS) ~ 0 + Intercept + GEN + spatial

Data: dat
AIC BIC logLik

NA NA NA

Parameters of special components:
spatial: rho: 0.8 0.8

Variance components:
Estimated variances S.E.

GEN.FZ                                       2954.6000 1.796e+02
GEN.FZ_GEN.AUSWFL                             -44.0750 3.237e+00
GEN.FZ_GEN.AUWP                              -355.5800 2.700e+01
GEN.FZ_GEN.AUWPS                               -9.6787 8.286e-01
GEN.AUSWFL                                      1.2611 7.668e-02
GEN.AUSWFL_GEN.AUWP                             9.9686 6.287e-01
GEN.AUSWFL_GEN.AUWPS                            0.3367 2.071e-02
GEN.AUWP                                       90.7380 5.517e+00
GEN.AUWP_GEN.AUWPS                              2.7607 1.728e-01
GEN.AUWPS                                       0.0940 5.718e-03
spatial.ar.FZ                                3535.2000 1.813e+02
spatial.ar.FZ_spatial.ar.AUSWFL               -43.4840 2.477e+00
spatial.ar.FZ_spatial.ar.AUWP                -544.1500 3.236e+01
spatial.ar.FZ_spatial.ar.AUWPS                -42.5070 2.386e+00
spatial.ar.AUSWFL                               0.7843 4.023e-02
spatial.ar.AUSWFL_spatial.ar.AUWP               9.4469 5.131e-01
spatial.ar.AUSWFL_spatial.ar.AUWPS              0.7002 3.716e-02
spatial.ar.AUWP                               141.3500 7.251e+00
spatial.ar.AUWP_spatial.ar.AUWPS                9.9345 5.123e-01
spatial.ar.AUWPS                                0.7128 3.657e-02
Residual.FZ                                  1012.6000 5.310e+01
Residual.FZ_Residual.AUSWFL                    21.1530 1.495e+01
Residual.FZ_Residual.AUWP                    -212.4600 1.950e+01
Residual.FZ_Residual.AUWPS                    345.6400 1.725e+01
Residual.AUSWFL                               -90.9040 1.490e+01
Residual.AUSWFL_Residual.AUWP                -420.0300 2.076e+01
Residual.AUSWFL_Residual.AUWPS                 14.4460 8.043e+00
Residual.AUWP                                -373.1400 1.842e+01
Residual.AUWP_Residual.AUWPS                    4.4864 8.923e+00
Residual.AUWPS                                  1.9950 6.365e+00

Fixed effects:
                  value   s.e.
Intercept.FZ     208.89 15.297
Intercept.AUSWFL   0.00 0.000
Intercept.AUWP     0.00 0.000
Intercept.AUWPS    0.00 0.000
> res <- remlf90(dat = dat, fixed = cbind(FZ, AUSWFL, AUWP, AUWPS) ~ 1,
+                random = ~ GEN, method = 'em', # ai is faster, and em does sometimes not converge
+                spatial = list(model = 'AR', coord = dat[, c('ROW','COL')],
+                               rho = c(.8,.8)))

Using default initial variances given by default_initial_variance()
See ?breedR.getOption.

> summary(res)

Formula: cbind(FZ, AUSWFL, AUWP, AUWPS) ~ 0 + Intercept + GEN + spatial

Data: dat
AIC BIC logLik

13676 unknown -6808

Parameters of special components:
spatial: rho: 0.8 0.8

Variance components:
$GEN
            FZ   AUSWFL     AUWP    AUWPS
FZ     2917.00 -48.7400 -539.300 -18.9100
AUSWFL -48.74   1.5060   15.920   0.6583
AUWP   -539.30 15.9200 213.500   8.6240
AUWPS   -18.91   0.6583    8.624   0.3666

$spatial
            ar.FZ ar.AUSWFL ar.AUWP ar.AUWPS
ar.FZ     2881.00 -31.7000 -380.400 -33.4500
ar.AUSWFL -31.70    0.5773    6.555   0.5299
ar.AUWP   -380.40    6.5550   96.590   7.3240
ar.AUWPS   -33.45    0.5299    7.324   0.5793

$Residual
              FZ     AUSWFL     AUWP     AUWPS
FZ     1154.0000 -10.410000 -18.1000 -0.575300
AUSWFL -10.4100   0.466600   1.2430 0.008631
AUWP    -18.1000   1.243000 10.4300 0.212700
AUWPS    -0.5753   0.008631   0.2127 0.836000

Fixed effects:
                    value    s.e.
Intercept.FZ     207.6091 13.8617
Intercept.AUSWFL   2.1997 0.2023
Intercept.AUWP    12.1687 2.5713
Intercept.AUWPS    2.0789 0.1976

Am Mittwoch, 17. Juni 2020 11:52:36 UTC+2 schrieb Facundo Muñoz:

Indeed! Thanks Richard for intervening so quickly.

Hermann, you are probably fitting a univariate model for a trait which is the sum of FZ and AUSWFL.

The syntax requires:

fixed = cbind(FZ, AUSWFL) ~ GEN

ƒacu.-

On 17/06/2020 11:43, Richard Whittet wrote:

Hello,

Could just be syntax..? Try replacing the "+" in the cbind with a ","

So cbind(var1, var2)

Best wishes

Richard

To unsubscribe from this group and stop receiving emails from it, send an email to bre...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/breedr/8e6c2329-f9ce-4b43-928a-25c3ff90eb84o%40googlegroups.com.

--
Report issues at https://github.com/famuvie/breedR/issues
---
You received this message because you are subscribed to the Google Groups "breedR" group.

To unsubscribe from this group and stop receiving emails from it, send an email to bre...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/breedr/CAM6Eiuo%3DAnZ4Zfnn99TENVraCTKoErK1%2BJ6DgMJBA0NdK8b32A%40mail.gmail.com.

jesus valdes hernandez

unread,

Jun 22, 2020, 1:43:00 PM6/22/20

to breedR

Dear Hermann,

Let me make other suggestions for your convergence problems.

A1:

-First make the combination of the bivariate analyzes and check: if there is convergence, as well as the goodness-of-fit criteria of the model (AIC BIC logLik).

-To use the values of the above estimates "Estimated variances" for multi-trait analysis (4 traits). So, these values can be incorporated via the argument var.ini as described in BreedR overview.

-So if the multi-character doesn't work for you, consider calibrating the rho argument. Here, the argument rho = c(rho_x, rho_y), for example rho = (.9,.9).

Best regards,

Jesús

H G D

unread,

Jun 22, 2020, 3:16:38 PM6/22/20

to breedR

Dear Jesus, thank you for the suggestions, I will look into them.

Best Regards,

Hermann

Facundo Muñoz

unread,

Jun 22, 2020, 3:45:17 PM6/22/20

to bre...@googlegroups.com

Dear Hermann,

It can be due to numerical issues. I can see from your EM outputs that the various traits have very different orders of magnitude: the intercept for FZ is of ~ 200 while for AUWPS if of ~2.

That factor of x100 can squares with variance components. The estimated variance for the genetic effect of FZ is of ~ 3e3 while for AUWPS if of ~0.36. Thus, 4 orders of magnitude of difference. Estimating correlations with these variation in magnitude is very unstable.

I suggest scaling the traits to a common order of magnitude (preferably on the order of tens or hundreds). Let us know if that improves the results.

ƒacu.-

To unsubscribe from this group and stop receiving emails from it, send an email to breedr+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/breedr/8f17d96d-d036-43ab-a51d-3ce1e0d1050do%40googlegroups.com.

H G D

unread,

Jun 23, 2020, 2:08:54 AM6/23/20

to breedR

Dear Facundo,

I found that removing one of the response values (AUWP or AUWPS) fixes the problem, probably since one is a transformation of the other,

Thank you very much for your assistance

Best Regards,

Hermann

To view this discussion on the web visit https://groups.google.com/d/msgid/breedr/8f17d96d-d036-43ab-a51d-3ce1e0d1050do%40googlegroups.com.

Reply all

Reply to author

Forward