Multivariate estimation of genetic parameters Quo vadis? Karin - - PowerPoint PPT Presentation

Multivariate estimation of genetic parameters – Quo vadis?

Karin Meyer

Animal Genetics and Breeding Unit, University of New England, Armidale

AAABG 2011

REML - quo vadis?

“Quo vadis?” Is a Latin phrase meaning: "Where are you going?" or "Whither goest thou?"

Wikipedia

Disclaimer:

Any similarity to previous AAABG title(s) is pure coincidence!

Statistical methods I MIXED MODELS IN ANIMAL BREEDING: WHERE TO NOW? A.R. Gilmour Cargo Vale, CARGO, NSW 2800, formerly Orange Agricultural Institute, NSW Department of Primary Industries SUMMARY Over the past 60 years, mixed models have underpinned huge gains in plant and animal production through genetic improvement. Charles Henderson (1912-1989) established mixed models for estimating breeding values (BLUP) using the popularly called Henderson's Mixed Model and provided early methods (Henderson's Methods I, II and III) for estimating variance

• parameters. Robin Thompson then published the widely acclaimed REML method for variance
• K. M. |

| AAABG 2011 2 / 24

REML - quo vadis?

Outline

Introduction Penalized REML Improved estimator Penalties Tuning factor Simulation results Application Finale

• K. M. |

| AAABG 2011 3 / 24

www.wordle.net

REML - quo vadis? | Introduction

Motivation

Multivariate estimation: more than a few traits

increasingly common −→ more traits of interest −→ technical capability: Hard/Soft-ware

• no. of parameters ↑ quadratically

→ q traits −→ q(q + 1)/2 covariances SAMPLING VARIATION ↑ with no. of parameters/traits

• K. M. |

| AAABG 2011 4 / 24

REML - quo vadis? | Introduction

Motivation

Multivariate estimation: more than a few traits

increasingly common −→ more traits of interest −→ technical capability: Hard/Soft-ware

• no. of parameters ↑ quadratically

→ q traits −→ q(q + 1)/2 covariances SAMPLING VARIATION ↑ with no. of parameters/traits

Alleviate S.V.

Data: large −→ gigantic samples Model: parsimony −→ less parameters than covar.s

e.g. RR, Reduced rank, structured

Estimation −→ use additional information

Bayesian: Prior REML: Impose penalty P on likelihood

• K. M. |

| AAABG 2011 4 / 24

REML - quo vadis? | Introduction

Key message

Yes, we can: Reduce effects of sampling variation in multivariate analyses −→ ‘better’ estimates Achieve this through ‘simple’ extension of standard REML −→ penalty

• K. M. |

| AAABG 2011 5 / 24

Penalized REML

REML - quo vadis? | Penalized REML | Improved estimator

What is a ‘better’ estimator?

Quality Some measure of deviation from true value Loss e.g. MSE = Sampling variance + Bias2 Covariance matrix

Entropy loss

(James & Stein 1961)

L1

• Σ, ˆ

Σ

• = tr
• Σ−1 ˆ

Σ

• − log
• Σ−1 ˆ

Σ

• − q

Quadratic loss L2

• Σ, ˆ

Σ

• = tr
• Σ−1 ˆ

Σ − I

2

• K. M. |

| AAABG 2011 6 / 24

REML - quo vadis? | Penalized REML | Improved estimator

What is a ‘better’ estimator?

Quality Some measure of deviation from true value Loss e.g. MSE = Sampling variance + Bias2 Covariance matrix

Entropy loss

(James & Stein 1961)

L1

• Σ, ˆ

Σ

• = tr
• Σ−1 ˆ

Σ

• − log
• Σ−1 ˆ

Σ

• − q

Quadratic loss L2

• Σ, ˆ

Σ

• = tr
• Σ−1 ˆ

Σ − I

2

Improved Modify estimator to reduce loss Sampling variance Bias Loss

• K. M. |

| AAABG 2011 6 / 24

REML - quo vadis? | Penalized REML | Improved estimator

Penalized maximum likelihood

‘Standard’ (RE)ML −→ Maximize logL(θ) w.r.t. θ Penalized logLP(θ) = logL(θ) − 1

2 ψ

P

Tuning factor Penalty −→ f(θ)

Crucial questions:

What kind of penalty? −→ P should reduce S.V. How to determine the tuning factor? −→ relative emphasis: data ↔ penalty

• K. M. |

| AAABG 2011 7 / 24

REML - quo vadis? | Penalized REML | Improved estimator

Toy example

Penalty changes shape

• f likelihood!

PHS design s=100, n=10 q = 1, h2 = 0.6 Penalty

P = (h2 − 0.3)2

Tuning factor

Ψ = 0, . . . , 100

• K. M. |

| AAABG 2011 8 / 24

Heritability log L 0.2 0.4 0.6 0.8 −15 −10 −5

h2 = 0.6

ψ = 0

REML - quo vadis? | Penalized REML | Penalties

Choosing a penalty: rationale

Estimate genetic parameters −→ partition total variance −→ sampling correlations ˆ σG ij & ˆ σE ij Observation: ˆ

ΣP estimated more accurately than ˆ ΣG

Idea: ‘Borrow strength’ from ˆ

ΣP

• K. M. |

| AAABG 2011 9 / 24

REML - quo vadis? | Penalized REML | Penalties

Choosing a penalty: rationale

Estimate genetic parameters −→ partition total variance −→ sampling correlations ˆ σG ij & ˆ σE ij Observation: ˆ

ΣP estimated more accurately than ˆ ΣG

Idea: ‘Borrow strength’ from ˆ

ΣP

How:

1

Modify eigenvalues of ˆ

Σ

−1

P

ˆ ΣG

−→ λi

“canonical” eigenvalues [0, 1]

2

Assume ˆ

ΣG ∼ IW distribution; scale ˆ ΣP

−→ general principle: P ∝ minus log of prior density → empirical Bayes approach

• K. M. |

| AAABG 2011 9 / 24

REML - quo vadis? | Penalized REML | Penalties

Penalty on eigenvalues

Eig.values of ˆ

Σ overdispersed

(Lawley 1956)

Overdispersion −→ S.V. Modify eigenvalues to ↓ S.V. Regress ˆ

λi towards mean ¯ λ

1 2 3 4 5 6 7 8 0.8 1 1.2

• Distribution of sample eigenvalues

8 traits yi ∼ N(0, I) n=500 S =

i yiy′ i /(n − 1)

1000 replicates

• K. M. |

| AAABG 2011 10 / 24

REML - quo vadis? | Penalized REML | Penalties

Penalty on eigenvalues

Eig.values of ˆ

Σ overdispersed

(Lawley 1956)

Overdispersion −→ S.V. Modify eigenvalues to ↓ S.V. Regress ˆ

λi towards mean ¯ λ

1 2 3 4 5 6 7 8 0.8 1 1.2

• Distribution of sample eigenvalues

8 traits yi ∼ N(0, I) n=500 S =

i yiy′ i /(n − 1)

1000 replicates

• K. M. |

| AAABG 2011 10 / 24

λ∗= ¯ λ+0.4(λ − ¯ λ)

REML - quo vadis? | Penalized REML | Penalties

Penalty on eigenvalues

Eig.values of ˆ

Σ overdispersed

(Lawley 1956)

Overdispersion −→ S.V. Modify eigenvalues to ↓ S.V. Regress ˆ

λi towards mean ¯ λ

‘Bending’ (Hayes & Hill 1981)

ˆ Σ

∗

G = β ˆ

ΣG + (1 − β) ¯ λ ˆ ΣP

Equivalent for (RE)ML

Pλ ∝

• i( ˆ

λi − ¯ λ)2 P ℓ

λ ∝

• i(log( ˆ

λi) − log(λ))2

1 2 3 4 5 6 7 8 0.8 1 1.2

• Distribution of sample eigenvalues

8 traits yi ∼ N(0, I) n=500 S =

i yiy′ i /(n − 1)

1000 replicates

• K. M. |

| AAABG 2011 10 / 24

λ∗= ¯ λ+0.4(λ − ¯ λ)

REML - quo vadis? | Penalized REML | Penalties

Penalty on eigenvalues - cont.

Quadratic penalty implies

λi ∼ N( ¯ λ, σ2)

Alternative: allow for non-common mean −→ Assume λi ∼ Order statistics on unit interval → Beta

Pβ ∝

• i(i − 1) log(λi) + (q − i) log(1 − λi)

0.5 1 1 2 3

PDF: Order statistics for q=5

• K. M. |

| AAABG 2011 11 / 24

REML - quo vadis? | Penalized REML | Penalties

Penalty on matrix divergence

Assume ˆ

ΣG has Inverse Wishart prior

Substitute ˆ

Σ

P (unpenalized) for scale parameter

PΣ ∝ C log | ˆ ΣG| + tr(ˆ Σ

−1

G ˆ

Σ

P)

Shrink ˆ

ΣG towards ˆ ΣP By analogy: genetic correlation matrix ˆ RG

Shrink ˆ RG towards ˆ R

P

PR ∝ C log |ˆ

RG| + tr(ˆ R

−1

G ˆ

R

P)

Advantages

Easy to implement for standard REML parameterization Extension to penalize multiple matrices straightforward

• K. M. |

| AAABG 2011 12 / 24

C = (ψ+q+1)/ψ

≈ 1

REML - quo vadis? | Penalized REML | Tuning factor

How to choose that tuning factor?

Pick ψ a priori −→ ‘degree of belief’

e.g. based on sample size

Estimate from the data

Cross-validation

Split data into ‘training’ & ‘validation’ sets Obtain ˆ

ΣG and ˆ ΣE for range of ψ → training

Evaluate logL(ˆ

ΣG, ˆ ΣE) for all ψ → validation

Pick ψ which maximizes logL(ˆ

ΣG, ˆ ΣE) K− fold cross-validation

K data subsets, K analyses

• K. M. |

| AAABG 2011 13 / 24

REML - quo vadis? | Penalized REML | Tuning factor

How to choose that tuning factor?

Pick ψ a priori −→ ‘degree of belief’

e.g. based on sample size

Estimate from the data

Cross-validation

Split data into ‘training’ & ‘validation’ sets Obtain ˆ

ΣG and ˆ ΣE for range of ψ → training

Evaluate logL(ˆ

ΣG, ˆ ΣE) for all ψ → validation

Pick ψ which maximizes logL(ˆ

ΣG, ˆ ΣE) K− fold cross-validation

K data subsets, K analyses

Limit change in logL

Pick largest ψ so that ∆logL ≈ chosen limit Choice of limit? e.g. mild penalty: |∆logL| ≈ 1

2χ2

α,df=1

• K. M. |

| AAABG 2011 13 / 24

Simulation

REML - quo vadis? | Simulation results

Does it work?

Check by simulation∗

q=5 traits h2

1 ≥ . . . ≥ h2 5, MVN

• Pat. half-sib design

s=100, n=10 60 sets of pop. values

• Coeff. Var(λi)

0 to 170%

1000 replicates

Percentage Reduction In Average Loss PRIAL = 100

      

1 −

¯

L1(Σ, ˆ

Σ

ˆ ψ)

¯

L1(Σ, ˆ

Σ0)

      

¯

L1 mean entropy loss

unpenalized penalized

∗more in papers 2 & 3

• K. M. |

| AAABG 2011 14 / 24

L1(Σ, ˆ

Σ) = tr(Σ−1 ˆ Σ) − log |Σ−1 ˆ Σ| − q

REML - quo vadis? | Simulation results

PRIAL in ˆ

ΣG: optimal

PRIAL Σ ^

G

20 40 60 80 100

• 36

71 68 71

Pλ Pλ

l

Pβ PΣ

Use population values to determine ψ −→ ‘best possible’

Pλ worse

−→ over-shrink if λi ≈ ¯ λ Other P similar −→ PΣ least variable

Pλ shrink ˆ λi towards mean P ℓ

λ shrink log ˆ

λi towards mean Pβ assume ˆ λi ∼ Order statistics PΣ shrink ˆ ΣG towards ˆ Σ

P

• K. M. |

| AAABG 2011 15 / 24

REML - quo vadis? | Simulation results

PRIAL in ˆ

ΣG: estimate ψ

PRIAL Σ ^

G

20 40 60 80 100

• 23

56 61 55

Pλ Pλ

l

Pβ PΣ

PRIAL Σ ^

G

20 40 60 80 100

• 41

68 69 64

Pλ Pλ

l

Pβ PΣ

3−fold cross-validation

PRIAL ↓ by 10-15%

ψ too large

|∆logL| ≤ 1

2χ2 5%,1 = 1.92

ψ too small −→ mild P P ℓ

λ & PΣ ↓ 3-6%

• K. M. |

| AAABG 2011 16 / 24

REML - quo vadis? | Simulation results

Bias: Canonical eigenvalues

Relative bias (in %)†

ˆ λi

NoP

Pλ P ℓ

λ

Pβ PΣ

1 10

• 11
• 4
• 8

8

†using 3-fold cross-validation to estimate ψ

• K. M. |

| AAABG 2011 17 / 24

REML - quo vadis? | Simulation results

Bias: Canonical eigenvalues

Relative bias (in %)†

ˆ λi

NoP

Pλ P ℓ

λ

Pβ PΣ

1 10

• 11
• 4
• 8

8 2 27 16 17 21 26 3 17 22 27 26 25 4

• 19

11 53 28 42 5

• 79
• 26

107 35 87

†using 3-fold cross-validation to estimate ψ

• K. M. |

| AAABG 2011 17 / 24

REML - quo vadis? | Simulation results

Bias: Heritabilities

ˆ

h2

i

NoP

Pλ P ℓ

λ

Pβ PΣ

1

• 1
• 14
• 7
• 11

1 2 4

• 5

4 11 3 5

• 1

11 6 16 4 7 7 23 15 27 5 12 20 45 31 46 No P: Bias due only to constraints on param.s −→ REML is biased! Pattern similar to bias ( ˆ

λi) P on ˆ λi

high ↓ low ↑

• esp. log ˆ

λi

P on ˆ ΣG

high ↑ ↓ low ↑

• K. M. |

| AAABG 2011 18 / 24

Application

REML - quo vadis? | Application

Application: Data & Analysis

Carcass traits, Hereford cattle, CRC-I

(Reverter et al., 2000)

WT 1030

                    

30% all 6 traits 54% five traits EMA 864 IMF 992 RBY 370 P8 999 RIB 1014 Simple animal model, 2817 animals in pedigree

42 covariance components to be estimated Penalties: P ℓ

λ shrink canonical eigenvalues on log scale towards mean

PΣ shrink ˆ ΣG towards ˆ ΣP

Estimate ψ: 3−fold cross-validation

• K. M. |

| AAABG 2011 19 / 24

REML - quo vadis? | Application

Heritability estimates with penalty

P ℓ

λ

ψ = 9.5 ∆logL(θ) = −5.1

(ψ∆5% = 2.9)

PΣ

17.0

−3.2

(

9.5)

Heritability WT EMA IMF RBY P8 RIB 0.4 0.8

• Penalty

None Can.eigenvalues IW cov. matrix

• K. M. |

| AAABG 2011 20 / 24

REML - quo vadis? | Application

Genetic correlations . . .

Genetic correlation

• Penalty

None Can.eigenvalues IW cov. matrix

−1 −0.5 0.5 1 P8,RIB EMA,RIB WT,RBY IMF ,RIB EMA,P8 EMA,IMF IMF ,P8 EMA,RBY WT,EMA WT,IMF RBY,RIB IMF ,RBY WT,RIB RBY,P8 WT,P8

• K. M. |

| AAABG 2011 21 / 24

Finale

REML - quo vadis? | Finale

What’s next?

Investigate alternative penalties

e.g. P based on Beta prior for λi

Extend methodology:

shrinkage towards structured matrix random regression models differential degree of penalization for traits with many & few records

Apply same principles: combining estimates from separate analyses? Fine tune implementation in WOMBAT Convince people to use it!

• K. M. |

| AAABG 2011 22 / 24

REML - quo vadis? | Finale

Conclusions

Can improve multivariate estimates of genetic param.s Straightforward extension to standard REML −→ penalize logL(θ)

→ ‘borrow’ strength from ˆ ΣP → exploit prior knowledge

• K. M. |

| AAABG 2011 23 / 24

REML - quo vadis? | Finale

Conclusions

Can improve multivariate estimates of genetic param.s Straightforward extension to standard REML −→ penalize logL(θ)

→ ‘borrow’ strength from ˆ ΣP → exploit prior knowledge

Most useful for small/moderate sample sizes; ≥ 4 traits

current favourite: PR, use ∆logL

Disadvantages −→ some increase in bias −→ higher computational requirements

• K. M. |

| AAABG 2011 23 / 24

REML - quo vadis? | Finale

Conclusions

Can improve multivariate estimates of genetic param.s Straightforward extension to standard REML −→ penalize logL(θ)

→ ‘borrow’ strength from ˆ ΣP → exploit prior knowledge

Most useful for small/moderate sample sizes; ≥ 4 traits

current favourite: PR, use ∆logL

Disadvantages −→ some increase in bias −→ higher computational requirements Obvious with hindsight No excuses: make best possible use of precious data

• K. M. |

| AAABG 2011 23 / 24

REML - quo vadis?

Where are we going?

• K. M. |

| AAABG 2011 24 / 24