##-------------------------------------------------------------------------------------------------------- ## SCRIPT : TD Teneurs en protéines avec Jags ## ## Authors : ## Authors : Collectif BioBayes ## Last update : 2019-02-25 ## R version 3.5.1 (2018-07-02) -- "Feather Spray" ## Copyright (C) 2018 The R Foundation for Statistical Computing ## Platform: x86_64-w64-mingw32/x64 (64-bit) ##-------------------------------------------------------------------------------------------------------- library(rjags) library(coda) #Les données Y=c(11.9,10.3,12.4,12.1,12.4,12.7,11.5,13.6,9.9,11.1,10.6,12,12.8,11.6,11.2,10.3, 12.9,10.3,11.5,14.1,12.6,12.7,11.9,9.2,12.5,11.9,11.3,11,10,15.9,10.6,11.6, 11.4,10.9,14.6,10.7,9.4,13,12.5,12.9,11.1,12.9,13.2) INN=c(1,0.8,0.8,1,0.7,0.9,1,0.9,0.6,0.7,0.7,0.7,0.9,0.8,0.8,0.6,0.9,0.69,0.72,1.1,0.97,0.98,0.87, 0.62,0.73,1.09,0.76,0.86,0.75,1.03,0.75,0.86,0.82,0.77,0.88,0.84,0.72,0.81,0.74,0.94,0.78,0.92,1.07) SPAD=c(45.5,42.6,47.8,48.6,48.2,48.8,48.9,49.9,48.6,44.9,50.2,46.1,51.6,50.9,40.4, 47.6,49.9,37.2,44.87,48,46.03,45.8,46.87,38.67,42.67,44.23,39.67,39.4,39.43,44.73,38.27,44.6,45.5,46.1, 38,43.6,43.5,44,49.9,51,48.8,48.7,49.6) ###Covariables centrées-réduites INN_cr= (INN-mean(INN))/sd(INN) SPAD_cr= (SPAD-mean(SPAD))/sd(SPAD) ################################ #Modèle à 1 variable explicative ################################ modelString =" model{ # Y= teneur en protéines du blé # covar_cr = variable explicative (SPAD ou INN) centrée réduite # alpha=paramètres de la régression # prec=paramètre de précision = 1/variance # m, s= moyenne et écart-type empirique de covar non centrée réduite #Modèle d'observations for (i in 1:N) { Y[i] ~ dnorm(mu[i],prec) #mu[i]<-alpha_nonc[1]+alpha_nonc[2]*covar[i] mu[i]<-alpha[1]+alpha[2]*covar_cr[i] #Calcul de la log_vraisemblance de l'observation i loglik[i] <- logdensity.norm(Y[i], mu[i], prec) } alpha_nonc[1]<- alpha[1]-alpha[2]*(m/s) alpha_nonc[2]<- alpha[2]/s #Lois a priori for (j in 1:2) { alpha[j]~dnorm(0.0,1.0E-6) } sd ~ dunif(0,100) prec<- 1/pow(sd,2) } " writeLines(modelString, con="proteines_1var.bug") ###Ajustement modèle à 1 variable explicative = INN jags_Inn <- jags.model(file="proteines_1var.bug", data =list(Y=Y,covar_cr=INN_cr,N=43,m=mean(INN),s=sd(INN)), inits =list(list(alpha=c(0,0), sd=1), list(alpha = c(5,5), sd = 0.1), list(alpha = c(-5,-5), sd = 10)), n.chains =3,n.adapt = 100) #Vérifier tout d'abord que les programmes fonctionnent en réalisant 5000 itérations MCMC update(jags_Inn, 5000) #Stocker les valeurs échantillonnées pour trois chaînes de Markov posteriorjags_Inn <- coda.samples(model=jags_Inn,variable.names=c("alpha_nonc", "sd"), n.iter=15000, thin=1, progress.bar="gui") #Utiliser le diagnostic de convergence de Gelman et Rubin et des représentations #graphiques des chaînes pour déterminer le nombre d'itérations nécessaire #pour espérer avoir atteint la convergence. plot(posteriorjags_Inn) gelman.plot(posteriorjags_Inn) #Analyser les auto-corrélations des chaînes. autocorr.plot(posteriorjags_Inn) #Statistiques a posteriori print(summary(posteriorjags_Inn,dig=2)) #Calcul du DIC dic_Inn<- dic.samples(model=jags_Inn, n.iter=10000, thin = 1) dic_Inn #Calcul du WAIC loglike <- coda.samples(model=jags_Inn,variable.names=c("loglik"), n.iter=10000, thin=1, progress.bar="gui") M<- as.matrix(loglike) #Calcul du 1er terme T1=-2*[log(vraisemblance moyenne par individu) puis somme sur les individus] T1<- -2*sum(log(apply(exp(M),FUN=mean,MARGIN=2))) #Calcul de la pénalité T2=2*[variance sur les itérations MCMC de la log-vraisemblance par individu puis somme sur les individus] T2<- 2*sum(apply(M,FUN=var,MARGIN=2)) #Calcul du WAIC WAIC_Inn<- T1+T2 WAIC_Inn ###Ajustement modèle à 1 variables explicative = SPAD (centrée réduite) jags_SPAD <- jags.model(file="proteines_1var.bug", data =list(Y=Y,covar_cr=SPAD_cr,N=43,m=mean(SPAD),s=sd(SPAD)), inits =list(list(alpha=c(0,0), sd=1), list(alpha = c(5,5), sd = 0.1), list(alpha = c(-5,-5), sd = 10)), n.chains =3,n.adapt = 100) #Vérifier tout d'abord que les programmes fonctionnent en réalisant 5000 itérations MCMC update(jags_SPAD, 5000) #Lancer trois chaînes de Markov posteriorjags_SPAD<- coda.samples(model=jags_SPAD,variable.names=c("alpha_nonc", "sd"), n.iter=15000, thin=1, progress.bar="gui") #Utiliser le diagnostic de convergence de Gelman et Rubin et des représentations graphiques #des chaînes pour déterminer le nombre d'itérations nécessaire pour espérer avoir atteint la convergence. plot(posteriorjags_SPAD) gelman.plot(posteriorjags_SPAD) #Analyser les auto-corrélations des chaînes. autocorr.plot(posteriorjags_SPAD) #Statistiques a posteriori print(summary(posteriorjags_SPAD),dig=2) #Calcul du DIC dic_SPAD<- dic.samples(model=jags_SPAD, n.iter=10000, thin = 1) dic_SPAD #Calcul du WAIC loglike <- coda.samples(model=jags_SPAD,variable.names=c("loglik"), n.iter=10000, thin=1, progress.bar="gui") M<- as.matrix(loglike) #Calcul du 1er terme T1=-2*[log(vraisemblance moyenne par individu) puis somme sur les individus] T1<- -2*sum(log(apply(exp(M),FUN=mean,MARGIN=2))) #Calcul de la pénalité T2=2*[variance sur les itérations MCMC de la log-vraisemblance par individu puis somme sur les individus] T2<- 2*sum(apply(M,FUN=var,MARGIN=2)) #Calcul du WAIC WAIC_SPAD<- T1+T2 WAIC_SPAD ################################## #Modèle à 2 variables explicatives ################################## modelString =" model{ # Y= teneur en protéines du blé # INN_cr = Indice de nutrition azotée centré réduit # SPAD_cr = Mesure de transmittance centrée réduite # alpha=paramètres de la régression # prec=paramètre de précision = 1/variance # m_INN, s_INN= moyenne et écart-type empirique de INN # m_SPAD, s_SPAD= moyenne et écart-type empirique de SPAD #Modèle d'observations for (i in 1:N) { Y[i] ~ dnorm(mu[i],prec) #mu[i]<-alpha_nonc[1]+alpha_nonc[2]*SPAD[I]+ alpha_nonc[3]*INN[I] => Problème convergence mu[i]<-alpha[1]+alpha[2]*SPAD_cr[i]+alpha[3]*INN_cr[i] #Calcul de la log_vraisemblance de l'observation i loglik[i] <- logdensity.norm(Y[i], mu[i], prec) } alpha_nonc[1]<- alpha[1]-alpha[2]*(m_SPAD/s_SPAD)-alpha[3]*(m_INN/s_INN) alpha_nonc[2]<- alpha[2]/s_SPAD alpha_nonc[3]<- alpha[3]/s_INN #Lois a priori for (j in 1:3) { alpha[j]~dnorm(0.0,1.0E-6) } sd ~ dunif(0,100) prec<- 1/pow(sd,2) } " writeLines(modelString, con="proteines_2var.bug") ###Ajustement modèle à 2 variables explicatives jags_2var <- jags.model(file="proteines_2var.bug", data =list(Y=Y,SPAD_cr=SPAD_cr,INN_cr=INN_cr,N=43,m_INN=mean(INN),s_INN=sd(INN),m_SPAD=mean(SPAD),s_SPAD=sd(SPAD)), inits =list(list(alpha=c(0,0,0), sd=1), list(alpha = c(5,5,5), sd = 0.1), list(alpha = c(-5,-5,-5), sd = 10)), n.chains =3,n.adapt = 100) #Vérifier tout d'abord que les programmes fonctionnent en réalisant 5000 itérations MCMC update(jags_2var, 5000) #Lancer trois chaînes de Markov posteriorjags_2var <- coda.samples(model=jags_2var,variable.names=c("alpha_nonc", "sd"), n.iter=10000, thin=1, progress.bar="gui") #Utiliser le diagnostic de convergence de Gelman et Rubin et des représentations graphiques #des chaînes pour déterminer le nombre d'itérations nécessaire pour espérer avoir atteint la convergence plot(posteriorjags_2var) gelman.plot(posteriorjags_2var) #Analyser les auto-corrélations des chaînes. autocorr.plot(posteriorjags_2var) #Statistiques a posteriori print(summary(posteriorjags_2var),dig=3) #Calcul du DIC dic_2var<- dic.samples(model=jags_2var, n.iter=10000, thin = 1) dic_2var #Calcul du WAIC loglike <- coda.samples(model=jags_2var,variable.names=c("loglik"), n.iter=10000, thin=1, progress.bar="gui") M<- as.matrix(loglike) #Calcul du 1er terme T1=-2*[log(vraisemblance moyenne par individu) puis somme sur les individus] T1<- -2*sum(log(apply(exp(M),FUN=mean,MARGIN=2))) #Calcul de la pénalité T2=2*[variance sur les itérations MCMC de la log-vraisemblance par individu puis somme sur les individus] T2<- 2*sum(apply(M,FUN=var,MARGIN=2)) #Calcul du WAIC WAIC_2var<- T1+T2 WAIC_2var ################################ #Modèle sans variable explicative ################################ modelString =" model{ # Y= teneur en protéines du blé # mu=paramètre d'espérance # prec=paramètre de précision = 1/variance #Modèle d'observations for (i in 1:N) { Y[i] ~ dnorm(mu,prec) #Calcul de la log_vraisemblance de l'observation i loglik[i] <- logdensity.norm(Y[i], mu, prec) } #Lois a priori mu~dnorm(0.0,1.0E-6) sd ~ dunif(0,100) prec<- 1/pow(sd,2) } " writeLines(modelString, con="proteines_sansvar.bug") ###Ajustement modèle sans variable explicative jags_sansvar <- jags.model(file="proteines_sansvar.bug", data =list(Y=Y,N=43), inits =list(list(mu=0, sd=1), list(mu=5, sd = 0.1), list(mu=-5, sd = 10)), n.chains =3,n.adapt = 100) #Vérifier tout d'abord que les programmes fonctionnent en réalisant 5000 itérations MCMC update(jags_sansvar, 5000) posteriorjags_sansvar <- coda.samples(model=jags_sansvar,variable.names=c("mu", "sd"), n.iter=10000, thin=1, progress.bar="gui") #Utiliser le diagnostic de convergence de Gelman et Rubin et des représentations graphiques #des chaînes pour déterminer le nombre d'itérations nécessaire pour espérer avoir atteint la convergence. plot(posteriorjags_sansvar) gelman.plot(posteriorjags_sansvar) #Analyser les auto-corrélations des chaînes. autocorr.plot(posteriorjags_sansvar) #Statistiques a posteriori print(summary(posteriorjags_sansvar),dig=3) #Calcul du DIC dic_sansvar<- dic.samples(model=jags_sansvar, n.iter=10000, thin = 1) dic_sansvar #Calcul du WAIC loglike <- coda.samples(model=jags_sansvar,variable.names=c("loglik"), n.iter=10000, thin=1, progress.bar="gui") M<- as.matrix(loglike) #Calcul du 1er terme T1=-2*[log(vraisemblance moyenne par individu) puis somme sur les individus] T1<- -2*sum(log(apply(exp(M),FUN=mean,MARGIN=2))) #Calcul de la pénalité T2=2*[variance sur les itérations MCMC de la log-vraisemblance par individu puis somme sur les individus] T2<- 2*sum(apply(M,FUN=var,MARGIN=2)) #Calcul du WAIC WAIC_2var<- T1+T2 WAIC_2var ######################################################################## #Prediction a posteriori selon le modèle à 1 variable explicative = INN ########################################################################### ##Le script "proteines_1var.bug" du modèle univarié reste inchangé. INNpred_cr<- (0.95-mean(INN))/sd(INN) ###Ajustement + Prediction modèle à 1 variables explicative = INN jags_Inn_pred <- jags.model(file="proteines_1var.bug", data =list(Y=c(Y,NA),covar_cr=c(INN_cr,INNpred_cr),N=44,m=mean(INN), s=sd(INN)), inits =list(list(Y=c(rep(NA,times=43),5),alpha=c(0,0), sd=1), list(Y=c(rep(NA,times=43),15),alpha = c(5,5), sd = 0.1), list(Y=c(rep(NA,times=43),10),alpha = c(-5,-5), sd = 10)), n.chains =3,n.adapt = 100) update(jags_Inn_pred, 5000) posteriorjags_Inn_pred <- coda.samples(model=jags_Inn_pred,variable.names=c("Y[44]"), n.iter=10000, thin=1, progress.bar="gui") plot(posteriorjags_Inn_pred[,1]) print(summary(posteriorjags_Inn_pred[,1]),dig=3)