Modèles non-linéaire

Calcul de variance résiduelle par leave-one-out

Description

On s'intérese à la repousse des pâturages. Les variable sont

• Temps: écoulé après pâture
• Production: mesure de la repousse d'herbage

t = c(9,14,21,28,42,57,63,70,79)
p = c(8.93,10.8,18.59,22.33,39.35,56.11,61.73,64.62,67.08)

p=b₁-b₂*exp(-exp(b₃+b₄*log(t))) + bruit.

N'ayant que 9 individus, on estimera la variance du bruit par validation croisée (PRESS).

Analyse

1/ Estimer le modèle au maximum de vraisemblance.

Il faut utiliser la routine nls de la bibliothèque nls. La difficulté est de trouver des valeurs initiales pour (b₁,b₂,b₃,b₄). Pour cela, on commencera par proposer une valeur initiale pour (b₁,b₂), basée sur les données (interpréter ces deux paramètres). Ces deux valeurs étant fixées, on estimera (b₃,b₄) avec nls(.) en partant de la valeur initiale (0, 0.01). Le quadruplet formé des valeurs utilisées pour (b₁,b₂) et des valeurs obtenues pour (b₃,b₄) constituera un bon point initial pour le calcul du vecteur de paramètres complet. Pour réintroduire l'estimée de (b₃,b₄) dans les conditions initiales, on pourra convertir la solution obtenue avec getAllPars() en un simple vecteur de réels en utilisant as.real().

2/ Tracer les points et la courbe de régression.

Pour éviter l'effacement, on utilisera plot() d'abord pour les points de données, et ensuite lines() avec les 9 x_i et les predictions correspondantes.

3/ Estimées de sigma

Estimer les residus par validation croiseée et en déduire l'estimation CV de sigma.
La comparer à l'estimateur MV.
Indication : pour chacune des optimisations à faire, on donnera comme valeur initiale des paramètres celles obtenues à la première question.

Une solution ici

Source

L'âge des lapins

Description

The European rabbit Oryctolagus cuniculus is a major pest in Australia. A reliable method of age determination for rabbits caught in the wild would be of importance in ecological studies. In this study, the dry weight of the eye lens was measured for 71 free-living wild rabbits of known age. Eye lens weight tends to vary much less with environmental conditions than does total body weight, and therefore may be a much better indicator of age.

The rabbits were born and lived free in an experimental 1.7 acre enclosure at Gungahlin, ACT. The birth data and history of each individual were accurately known. Rabbits in the enclosure depended on the natural food supply. In this experiment, 18 of the eye lenses were collected from rabbits that died in the course of the study from various causes such as coccidiosis, bird predation or starvation. The remaining 53 rabbits were deliberately killed, immediately after being caught in the enclosure or after they had been kept for some time in cages. The lenses were preserved and their dry weight determined.

• Age: Age of rabbit in days
• Lens:Dry weight of eye lens in milligrams

Analyse

1/ Lire les données.

2/ Choix du modèle.

Au vu de tracés des données, discuter les modèles

Lens = aexp{ -b/ (Age + c) } + bruit

et

log(Lens) =a-b/ (Age + c) + bruit

3/  Estimation.

Pour l'estimation avec la fonction nls, on observera sans doute que l'algorithme diverge. Pour remédier à cela, on pourra soit

(a) deviner des valeurs initiales raisonnables pour a, b et c au vu des données (penser à la signification de ces trois quantités) : mod= nls(y~a-b/(Age+c),start=c(a=?,b=?,c=?))

(b) soit utiliser l'option algorithm="plinear". Par exemple un modèle y=b₁+b₂ x+b₃ exp(b₄ x) se définira par la formule

mod =nls(y~cbind(1,x,exp(b4 x)),start=c(b4=?),algorithm="plinear")

où l'on ne donne de valeur initiale que pour b₄. Ainsi l'algorithme exploitera la linéarité en (b₁,b₂,b₃) et aura davantage de chances de converger.

Tracer ensuite la courbe de régression et les données sur un même graphique. Pour éviter l'effacement, on utilisera plot() d'abord pour les points de données, et ensuite lines() avec les régresseurs des données et les predictions correspondantes (utiliser fitted(mod)).

4/ Calcul de la variance MV.

Comparer les résultats de sum(mod$m$resid()^2) et mod$m$deviance(). En déduire une estimée de la variance.

5/ Calcul de la variance par validation croisée.

Il faut faire une boucle.
On ôtera la donnée i par la commande lapi=lap[row.names(lap)!=i,].
Pour l'estimation on donnera le modèle précédemment estimé comme valeur initiale (utiliser mod$m$getAllPars()).
On prédira avec la fonction predict avec comme argument newdata la liste list(Age=lap$Age[i]).

Sources

Available Chlorine (Bootstrap)

Description

Chlorine = a + (0.49 - a) exp{ -b (Weeks - 8) }

• Weeks: Time in weeks since manufacture
• Chlorine: Available chlorine

Analyse

Pour cela, on identifiera le modèle y=b₁+b₂*exp(-b₃*x) avec l'option algorithm=«plinear» (cf.TP «lapins» point 3.(b) ci-dessus). On vérifiera ensuite que la contriainte donnée a priori b₂*exp(-8*b₃)+b₁=0,49 est réaliste au vu des intervalles de confiance sur les estimées.

On peut travailler ensuite directement sur le modèle un peu plus compliqué à estimer: y=a+(0,49-a)*exp(-b*(x-8)).

2/ Tracer les points et la courbe de régression.

Pour éviter l'effacement, on utilisera plot() d'abord pour les points de données, et ensuite lines() avec les x_i et les predictions correspondantes.

3/ Estimer la variance d'estimation de chaque b_i par la méthode de boostrap sur les résidus; comparer à celles données par R. Calculer par la même méthode un intervalle de confiance à 95% pour chaque b_i.

Sources

• OzDASL
• Jennrich, R. I., and Bright, P. B. (1976). Fitting systems of linear differential equations using computer generated exact derivatives. Technometrics, 18, 385-392.
• Jennrich, R. I. (1995). An Introduction to Computational Statistics. Prentice-Hall, Englewood Cliffs, New Jersey. Section 8.2.1.