Background You are interested in fecundity of an annual plant (as - - PDF document

SLIDE 1

BIOL 701 Likelihood Methods in Biology: Homework #7 Due Wednesday, May 11th

Background

You are interested in fecundity of an annual plant (as assessed by the total mass of seeds produced measured in g) as a function of fertilization regimes. You have collected data from several years in eight different fields. Within each field season, you examined plants four different treatments in each field. These treatments correspond to growing the plants with (coded as 1) and without (coded as 0) fertilization with a nitrogen source and a phosphorus source. Thus, the four treatments are: no fertilization, +N, +P, and +N,P. For each treatment×field×year combination, you recorded data for four plants. You are interested in inferring the effectiveness of different treatments in general (e.g. the expected effect of adding Nitrogen to some unspecified field), as well as learning about which fields have the highest yield. Consider a model in which each of the following effects contribute additively to the expected mass for an individual:

• 1. a year effect (centered around 0),
• 2. a field-specific expected mass without fertilization (this effect is the same for a field across all

years),

• 3. a field-specific effect of nitrogen fertilization (this effect is the same for a field across all years),
• 4. a field-specific effect of phosphorus fertilization (this effect is the same for a field across all

years)

• 5. a field-specific effect of adding both N and P (this effect is the same for a field across all

years). Of course, you should also expect some variability around this expected value (not all individuals from the same treatment, field, and year will have exactly the same mass). Note that you would like to make some general conclusions about the effect of different fertilizer treatments on some hypothetical field. So, you should structure the model so that you can learn about expected effects across fields, while accounting for the fact that their may be variability in the response of a particular field to a particular treatment. 1

SLIDE 2

• 1. Write down the likelihood for your model.

Begin Answer: i indexes the year. j indexes the field. k indexes the Nitrogen (0= no Nitrogen, 1 = Nitrogen). m indexes the Phosphorus (0= no Phosphorus, 1 = Phosphorus). q indexes the individual in that year×field×treatment. yijkmq ∼ N(αi + βj + kγj + mδj + kmρj, σ2

e)

a is the number of years b is the number of fields nijkm is the number of individuals in the specified year×field×treatment f(yijkmq|αi, βj, γj, δj, ρj, σe) = 1

• 2πσ2

e

e

−(yijkmq−αi−βj−kγj−mδj−kmρj)2 2σ2 e

(1) f(Y |α, β, γ, δ, ρ, σe) =

a

• i=1

b

• j=1

1

• k=0

1

• m=0

nijkm

• q=1

f(yijkmq|αi, βj, γj, δj, ρj, σe) (2) Using the priors introduced below (next section), we can write the likelihood in terms of the “highest level” parameters:

f(Y |σα, µβ, σβ, µγ, σγ, µδ, σδ, µρ, σρ, σe) = Z Z Z Z Z f(Y |α, β, γ, δ, ρ, σe)P(α|σα)P(β|µβ, σβ)P(γ|µγ, σγ)P(δ|µδ, σδ)P(ρ|µρ, σρ)dαdβdγdδdρ (3)

(if we take the integration symbol over a vector of latent variables to mean that we perform the integration with respect to each variable in the vector, and we calculate the definite integral from −∞ to ∞). Where P(α|σα) =

a

• i

1

• 2πσ2

α

e

−α2 i 2σ2 α

P(β|µβ, σβ) =

b

• j

1

• 2πσ2

β

e

−(βj−µβ)2 2σ2 β

P(γ|µγ, σγ) =

b

• j

1

• 2πσ2

γ

e

−(γj−µγ)2 2σ2 γ

P(δ|µδ, σδ) =

b

• j

1

• 2πσ2

δ

e

−(δj−µδ)2 2σ2 δ

P(ρ|µρ, σρ) =

a

• j

1

• 2πσ2

ρ

e

−(ρj−µρ)2 2σ2 ρ

2

SLIDE 3

2 List each parameter and the prior distribution that you have chosen to use for the parameter. Begin Answer: (there are lots of possible ways to answer - priors are up to the person analyzing the data!) a year effects (latent variables) αi ∼ Normal(0, σ2

α)

The variance of the year effects σ2

α ∼ Exponential(λ = .1)

b field effects (latent variables) βj ∼ Normal(µβ, σ2

β)

The expected value of the field effects (no fert) µβ ∼ Gamma(mean = 250, variance = 50) The variance of the field effects(no fert) σ2

β ∼ Exponential(λ = .05)

b field×nitrogen effects (latent variables) γj ∼ Normal(µγ, σ2

γ)

The expected value of the nitrogen effects µγ ∼ Normal(mean = 5, variance = 10) The variance of the nitrogen effects σ2

γ ∼ Exponential(λ = .1)

b field×phosphorus effects (latent variables) δj ∼ Normal(µδ, σ2

δ)

The expected value of the phosphorus effects µδ ∼ Normal(mean = 5, variance = 10) The variance of the phosphorus effects σ2

δ ∼ Exponential(λ = .1)

b field×nitrogen×phosphorus interaction effects (la- tent variables) ρj ∼ Normal(µρ, σ2

ρ)

The expected value of the nitrogen×phosphorus inter- action effect µρ ∼ Normal(mean = 0, variance = 10) The variance of the nitrogen×phosphorus interaction effects σ2

ρ ∼ Exponential(λ = .1)

The environmental/error standard deviation σ2

e ∼ Exponential(λ = .1)

2.5 Deriving the likelihoods - and it is most helpful to focus on the likelihood ratios for the updates: Begin Answer 3

SLIDE 4

Updating αi

When updating the a latent variables (αi), we have to consider a likelihood (the probability of the data for year i) and a prior (the probability of drawing αi from a Normal(0, σ2

α) distribution):

P(Yi | ln P(Yi | αi, β, γ, δ, ρ, σe)αi, β, γ, δ, ρ, σe) =

b

• j=1

1

• k=0

1

• m=0

nijkm

• q=1

f(yijkmq|αi, βj, γj, δj, ρj, σe) (4) Notice that this looks just like Eqn (2) except that instead of taking the product over all values

• f i, we are only doing the calculation for one value of i.

Taking the log, and substituting in for f(yijkmq|αi, βj, γj, δj, ρj, σe), we can express the likelihood and the log-likelihood ratio for an update: ln P(Yi | αi, β, γ, δ, ρ, σe) =

b

• j=1

1

• k=0

1

• m=0

nijkm

• q=1
• −1

2 ln

• 2πσ2

e

• − −(yijkmq − αi − βj − kγj − mδj − kmρj)2

2σ2

e

• =

−ni 2 ln

• 2πσ2

e

• −

1 2σ2

e b

• j=1

1

• k=0

1

• m=0

nijkm

• q=1

(yijkmq − αi − βj − kγj − mδj − kmρj)2 ln LR(αi → α⋆

i )

= ln P(Yi | α⋆

i , β, γ, δ, ρ, σe) − ln P(Yi | αi, β, γ, δ, ρ, σe)

= − 1 2σ2

e b

• j=1

1

• k=0

1

• m=0

nijkm

• q=1

(yijkmq − α⋆

i − βj − kγj − mδj − kmρj)2

. . . + 1 2σ2

e b

• j=1

1

• k=0

1

• m=0

nijkm

• q=1

(yijkmq − αi − βj − kγj − mδj − kmρj)2 = − 1 2σ2

e b

• j=1

1

• k=0

1

• m=0

nijkm

• q=1

(Zα(i, j, k, m, q) − α⋆

i )2 +

1 2σ2

e b

• j=1

1

• k=0

1

• m=0

nijkm

• q=1

(Zα(i, j, k, m, q) − αi)2 Zα(i, j, k, m, q) = yijkmq − βj − kγj − mδj − kmρj ln LR(αi → α⋆

i )

= 1 2σ2

e

 

b

• j=1

1

• k=0

1

• m=0

nijkm

• q=1
• (Zα(i, j, k, m, q) − αi)2 − (Zα(i, j, k, m, q) − α⋆

i )2

  = 1 2σ2

e

 

b

• j=1

1

• k=0

1

• m=0

nijkm

• q=1
• Zα(i, j, k, m, q)(α⋆

i − αi) + α2 i − (α⋆ i )2

  = ni(α2

i − (α⋆ i )2)

2σ2

e

+ (α⋆

i − αi)

2σ2

e

 

b

• j=1

1

• k=0

1

• m=0

nijkm

• q=1

Zα(i, j, k, m, q)   = ni(α2

i − (α⋆ i )2)

2σ2

e

+ (α⋆

i − αi)

2σ2

e

Wα(i) (5) Wα(i) =

b

• j=1

1

• k=0

1

• m=0

nijkm

• q=1

Zα(i, j, k, m, q) ni =

b

• j=1

1

• k=0

1

• m=0

nijkm

• q=1

1 4

SLIDE 5

This is just one example of a log-likelihood ratio used in the update of a latent variable. Note that the relevant sums can be calculated using the: l a t e n t v a r i a b l e . value = 0.0 sum ot her ef fec ts = d a t a f o r e f f e c t . c a l c u l u a t e s u m o v e r e f f e c t s () l a t e n t v a r i a b l e . value = current value trick shown in the template1 to get the sum of effects shown that is one part of Wα(i). The data set’s by year attribute is such a list - each element represents one year’s worth of data (the SameFieldTreatmentYearGroup for the year are stored in a list that is an instance of SharedEffectGroupList; in essence this is just a list with a few helper functions to let you Note that we can use the same trick with a call to calc sum effects on a SameFieldTreatmentYearGroup

• bject if we wanted to calculate the sum of effects in the process of calculating something like

nijkm

q=1

Zα(i, j, k, m, q) for some particular group of individuals with that experienced the same year, field and treatment.

Updating βj, γj, δj, ρj

The updates for the βj latent variables will be just like update of αi except the summations will need to be over a

i instead of b j, and you will need a Zβ, Wβ, and nj sums that remove the

effects of βj. Plus you’ll have to be walking over the relevant “slice” of the data when you do the update. Similar modifications of the algorithm are required for updating γj, δj, and ρj.

Updating σe

When updating σe, the relevant summary statistic that we can calculate by summing over individ- uals is the sum of squared errors: SE(i, j, k, m, q) = (yijkmq − αiβj − kγj − mδj − kmρj)2 SSE =

a

• i=1

b

• j=1

1

• k=0

1

• m=0

nijkm

• q=1

SE(i, j, k, m, q) ln LR(σe → σ⋆

e)

= n 2 ln

• 2π(σe)2

+ SSE 2(σe)2 − n 2 ln

• 2π(σ⋆

e)2

− SSE 2(σ⋆

e)2

= n 2 ln σ2

e

(σ⋆

e)2

• + SSE

2 1 σ2

e

− 1 (σ⋆

e)2

• (6)

1sorry for the “calculuate” rather than “calculate” typo in the function name

5

SLIDE 6

Updating µβ

In a hierarchical model with latent variables (like ours), updating µβ does not change the actual values of β. Thus changing µβ does not affect Equation (2), but it does change Eqn (3). In particular, the only factor that changes in Eqn (3) is P(β|µβ, σβ). So the likelihood ratio of this update is based on the probability statement that acts as a prior on β. Of course there is still a prior on µβ, so the move have prior ratio and a likelihood ratio: P(β|µβ, σβ) =

b

• j

1

• 2πσ2

β

e

−(βj−µβ)2 2σ2 β

ln P(β|µβ, σβ) =

b

• j
• −1

2 ln 2πσ2

β − (βj − µβ)2

2σ2

β

• ln LR(µβ → µ⋆

β)

= 1 2σβ  

b

• j
• (βj − µβ)2 − (βj − µ⋆

β)2

  = 1 2σβ  

b

• j
• 2βj(µ⋆

β − µβ) + µ2 β − (µ⋆ β)2

  = b

• µ2

β − (µ⋆ β)2

+ 2(µ⋆

β − µβ)Zµβ

2σβ (7) Zµβ =

b

• j

βj Updates to the other parameters that place priors on the latent variables will use similar logic. 6