Monitoring and data filtering II. Dynamic Linear Models Advanced - - PDF document

1

Dias 1

Monitoring and data filtering

• II. Dynamic Linear Models

Advanced Herd Management Cécile Cornou, IPH

Dias 2

Outline

Introduction to the DLM (West and Harrison, chapter 2)

• Updating equations: Kalman Filter
• Discount factor as an aid to choose W
• Incorporate external information: Intervention

General form of the DLM Examples Concluding remarks

2

Dias 3

Introduction to the DLM

A Simple DLM Time series Yt = (y1, ... , yn) Observation equation: yt = µt + vt, vt ∼ N(0, Vt)

• Like before: vt = es + eo
• The symbol µt is the underlying true value at time t.

System equation: µt = µt-1 + wt, wt ∼ N(0, Wt) The true value is not any longer assumed to be constant. A fair assumption in animal production Basically, we wish to detect ”large” changes in µt

Dias 4

Introduction to the DLM

A DLM with a trend Time series Yt = (y1, ... , yn) Observation equation: yt = µt + vt, vt ∼ N(0, Vt) System equation: µt = µt-1 + βt-1 + w1t, w1t ∼ N(0, W1t) βt = βt-1 + w2t, w2t ∼ N(0, W2t)

3

Dias 5

Introduction to the DLM

Updating equations: Kalman Filter (a) Posterior for µt-1 : (µt-1 | Dt-1) ∼ N(mt-1, Ct-1) (b) Prior for µt : (µt | Dt-1) ∼ N(mt-1, Rt) where Rt = Ct-1 + Wt (c) 1-step forecast: (Yt | Dt-1) ∼ N(ft, Qt) where ft = mt-1 and Qt = Rt + Vt (d) Posterior for µt : (µt | Dt) ∼ N(mt, Ct) with mt = mt-1 + At.et and Ct = At.Vt where At = Rt / Qt and et = Yt - ft Initial Information: (µ0 | D0) ∼ N(m0, C0)

Dias 6

Introduction to the DLM

Discount factor as an aid to choosing Wt To run the model (assume with constant parameters) we need: m0 , C0 , V, W Discount factor can be used if W is unknown we know that W is a fixed proportion of C (West & Harrison 2.4.2) Rt = Ct-1 + W → Rt = Ct-1 / δ Typically 0.8 < δ < 1

4

Dias 7

Introduction to the DLM

Incorporate external information: intervention Types of external information: 1.

• 1. Known effect, experienced before (ex: change in breed for

which we know the different performances) → We want the model to adapt to the new known conditions

• 2. Unknown effect (ex: wave of heat, introduction of new animals

in a group) → We want the model to adapt to the new unknown conditions

• 3. Unknown effect we want to measure (ex: change of feed

composition, new veterinary treatments) → We want to measure the effect of a volontary change

Dias 8

Introduction to the DLM

Intervention - 1. Known effect We want the model to adapt to the new known conditions Ex: Kurrit example from West and Harrison (2.3.2) Estimated mean after change: 286 (vs. 143) Expected change = 143 (286 – 143) Uncertainty: from 80 (pessimistic) to 200 (optimistic) σ = 30 ( = (200 – 80) / 4) : 4 st.dev (95% interval) Variance associated (Uncertainty) = 302 = 900 (ω10 | D9, S9) ∼ N (143, 900) Revised one-step ahead forecast: (µt | Dt-1) ∼ N(mt-1, Rt) (µ10 | D9, S9) ∼ N (286, 920) m9 = 286 and R10 = C9 + W10 = 20 + 900 = 920

5

Dias 9

Introduction to the DLM

Intervention - 2. Unknown effect (1/2) We want the model to adapt to the new unknown conditions Ex 1: wave of heat We can not adjust because we do not know the exact effect Ex 2: introduction of new animals in a group Consider a method aimed to detect oestrus by monitoring animal behaviour Incoming animals may modify the behaviour of the group Here, intervention aimed to increase model adaptation to new behaviour so to avoid ”alarms” due to a known event

Dias 10

Introduction to the DLM

Intervention - 2. Unknown effect (2/2) In practice we can temporarily reduce the value of the discount factor so the evolution variance increases We put more weight on the new observations and ”forget about the past”

See also eating rank

6

Dias 11

Introduction to the DLM

Intervention - 3. Unknown effect We want to measure the effect of a volontary change Ex: A new feed is used and we want to estimate the associated change in daily gain We know that the new feed is used from time τ (0 < τ < n) yt = µt + λt It + vt, vt ∼ N(0, Vt) µt = µt-1 + wt, wt ∼ N(0, Wt) With: It : intervention effect that we want to measure λt = 0 when t < τ λt = 1 when t > τ

Dias 12

The general DLM

Generalisation from the 1.order pol. Model

• Simple, most widely used DLM
• Matrix notation allows to present the DLM in a general form and to treat

more complex cases Three examples of application

• Monitoring activity level
• Monitoring activity types (MPKF)
• Monitoring eating behaviour

7

Dias 13

Modeling of the variable

Dynamic Linear Models (DLMs) combined with Kalman Filter (KF) Let Yt = (y1, … , yn)’ be a vector of key figures observed at time t. Let θt = (θ1, … , θm)’ be a vector of parameters describing the system at time t. General form of the DLM Observation Equation: Yt = F’t θt + νt , νt ~ N(0,Vt) System Equation: θt = Gt θt-1 + ωt, ωt ~ N(0,Wt) DLM combined with Kalman Filter: estimate the underlying state vector θt by its mean vector mt and its variance-covariance matrix Ct. Elements from KF used in monitoring deviations:

• ft : One step forecast mean
• et : One step forecast error (et = Yt – ft )
• Qt : One step forecast variance

Dias 14

Monitoring Deviations from the model

V-mask (parameters d and Ψ) Applied on the cumulative sum (cusum) of the standardized errors ut = et/√Qt

1 1 − =

+ = =∑

t t t t t t

c u u C

Tabular Cusum (parameters K and H) Create a cusum: accumulate ui, using a reference value (K) Alarm when cusum exceeds a decision interval (H)

8

Dias 15

Example 1. Monitoring activity level

Context Development of Group housing in EU results of Council Directive 2001/88/EEC Difficulties identifying and accessing individual sow Idea Store data in a chip and transmit info to the farmer’s PC Sensor in the chip allows to monitor activity of the sow Assumption Body Activity of sows is expected to change around the onset of oestrus Objective Develop an automated oestrus detection method for group housed sows using sows’ acceleration measurements Method Use of Dynamic Linear Models to model the sows’ activity Use of control methods that detects model deviations at the onset of oestrus

Dias 16

Oestrus Detection

Oestrus Detection I (from day 4)

• BPT (3 x / day) in the mating section

Sow not inseminated Transfered to gestation section Oestrus Detection II (from day 21)

• BPT (3 x / day) in the gestation section

Golden Standard Detect whether activity pattern changes at onset of

• estrus

Weaning d0 d7 Transfer d10 d30 d21 Acceleration measurements

BPT I BPT II

9

Dias 17

Data collection

Place, Animals, Housing and Feeding

• 1 production herd, March 2005
• 5 sows in group of 100
• 20 days

Activity Measurements

• Acceleration in 2 and 3 dimensions
• Four measurements per second
• Transfer PC via Blue Tooth

Video Recordings

• Four cameras used as web cam

Oestrus Detection

• Golden standard

Detect whether activity pattern changes at onset of oestrus

Dias 18

Definition of the DLM

Use hourly averages of the length of the acceleration vector Yt = acc = √ (accx

2 + accy 2 + accz 2)

Vt = unknow and constant Wt = 0 (In normal condition: no change in activity) Model initialized by mean of Reference Analysis Model observations (Yt) weighted by number of observations per hour Missing observation: et=0

        =

t t

µ θ

( )

, 1 ' =

t

F

I Gt =

10

Dias 19

Tabular Cusum V-mask Cusum

Illustration

Model

Dias 20

Example 2. Monitoring sows’ activity types

Assumption Sow’s behaviour is affected by physiological state / illness

• Oestrus: increase in activity
• Lameness: walking

Objective Develop a method that automatically classify sows’ activity types Model selected activity types using DLM Classify each activity type using a Multi Process Kalman Filter

Daily Anoestrus Daily Oestrus

Accelerometer: measured any time / during whole reproductive cycle

11

Dias 21

Time series and activity types

Activity types Extracts from time series of acceleration are associated to five activity types

• Feeding (FE)
• Rooting (RO)
• Walking (WA)
• Lying sternally (LS)
• Lying Laterally (LL)

Two data sets Activity filled whole data set / no overlapping Learning data set: 10 minutes of each activity type Estimate the model parameters Test data set: 10 x 2 minutes of each activity type Implement the classification method

X Y Z ACC

3 dimensions: X,Y, Z ACC = √ (accx

2 + accy 2 + accz 2)

Dias 22

Acceleration data

12

Dias 23

Modeling each activity type

Use averages per second of acceleration data

          =

t t t t

c s µ θ

                    = t T t T F t π π 2 cos , 2 sin , 1 '

I Gt =

Model includes a periodic movement – cyclic components V and W: estimated using the EM algorithm Learning data set 20 DLMs 5 activities x 4 axes (X, Y, Z, ACC)

Dias 24

Classification method

Multi Process Kalman Filter of class I

Each DLM is analysed using the updating equations of the Kalman Filter:

• One step forecast mean ft
• One step forecast variance Qt

At time t: Posterior Probabilities are estimated for each DLM

) ( ) ( ) (

1 i

p i i p

t t t −

× ∝φ

13

Classification method

Sted og dato (Indsæt --> Diasnummer) Dias 25

Illustration: walking activity

FE: best recognized WA: axis Z better RO: slow recognition LL: axis Y better LS: well recognized Axis z

Results

Application? Active vs. Passive Illness Parturition

Perspectives

Dias 26

Example 3. Modeling Eating Behaviour

Assumption Sow’s feeding behaviour is affected by oestrus and illness Currently: list of sows that have not eaten is used to identify individuals Objective Develop a method that automatically detect oestrus, lameness and

• ther health disorders for sows fed by ESF
• Model feeding behaviour (Feeding rank) using DLM
• Detect deviations by mean of control chart

14

Dias 27

Data collection

Place, Housing and Feeding

• 3 production herds, January 2005 – January 2006
• Herds 1 and 2: dynamic groups
• Herd 3: static groups
• Electronic Sow Feeders (ESF)

Registration

• ESF
• Oestrus (BPT)
• Lameness and Health disorders

Dias 28

Modeling of the variable

Use daily Feeding Rank (Yt) Vt assumed unknow and constant Wt estimated by discount factor Missing observation : et=0 External information: subgroup of sows enters or leaves a group or both : Lower discounting

        =

t t t

β µ θ

( )

, 1 ' =

t

F

        = 1 1 1

t

G

15

Dias 29

Detection method

Optimization of V-mask parameters for 3 conditions i) oestrus ii) lameness iii) Other health disorders Criteria

• Sensitivity of at least 50%
• Number of FP is minimum

Dias 30

Illustration - Intervention

Individual eating rank Model forecast Subgroup in Subgroup in + out Increase Adaptive Coefficient

16

Dias 31

Results

i) Oestrus detection Sensitivity ranges from 59 to 75% (vs. List of sows: 9 to 20%) ii) Lameness and iii) Other health disorders Sensitivity ranges from 41 to 70% (vs. List of sows: 22 to 39%) Too many false alarms

Perspectives

• Include other variables: e.g. ear base temperature, activity

Multivariate model

Dias 32

Example Daily Gain (from first lecture)

Matrices Specification

Include Seasonal components

17

Dias 33

Example Daily Gain (from first lecture)

Daily gain, slaughter pigs

600 650 700 750 800 850 900 950 2 . k v a r t a l 9 7 3 . k v a r t a l 9 7 4 . k v a r t a l 9 7 1 . k v a r t a l 9 8 2 . k v a r t a l 9 8 3 . k v a r t a l 9 8 4 . k v a r t a l 9 8 1 . k v a r t a l 9 9 2 . k v a r t a l 9 9 3 . k v a r t a l 9 9 4 . k v a r t a l 9 9 1 . k v a r t a l 2 . k v a r t a l 3 . k v a r t a l 4 . k v a r t a l 1 . k v a r t a l 1 Period g Observed gain Predicted gain

Daily gain, slaughter pigs

600 650 700 750 800 850 900 950 2 . k v a r t a l 9 7 3 . k v a r t a l 9 7 4 . k v a r t a l 9 7 1 . k v a r t a l 9 8 2 . k v a r t a l 9 8 3 . k v a r t a l 9 8 4 . k v a r t a l 9 8 1 . k v a r t a l 9 9 2 . k v a r t a l 9 9 3 . k v a r t a l 9 9 4 . k v a r t a l 9 9 1 . k v a r t a l 2 . k v a r t a l 3 . k v a r t a l 4 . k v a r t a l 1 . k v a r t a l 1 Period g Observed gain Predicted gain Level Season 1 Season 2 Season 3 Season 4

Daily gain, slaughter pigs

• 100
• 80
• 60
• 40
• 20

20 40 60 80 100

• 2. kvartal 97
• 3. kvartal 97
• 4. kvartal 97
• 1. kvartal 98
• 2. kvartal 98
• 3. kvartal 98
• 4. kvartal 98
• 1. kvartal 99
• 2. kvartal 99
• 3. kvartal 99
• 4. kvartal 99
• 1. kvartal 00
• 2. kvartal 00
• 3. kvartal 00
• 4. kvartal 00
• 1. kvartal 01

Period g Season 1 Season 2 Season 3 Season 4

Daily gain, slaughter pigs

• 100
• 80
• 60
• 40
• 20

20 40 60 80 100 2 . k v a r t a l 9 7 3 . k v a r t a l 9 7 4 . k v a r t a l 9 7 1 . k v a r t a l 9 8 2 . k v a r t a l 9 8 3 . k v a r t a l 9 8 4 . k v a r t a l 9 8 1 . k v a r t a l 9 9 2 . k v a r t a l 9 9 3 . k v a r t a l 9 9 4 . k v a r t a l 9 9 1 . k v a r t a l 2 . k v a r t a l 3 . k v a r t a l 4 . k v a r t a l 1 . k v a r t a l 1 Period g Forecast error Lower limit Upper limit

Trend Seasons Errors

Dias 34

Concluding remarks

Differents Models were presented

• Simple local level model
• DLM in its general form
• Examples

The general form of the model allow to include cyclic pattern (as for eating activity, daily gain) Thomas Nejsum Madsen will present an approach based on sine functions to incorporate a diurnal pattern. Not necessarily as graphs – automatic alarms (as V mask). Many handles to adjust – dangerous Always combine with your knowledge on animal production