Modélisation du coefficient apparent d'utilisation de l'azote issu d'un engrais minéral apporté sur blé tendre d'hiver

Annexe V : Algorithme de Backfitting

1. Initialisation :

^àâ0 =

Y,

^àfj =

àf(^0) j

2. Itération : estimer tout àtour fk par le modèle

??Y - ^àâ0 -

Xp j=1,j6=k

? ( àfj (xj) x(k)) ? = fk + å

3. Arrêt : l'itération laisse inchangée l'estimation.

vi

Annexe VI : Algorithme de Newton, annexe du cours de P.A. Cornillon

Pierre AndréCornillon écrit :

»Nous rappelons brièvement le principe de cet algorithme. L'étude de sa convergence

et des vitesses associées peuvent être trouver dans un livre d'optimisation. Soit g une fonction àminimiser, par exemple l'opposéde la log-vraisemblance.

g : R^s -> R

0 i-> g(0).

On cherche la valeur de qui réalise le minimum de g, valeur que nous noterons ^à0. Pour cela on utilise une procédure itérative qui va procéder par approximations successives. D'abord, on prend un point de départ 00 puis on chercher autour de celui ci une valeur 01 qui fasse diminuer la fonction g (ie g (01) < g (00)).

Plaçons nous à l'étape k + 1. Nous avons donc 0k qui est connu et nous souhaitons trouver 0k+1 au voisinage de. Comme nous ne savons pas minimiser g, on va l'approxi-mer par une fonction au voisinage de 0k. Cette approximation est choisie de telle sorte qu'un minimum est directement calculable. On va donc approximer la fonction g par une fonction quadratique, grâce àune développement de Taylor :

g (0) ^ g (0k) + (0 - 0k)^' Vg + ₂ (0 - 0k)' V² (0 - 0k)

1

oùles vecteur des dérivées sont évalués au point 0k.

Remplaçons maintenant g(0) par h(0) = g (0k)+(0 - 0k)^' Vg+¹₂ (0 - 0k)' V² (0 - 0k) et cherchons en donc le minimum. Comme il s'agit d'une fonction quadratique, l'unique minimum est le point 0k+1 qui annule la dérivée de h par rapport à 0. La dérivée est donc :

h^'(0) = Vg + V² (0 - 0k)

et elle s'annule au point 0k+1 :

Vg + V² (0 - 0k) = 0

0k+1 = 0k - (V2)-1 Vg

oùchaque dérivée est évaluée au point 0k.»

web.supagro.inra.fr/partage/cornillo/COURS/MODELES-MIXTES/mixtes.pdf

consultéle 23 août 2012.

vii

Annexe VII : Déterminer la structure de variance-covariance d'un effet aléatoire : réponse de Douglas Bates à Michael Ku-bovy sur le R-SIG

Réponse de Douglas Bates, le 7 avril 2008,

https://stat.ethz.ch/pipermail/r-sig-mixed-models/2008q2/000789.html

On Mon, Apr 7, 2008 at 12 :34 PM, Michael Kubovy <kubovy at virginia.edu> wrote :

> Dear lme4 folk

> The lmer help page gives two examples :

> (fm1 <- lmer(Reaction Days + (Days|Subject), sleepstudy))

> (fm2 <- lmer(Reaction Days + (1|Subject) + (0+Days|Subject), sleepstudy))

> How is the following different, in principle, from the above? Is it that the

> above treats (Intercept) and Days as orthogonal, whereas the latter checks

> to see if they are? What would be appropriate if the correlation between

> Days and Intercept (here 0.067, apparently) were large?

> (fm3 <- lmer(Reaction Days + (1 + Days | Subject), sleepstudy))

Model fm3 is equivalent to model fm1. In the linear model formula language used in the S language, the intercept term is implicit so the random-effects term (Days|Subject) is equivalent to (1+Days|Subject). Some authors, notably Gelman and Hill in their 2007 book, prefer to use the second form so that the presence of the intercept is explicit. I can see the point of that.

Every random effect is associated with one and only one random-effects term in the model formula and with one and only one level of the grouping factor for that random-effects term. The general rules for determining the variance-covariance of the random effects (as fit in lmer) are :

- random effects associated with different terms are independent

- random effects associated with the same term but with different levels of the grouping factor are independent

- within a term the random effects may be partitioned according to the levels of the

grouping factor. The variance-covariance matrix of the vector of random effects associated with each of these levels of the grouping factor is a constant, symmetric, positive semidefinite matrix. It has no additional constraints other than being symmetric and positive semidefinite. (In SAS-speak this is called an »unstructu-red»variance-covariance matrix but the mathematician in me refuses to accept the concept of an unstructured, symmetic, positive semidefinite matrix.)

(Note that when I refer to »levels» in the above description I am referring to the S-

language concept of the levels of a factor, not levels of random effects in the sense of

multilevel models.)

In practice the difference between the two models is that fm2 is a restricted form of

fm1/fm3 in which the correlation of the random effects has been set to zero.

> Random effects :

> Groups Name Variance Std.Dev. Corr

> Subject (Intercept) 610.8 24.72

> Days 35.1 5.92 0.067

> Residual 655.1 25.59

> Professor Michael Kubovy > University of Virginia

viii

> Department of Psychology

> USPS : P.O.Box 400400 Charlottesville, VA 22904-4400

> Parcels : Room 102 Gilmer Hall

> McCormick Road Charlottesville, VA 22903

> Office : B011 [phone]

> Lab : B019 [phone]

> Fax: [phone]

> WWW : http ://www.people.virginia.edu/ mk9y/

ix

Annexe VIII : Fonctions d'estimation des RMSEP pour des modèles de classe R nls ou mer

function (lmer)

{

system.time({

RMSEP <- data.frame(iIndex = 1, iDonnees = 1, realiseDeI = 1,

predictionMoinsI = 1, ecartPrediction = 1)

RMSEP <- cbind(RMSEP, t(data.frame(lmer@fixef)))

jdd <- lmer@frame

X0 <- lmer@X

nlignes <- nrow(jdd)

for (i in 1:nlignes) {

updateLmer <- update(lmer, data = jdd[-i, ], REML = T)

predit <- as.numeric(X0[i, ] %*% updateLmer@fixef)

realise <- lmer@y[i]

RMSEP[i, ] <- c(i, rownames(jdd[i, ]), realise, predit,

(realise - predit)^2, updateLmer@fixef)

cat(round(i/nlignes * 100))

}

res.EQMP = sum(as.numeric(RMSEP$ecartPrediction))/nlignes

res.REQMP = sqrt(res.EQMP)

cat("----")

res = list(simulations = RMSEP, estimation = c(EQMP = res.EQMP,

REQMP = res.REQMP))

print(res$estimation)

})

return(res)

}

function (nls, reponse, data)

{

resultat = cbind(data.frame(ind = 1, predit = 1, realise = 1,

ecartSq = 1), t(data.frame(coefficients(nls))))

for (i in 1:nrow(data)) {

modeleMoinsUn <- update(nls, data = data[-i, ])

predit = predict(modeleMoinsUn, newdata = data[i, ])

realise = data[i, reponse]

ecartSq = (predit - realise)^2

resultat[i, ] = c(rownames(data[i, ]), predit, realise,

ecartSq, coefficients(modeleMoinsUn))

cat(round(i/nrow(data) * 100), "")

}

cat("%\n")

EQMP = sum(as.numeric(resultat$ecartSq))/nrow(data)

REQMP = sqrt(EQMP)

res = list(estimations = resultat, validation = c(BIC = BIC(nls),

AIC = AIC(nls), EQMP = EQMP, REQMP = REQMP))

return(res)

}