La rgression sous WinBugs: une vieille mthode revisite partir d'un - PowerPoint PPT Presentation

JOURNEES Club WinBUGS 7 décembre 2006 La régression sous WinBugs: une vieille méthode revisitée à partir d'un exemple environnemental. ������������������� ����������� ������������������ Eric PARENT, Etienne RIVOT MOdélisation, Risque, Statististique, Environnement de l�’UMR MIA 518 INRA/ENGREF/INAPG 1 (Math. Info. App. 518)

Basics for regression • Population vs sample • True parameters vs estimates • Hypotheses

Représentation graphique des notions suivantes : population (bleues), modèle (droite bleue), échantillon (croix rouges), modèle ajusté (droite rouge)

Représentation graphique de l'échantillon (croix rouges) et du modèle ajusté (droite rouge). La droite rouge obtenue est conditionnelle à l'échantillon en main (4 échantillons � 4 droites différentes). Les coefficients (pente, ordonnée à l'origine et r 2 ) sont donc incertains

Hypotheses • Residuals are independant • Residuals are from the same pdf • Residuals are normal • model structure

= − � en fonction de Distribution des résidus e y y i i i = + y a bx � . i i

− µ X µ σ = known Z ( , ) ; σ Normal PDF

Bayes theory for the linear model

Remember Bayesian statistics ? ����� ���������� [ ] [ ] [ ] θ × θ y θ = ��������� y [ ] y ������������������������������������ DATA (y) Prior information Posterior information [ ] [ ] θ θ y

Bayesian analysis of the linear model = α + β + ε y x ε σ N 2 ~ ( 0 , ) where t t t t n [ ]   { ( )( ) }   1 1 α β σ = × − − θ − θ y X   Y X Y X   2 , , , exp ’ σ π σ     2 2 2 − a 1 b       1 1 1 ∑ − 1 θ ∝ − − θ − θ θ − θ H     , exp exp ( )’ ( )   σ σ σ 0 0       2 2 2 0 2 y x 1 θ θ 0 Σ N 1 1 ~ ( , ) 1/ σ 2 ~Gamma(a,b) 0 ... 1 ... = = θ = α β Y y X x ; ... ; t t ... ... ... y x 1 n n

Bayesian analysis of the linear model = X θ + ε ε σ Y N 2 I ~ ( 0 , ) where ( ) − θ = ˆ X X 1 X Y ’ ’ Least square estimate n { ( )( ) } [ ]     1 1 α β σ = × − − θ − θ ˆ + θ ˆ − θ − θ ˆ + θ ˆ y X   Y X Y X   2 , , , exp ( ' ( ) σ π σ  2    2 2 Pythagore�’ theorem { ( ) ( ) ( ) ( )( ) } [ ]   1 α β σ ∝ − θ − θ ˆ θ − θ ˆ + − θ ˆ − θ y X X X Y X Y X ˆ   2 , , , exp ' ' ' σ   2 2 θ only! − a 1 b       1 1 1 ∑ − 1 θ ∝ − θ − θ θ − θ − H     ( ) exp , exp ( )’   σ σ σ 0 0       2 2 2 0 2

Bayesian analysis of the linear model = X θ + ε ε σ Y N 2 I ~ ( 0 , ) where ( ) − θ = ˆ X X 1 X Y ’ ’ Least square estimate ( ) ( )       X X  1 1 ' ∑ − 1 θ ∝ −  θ − θ θ − θ + θ − θ ˆ θ − θ ˆ  Y X   , , exp ( )' ( ) '    σ σ 0 0  2   2   0  2 θ σ θ Σ N 2 If σ were known, θ would be normal ~ ( , ) 1 1 X X   ' ∑ ∑ − − 1 1 = +   Addition of the precision matrices σ   2 1 0  X X  ' ∑ ∑ − − 1 θ = 1 θ + θ ˆ   Linear combination between prior σ 1 0   2 1 0 mean and least square estimates

Bayesian analysis of the linear model = X θ + ε ε σ Y N 2 I ~ ( 0 , ) where Gamma (a,b) n − a 1     b   ( ) 1 1 1 1 2 θ ∝ − − − θ − θ Y X Y X Y X     , , exp exp ( )’ ( )   σ σ σ πσ σ       2 2 2 2 2 2 2 If θ were known, σ 2 would be inverse Gamma�(a’,b’) n = + a a 1 ’ θ gamma a b ) ~ ( ’ , ’ 2 σ 2 − θ − θ Y X Y X ( )’ ( ) = + b b ’ 2

Stock Recruitment issues • Fishing should not endanger population conservation • How many spawners are necessary to guarantee salmon population conservation ? • How does the number of spawners influence on the recruitment ? • Models of stock/recruitment relationship wanted

Margaree River

Margaree River SALMON: The Salmon fishing season begins June 1 and ends October 31 ( fly only ).Resident fees are set at $36.06. Non- resident anglers must pay $133.19 for a seasonal license or $54.26 for a 7 day license. Bag limit: 2 per day and 8 per season ( only grilse up to 63 cm 24.8 inches )

Cohorte Géniteurs Recrues 1947 1685 4852 1948 3358 7204 1949 1839 5716 1950 1744 4000 Margaree River 1951 2093 2440 1952 969 2833 1956 486 2616 1957 822 4534 1961 344 3620 1962 1306 3850 1963 887 3538 1964 1053 2515 1965 993 3694 1966 727 1393 1967 1009 2083 1968 828 2378 1969 488 3394 1970 901 2702 1971 351 2630 1972 373 3261 MARGAREE River 1973 393 3131 1974 436 1066 1975 293 2813 8000 1976 366 1819 1977 538 2909 7000 1978 699 3292 6000 1979 363 1868 1980 681 1462 5000 1981 618 3616 Recrues 1982 760 4015 4000 1983 657 1688 3000 1984 381 2289 1985 1378 5156 2000 1986 3461 3484 1987 3899 6375 1000 1988 1545 3358 0 1989 2164 2900 1990 5022 2365 0 1000 2000 3000 4000 5000 6000 Géniteurs

A look at the data of SR Margaree river 8000 Recrutement (nb retours adultes) 7000 6000 5000 4000 3000 2000 1000 0 0 1000 2000 3000 4000 5000 6000 Stock (nb reproducteurs)

Pick your own model!

Searching for a model • R as a function of S : a smooth idealized phenomenon R= µ R (S), – when S ↑ , R ↑ (at least at the beginning) – ⇒ R est fonction de S • Indeed R is uncertain for many reasons (S is the only explanation available) ⇒ stochastic model • Math. translation : R ~ pdf with caracteristics ( µ R (S), σ ) position parameter µ R (S), dispersion parameter σ – dispersion parameter σ is not a fonction of S – µ R (S) ? Which pdf ? •

Searching for µ R (S) µ R (S) ? • – Ecological knowledge : competition for territory R ↑ then S ↑ , then R stops or decreases – ⇒ Feedback control by density • 2 classical functions : – Beverton-Holt µ R (S) = aS/(1 + bS) reproduction ratio µ R (S)/S = a/(1 + bS) – Ricker µ R (S) = aSe -bS reproduction rate µ R (S)/S = ae -bS Log( µ R (S)/S)=Log(a)-bS –

Fonctions Ricker et BH 60 3,5 3 50 2,5 40 muR(S)/S muR(S) 2 BH BH 30 Ricker Ricker 1,5 20 1 10 0,5 0 0 0 0 50 50 100 100 150 150 200 200 Beverton-Holt S S µ R (S)/S = a/(1 + bS) Ricker µ R (S)/S = ae -bS

Ricker = a hidden linear model • Summing up R ~ logN( µ R (S), σ ) – µ R (S) = aSe -bS (Ricker function) – – deterministic explanatory structure + random effect • ε − = µ ε σ µ = bS R S e N S aSe 2 ( ) ~ ( 0 , ) ( ) R R = + − + ε R a S bS log( ) log( ) log( )  R  = = α = β = − θ = y x S a b a b   log log( ) (log( ), ) S   = α + β + ε ε σ y x N 2 ~ ( 0 , )

DAG of the Ricker SR model: classical formulation a ~ unif ß ~ unif Ricker function : s t R(s t ) = s t exp( a - ß s t ) µ R t = log(R( s t )) σ ~ uninf r t, ~ LogN( µ R t , σ ) log( σ 2 ) ~ unif t = 1 to T observations

DAG of the Ricker SR model: classical formulation Slope = exp(a) S max = 1/ß a ~ unif ß ~ unif Ricker function : s t R( s t ) = s t exp( a - ß s t ) µ R t = log(R( s t )) σ ~ uninf r t, ~ LogN( µ R t , σ ) log( σ 2 ) ~ unif t = 1 to T observations

15.0 10.0 tau 0.05.0 9999 20000 30000 iteration Posterior inference…. alpha sample: 30001 4.0 model fit: r.star P(alpha) 2.0 1.5E+4 0.0 1.25 1.5 1.75 2.0 2.25 2.5 alpha beta sample: 30001 1.0E+4 6.0E+3 r.star P(beta) 0.0 2.0E-4 4.0E-4 6.0E-4 beta 5.0E+3 tau sample: 30001 0.4 P(tau) 0.2 0.0 0.0 0.0 5.0 10.0 tau 0.0 2.0E+3 4.0E+3 6.0E+3

Reference points of the Ricker function 12 10 8 MSY oeufs 6 recrutés/m² 4 S MSY 2 0 0 5 10 15 oeufs déposés/m² production remplacement excédent

CHANGE of PARAMETERS Schnute & Kronlund formulation C opt ~ ? S opt ~ ? a ß = log((C opt + S opt )/S opt ) + (C opt /(S opt + C opt )) = C opt /(S opt (S opt + C opt )) Ricker function : s t R( s t ) = s t exp( a - ß s t ) µ R t = log(R( s t )) σ ~ uninf r t, ~ LogN( µ R t , σ ) log( σ 2 ) ~ unif t = 1 to T observations

Which priors to assign to C opt and S opt ? - Option 0 : uniformative priors to and - Option 1: assign a prior equivalent to the reference prior assigned to and - Option 2: propose a little informative prior(s) based on some rationale