L’une des hypothèses de base de la modélisation linéaire est la variance constante, ou homogène. Qu’est-ce que cela signifie exactement ? Simulons quelques données qui satisfont à cette condition pour illustrer le concept.
Ci-après, nous créons un vecteur trié de nombres allant de 1 à 10 appelé x, puis nous créons un vecteur de nombres appelé y qui est une fonction de x. Lorsque nous traçons x vs y, nous obtenons une ligne droite avec un intercept de 1,2 et une pente de 2,1.
x <- seq(1,10, length.out = 100)y <- 1.2 + 2.1 * xplot(x, y)
Maintenant, ajoutons un peu de « bruit » à nos données afin que y ne soit pas complètement déterminé par x. Nous pouvons le faire en tirant au hasard des valeurs d’une distribution normale théorique avec une moyenne de 0 et une certaine variance définie, puis en les ajoutant à la formule qui génère y. La fonction rnorm de R nous permet de le faire facilement. Ci-dessous, nous tirons 100 valeurs aléatoires d’une distribution normale avec une moyenne de 0 et un écart-type de 2 et les enregistrons comme un vecteur appelé noise. (Rappelez-vous que l’écart-type est simplement la racine carrée de la variance.) Ensuite, nous générons y avec le bruit ajouté. La fonction set.seed(222) vous permet d’obtenir les mêmes données « aléatoires » au cas où vous voudriez suivre le mouvement. Enfin, nous re-tracons les données.
set.seed(222)noise <- rnorm(n = 100, mean = 0, sd = 2)y <- 1.2 + 2.1 * x + noiseplot(x, y)
Nous avons maintenant des données qui sont une combinaison d’une composante linéaire, déterministe (\(y = 1,2 + 2,1x\\)) et d’un bruit aléatoire tiré d’une distribution \(N(0, 2)\). Ce sont les hypothèses de base que nous faisons sur nos données lorsque nous ajustons un modèle linéaire traditionnel. Ci-dessous, nous utilisons la fonction lm pour « récupérer » les « vraies » valeurs que nous avons utilisées pour générer les données, qui sont au nombre de trois :
- L’intercept : 1,2
- La pente : 2.1
- L’écart type : 2
m <- lm(y ~ x)summary(m)## ## Call:## lm(formula = y ~ x)## ## Residuals:## Min 1Q Median 3Q Max ## -5.5831 -1.2165 0.3288 1.3022 4.3714 ## ## Coefficients:## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 1.27426 0.44720 2.849 0.00534 ** ## x 2.09449 0.07338 28.541 < 2e-16 ***## ---## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## ## Residual standard error: 1.926 on 98 degrees of freedom## Multiple R-squared: 0.8926, Adjusted R-squared: 0.8915 ## F-statistic: 814.6 on 1 and 98 DF, p-value: < 2.2e-16Les estimations de (Intercept) et x (c’est-à-dire la pente) dans la section Coefficients, 1,27 et 2,09, sont assez proches de 1,2 et 2,1, respectivement. L’erreur standard résiduelle de 1,926 est également assez proche de la valeur constante de 2. Nous avons produit un « bon » modèle parce que nous savions comment y était généré et que nous avons donné à la fonction lm le « bon » modèle à ajuster. Bien que nous ne puissions pas faire cela dans la vie réelle, c’est un exercice utile pour nous aider à comprendre les hypothèses de la modélisation linéaire.
Maintenant, et si la variance n’était pas constante ? Et si nous multipliions l’écart-type de 2 par la racine carrée de x ? Comme nous le voyons dans le graphique ci-dessous, la dispersion verticale des points augmente lorsque x augmente.
set.seed(222)noise <- rnorm(n = 100, mean = 0, sd = 2 * sqrt(x))y <- 1.2 + 2.1 * x + noiseplot(x, y)
Nous avons multiplié 2 par sqrt(x) parce que nous spécifions l’écart-type. Si nous pouvions spécifier la variance, nous aurions multiplié 4 par seulement x.
Si nous ajustons le même modèle en utilisant lm, nous obtenons une erreur standard résiduelle de 4,488.
m2 <- lm(y ~ x)summary(m2)## ## Call:## lm(formula = y ~ x)## ## Residuals:## Min 1Q Median 3Q Max ## -15.0460 -2.4013 0.5638 2.8734 10.2996 ## ## Coefficients:## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 1.313 1.042 1.26 0.211 ## x 2.096 0.171 12.26 <2e-16 ***## ---## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## ## Residual standard error: 4.488 on 98 degrees of freedom## Multiple R-squared: 0.6051, Adjusted R-squared: 0.6011 ## F-statistic: 150.2 on 1 and 98 DF, p-value: < 2.2e-16Nous savons que ce n’est pas correct, car nous avons simulé les données. Il n’y avait pas d’écart-type constant lorsque nous avons créé le « bruit ». Chaque valeur aléatoire a été tirée d’une distribution normale différente, chacune ayant une moyenne de 0 mais un écart-type qui variait selon x. Cela signifie que notre hypothèse de variance constante est violée. Comment pourrions-nous détecter cela dans la vie réelle ?
La façon la plus courante est de tracer les résidus par rapport aux valeurs ajustées. Ceci est facile à faire dans R. Il suffit d’appeler plot sur l’objet modèle. Cela génère quatre tracés différents pour évaluer les hypothèses de modélisation traditionnelles. Voir cet article de blog pour plus d’informations. Les tracés qui nous intéressent sont les 1er et 3e tracés, que nous pouvons spécifier avec l’argument which.
plot(m2, which = 1)
plot(m2, which = 3)
Dans le premier tracé, nous voyons la variabilité autour de 0 augmenter à mesure que nous nous déplaçons vers la droite avec des valeurs ajustées plus grandes. Dans le troisième tracé, nous voyons également une variabilité croissante à mesure que nous nous déplaçons vers la droite, bien que cette fois les résidus aient été normalisés et transformés à l’échelle de la racine carrée. La trajectoire positive de la ligne rouge lisse indique une augmentation de la variance.
Alors, maintenant que nous avons confirmé que notre hypothèse de variance non constante est invalide, que pouvons-nous faire ? Une approche consiste à transformer les données en logarithme. Cela fonctionne parfois si votre variable de réponse est positive et fortement asymétrique. Mais ce n’est pas vraiment le cas ici. y n’est que légèrement asymétrique. (Appelez hist() sur y pour vérifier.) De plus, nous savons que nos données ne sont pas sur une échelle logarithmique.
Une autre approche consiste à modéliser la variance non constante. C’est le sujet de ce billet de blog.
Pour ce faire, nous allons utiliser les fonctions du paquet nlme qui est inclus dans l’installation de base de R. La fonction cheval de bataille est gls, qui signifie « moindres carrés généralisés ». Nous l’utilisons comme nous le ferions avec la fonction lm, sauf que nous utilisons également l’argument weights ainsi qu’une poignée de fonctions de variance. Faisons une démonstration avec un exemple.
Ci-après, nous ajustons le modèle « correct » à nos données qui présentaient une variance non constante. Nous chargeons le package nlme pour pouvoir utiliser la fonction gls1. Nous spécifions la syntaxe du modèle comme précédemment, y ~ x. Puis nous utilisons l’argument weights pour spécifier la fonction de variance, dans ce cas varFixed, qui fait également partie du package nlme. Cela indique que notre fonction de variance n’a pas de paramètres et une seule covariable, x, ce qui est précisément la façon dont nous avons généré les données. La syntaxe, ~x, est une formule unilatérale qui peut être lue comme « variance du modèle en fonction de x. »
library(nlme)vm1 <- gls(y ~ x, weights = varFixed(~x))summary(vm1)## Generalized least squares fit by REML## Model: y ~ x ## Data: NULL ## AIC BIC logLik## 576.2928 584.0477 -285.1464## ## Variance function:## Structure: fixed weights## Formula: ~x ## ## Coefficients:## Value Std.Error t-value p-value## (Intercept) 1.369583 0.6936599 1.974431 0.0511## x 2.085425 0.1504863 13.857900 0.0000## ## Correlation: ## (Intr)## x -0.838## ## Standardized residuals:## Min Q1 Med Q3 Max ## -2.8942967 -0.6293867 0.1551594 0.6758773 2.2722755 ## ## Residual standard error: 1.925342 ## Degrees of freedom: 100 total; 98 residualLe résultat produit une erreur standard résiduelle de 1,925 qui est très proche de 2, la « vraie » valeur que nous avons utilisée pour générer les données. Les valeurs de (Intercept) et de pente, 1,37 et 2,09, sont également très proches de 1,2 et 2,1.
Une fois encore, nous pouvons appeler plot sur l’objet modèle. Dans ce cas, un seul graphique est généré : résidus standardisés en fonction des valeurs ajustées :
plot(vm1)
Ce graphique semble bon. Tant que nous modélisons notre variance en fonction de x, le modèle ajusté ne surajuste ni ne sous-ajuste d’aucune sorte systématique (contrairement à ce qui se passait lorsque nous utilisions lm pour ajuster le modèle et supposions une variance constante.)
La fonction varFixed crée des poids pour chaque observation, symbolisés par \(w_i\). L’idée est que plus le poids d’une observation donnée est élevé, plus la variance de la distribution normale dont a été tirée sa composante de bruit est faible. Dans notre exemple, la variance augmente à mesure que x augmente. Par conséquent, des poids plus élevés devraient être associés à des valeurs x plus petites. Ceci peut être accompli en prenant l’inverse de x, c’est-à-dire \(w_i = 1/x\). Donc, lorsque \(x = 2\), \(w = 1/2\). Quand \(x = 10\), \(w = 1/10\). Les plus grands x obtiennent des poids plus petits.
Enfin, pour s’assurer que les plus grands poids sont associés à de plus petites variances, nous divisons la variance constante par le poids. Exprimé mathématiquement,
\
Ou en prenant la racine carrée pour exprimer l’écart type,
\
Donc, plus le dénominateur est grand (c’est-à-dire, plus le poids est grand et donc, plus le x est petit), plus la variance est petite et plus la y observée est précise.
Incidemment, nous pouvons également utiliser lm pour pondérer les observations. Elle aussi possède un weights argument comme la fonction gls. La seule différence est que nous devons être plus explicites sur la façon dont nous exprimons les poids. Dans ce cas, nous devons spécifier la réciproque de x. Remarquez que le résultat est presque identique à celui que nous obtenons en utilisant gls et varFixed.
m3 <- lm(y ~ x, weights = 1/x)summary(m3)## ## Call:## lm(formula = y ~ x, weights = 1/x)## ## Weighted Residuals:## Min 1Q Median 3Q Max ## -5.5725 -1.2118 0.2987 1.3013 4.3749 ## ## Coefficients:## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 1.3696 0.6937 1.974 0.0511 . ## x 2.0854 0.1505 13.858 <2e-16 ***## ---## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## ## Residual standard error: 1.925 on 98 degrees of freedom## Multiple R-squared: 0.6621, Adjusted R-squared: 0.6587 ## F-statistic: 192 on 1 and 98 DF, p-value: < 2.2e-16La puissance du package nlme est qu’il permet une variété de fonctions de variance. La fonction varFixed que nous venons d’illustrer est la plus simple et quelque chose qui peut être fait avec lm. Les autres fonctions de variance comprennent :
varIdentvarPowervarExpvarConstPowervarComb
Explorons chacune d’entre elles en utilisant des données simulées.
varIdent
La fonction varIdent nous permet de modéliser différentes variances par strate. Pour voir comment cela fonctionne, nous allons d’abord simuler des données avec cette propriété. Ci-dessous, nous utilisons set.seed(11) au cas où quelqu’un voudrait simuler les mêmes données aléatoires. Puis nous mettons n égal à 400, le nombre d’observations que nous allons simuler. x est généré de la même manière que précédemment. Dans cet exemple, nous incluons un prédicteur supplémentaire appelé g qui peut être considéré comme le sexe. Nous échantillonnons aléatoirement 400 valeurs de « m » et « f ». Ensuite, nous simulons un vecteur de deux écarts types et l’enregistrons comme msd. Si g == "m", alors l’écart-type est \(2 \times 2.5\). Sinon, c’est \(2\) pour g == "f". Nous utilisons ce vecteur dans la ligne suivante pour générer y. Remarquez où msd est branché dans la fonction rnorm. Remarquez également comment nous générons y. Nous incluons une interaction entre x et y. Lorsque g == "f", l’intercept et la pente sont de 1,2 et 2,1. Lorsque g == "m", l’intercept et la pente sont (1,2 + 2,8) et (2,1 + 2,8), respectivement. Enfin, nous plaçons nos données simulées dans un cadre de données et traçons avec ggplot2.
set.seed(11)n <- 400x <- seq(1,10, length.out = n)g <- sample(c("m","f"), size = n, replace = TRUE)msd <- ifelse(g=="m", 2*2.5, 2)y <- 1.2 + 2.1*x - 1.5*(g == "m") + 2.8*x*(g == "m") + rnorm(n, sd = msd)d <- data.frame(y, x, g)library(ggplot2)ggplot(d, aes(x, y, color = g)) + geom_point()
Notez les différentes variances pour chaque niveau de g. La variance pour « m » est beaucoup plus grande que la variance pour « f ». Elle a beaucoup plus de dispersion que « f ». Nous avons simulé les données de cette façon. Nous avons fixé la variance de « m » à 2,5 fois celle de « f ».
Maintenant, utilisons la fonction gls avec varIdent pour tenter de retrouver ces vraies valeurs. Nous utilisons la même méthode que précédemment : définir notre fonction de variance dans l’argument weights. Ci-dessous, nous précisons dans l’argument form que la formule de modélisation de la variance est conditionnelle à g. L’expression ~ 1|g est une formule unilatérale qui dit que la variance diffère entre les niveaux de g. L’expression 1 signifie simplement que nous n’avons pas de covariables supplémentaires dans notre modèle de variance. Il est possible d’inclure une formule telle que ~ x|g, mais ce serait incorrect dans ce cas puisque nous n’avons pas utilisé x pour générer la variance. L’examen du graphique montre également que si la variabilité de y est différente entre les groupes, la variabilité n’augmente pas lorsque x augmente.
vm2 <- gls(y ~ x*g, data = d, weights = varIdent(form = ~ 1|g))summary(vm2)## Generalized least squares fit by REML## Model: y ~ x * g ## Data: d ## AIC BIC logLik## 2081.752 2105.64 -1034.876## ## Variance function:## Structure: Different standard deviations per stratum## Formula: ~1 | g ## Parameter estimates:## f m ## 1.00000 2.58127 ## ## Coefficients:## Value Std.Error t-value p-value## (Intercept) 0.9686349 0.3237724 2.99172 0.0029## x 2.1222707 0.0525024 40.42239 0.0000## gm -1.9765090 0.9352152 -2.11343 0.0352## x:gm 2.9957974 0.1553551 19.28355 0.0000## ## Correlation: ## (Intr) x gm ## x -0.901 ## gm -0.346 0.312 ## x:gm 0.304 -0.338 -0.907## ## Standardized residuals:## Min Q1 Med Q3 Max ## -2.74115039 -0.67013954 0.01619031 0.69793776 2.72985748 ## ## Residual standard error: 2.005397 ## Degrees of freedom: 400 total; 396 residualPremièrement, remarquez au bas de la sortie que l’erreur standard résiduelle estimée est de 2.0053974, très proche de la « vraie » valeur de 2. Remarquez également dans la section « Fonction de variance » que nous obtenons une valeur estimée de 2,58127 pour le groupe « m », ce qui est très proche de la « vraie » valeur de 2,5 que nous avons utilisée pour générer la variance différente pour le groupe « m ». En général, lorsque vous utilisez la fonction varIdent pour estimer des variances différentes entre les niveaux de strates, l’un des niveaux sera fixé à la ligne de base, et les autres seront estimés comme des multiples de l’erreur standard résiduelle. Dans ce cas, puisque « f » vient avant « m » alphabétiquement, « f » a été fixé à la ligne de base, ou 1. L’erreur type résiduelle estimée pour le groupe « f » est \(2,005397 \times 1\). L’erreur standard résiduelle estimée pour le groupe « m » est \(2,005397 \times 2,58127\).
Il est important de noter que ce ne sont que des estimations. Pour avoir une idée de l’incertitude de ces estimations, nous pouvons utiliser la fonction intervals du package nlme pour obtenir des intervalles de confiance à 95%. Pour réduire la sortie, nous enregistrons le résultat et visualisons les éléments sélectionnés de l’objet.
int <- intervals(vm2)L’élément varStruct contient l’intervalle de confiance à 95% pour le paramètre estimé dans la fonction de variance. Dans ce cas, le paramètre est le multiplicateur d’erreur standard résiduelle pour le groupe des hommes.
int$varStruct## lower est. upper## m 2.245615 2.58127 2.967095## attr(,"label")## "Variance function:"L’élément sigma contient l’intervalle de confiance à 95% pour l’erreur standard résiduelle.
int$sigma## lower est. upper ## 1.818639 2.005397 2.211335 ## attr(,"label")## "Residual standard error:"Les deux intervalles sont assez petits et contiennent la « vraie » valeur que nous avons utilisée pour générer les données.
varPower
La fonction varPower nous permet de modéliser la variance comme une puissance de la valeur absolue d’une covariable. Encore une fois, pour aider à expliquer cela, nous allons simuler des données avec cette propriété. Ci-dessous, la principale chose à remarquer est l’argument sd de la fonction rnorm. C’est là que nous prenons l’écart type de 2, puis le multiplions par la valeur absolue de x élevée à la puissance 1,5. C’est similaire à la façon dont nous avons simulé des données pour démontrer la fonction varFixed ci-dessus. Dans ce cas, nous avons simplement supposé que l’exposant était de 0,5. (Rappelez-vous que prendre la racine carrée d’un nombre est équivalent à l’élever à la puissance 0,5). Ici, nous avons arbitrairement choisi une puissance de 1,5. Lorsque nous utiliserons gls avec varPower, nous tenterons de récupérer la « vraie » valeur de 1,5 en plus de 2.
set.seed(4)n <- 1000x <- seq(1,10,length.out = n)y <- 1.2 + 2.1*x + rnorm(n, sd = 2*abs(x)^1.5)d <- data.frame(y, x)ggplot(d, aes(x, y)) + geom_point()
Nous voyons que les données présentent la forme classique d’une variance qui augmente lorsque le prédicteur augmente. Pour modéliser correctement ces données en utilisant gls, nous lui fournissons la formule y ~x et utilisons l’argument weights avec varPower. Remarquez que nous spécifions la formule unilatérale tout comme nous l’avons fait avec la fonction varFixed. Dans ce modèle, cependant, nous obtiendrons une estimation de la puissance.
vm3 <- gls(y ~ x, data = d, weights = varPower(form = ~x))summary(vm3)## Generalized least squares fit by REML## Model: y ~ x ## Data: d ## AIC BIC logLik## 8840.188 8859.811 -4416.094## ## Variance function:## Structure: Power of variance covariate## Formula: ~x ## Parameter estimates:## power ## 1.520915 ## ## Coefficients:## Value Std.Error t-value p-value## (Intercept) 2.321502 0.4750149 4.887218 0## x 1.582272 0.2241367 7.059404 0## ## Correlation: ## (Intr)## x -0.845## ## Standardized residuals:## Min Q1 Med Q3 Max ## -2.892319514 -0.655237190 0.001653178 0.691241690 3.346273816 ## ## Residual standard error: 1.872944 ## Degrees of freedom: 1000 total; 998 residualint <- intervals(vm3)int$varStruct## lower est. upper## power 1.445064 1.520915 1.596765## attr(,"label")## "Variance function:"int$sigma## lower est. upper ## 1.650987 1.872944 2.124740 ## attr(,"label")## "Residual standard error:"La puissance est estimée à 1,52, ce qui est très proche de la « vraie » valeur de 1,5. L’erreur standard résiduelle estimée de 1,8729439 est également assez proche de 2. Les deux intervalles capturent les valeurs que nous avons utilisées pour simuler les données. Les coefficients du modèle, en revanche, sont quelque peu mal estimés. Cela n’est pas surprenant étant donné la quantité de variabilité dans y, en particulier pour \(x > 2\).
varExp
La fonction varExp nous permet de modéliser la variance comme une fonction exponentielle d’une covariable. Encore une fois, nous allons expliquer cette fonction de variance à l’aide de données simulées. Le seul changement concerne l’argument sd de la fonction rnorm. Nous avons une valeur fixe de 2 que nous multiplions par x qui est elle-même multipliée par 0,5 et exponentiée. Remarquez la rapidité avec laquelle la variance augmente lorsque nous augmentons x. Ceci est dû à la croissance exponentielle de la variance.
set.seed(55)n <- 400x <- seq(1, 10, length.out = n)y <- 1.2 + 2.1*x + rnorm(n, sd = 2*exp(0.5*x))d <- data.frame(y, x)ggplot(d, aes(x, y)) + geom_point()
Pour travailler à rebours et récupérer ces valeurs, nous utilisons la fonction varExp dans l’argument weights de gls. La formule unilatérale ne change pas dans ce cas. Elle dit de modéliser la variance comme une fonction de x. La fonction varExp indique que x a été multipliée par une certaine valeur et exponentielle, donc gls va essayer d’estimer cette valeur.
vm4 <- gls(y ~ x, data = d, weights = varExp(form = ~x))summary(vm4)## Generalized least squares fit by REML## Model: y ~ x ## Data: d ## AIC BIC logLik## 3873.153 3889.098 -1932.576## ## Variance function:## Structure: Exponential of variance covariate## Formula: ~x ## Parameter estimates:## expon ## 0.4845623 ## ## Coefficients:## Value Std.Error t-value p-value## (Intercept) 1.198525 1.1172841 1.072713 0.284## x 2.073297 0.4933502 4.202486 0.000## ## Correlation: ## (Intr)## x -0.892## ## Standardized residuals:## Min Q1 Med Q3 Max ## -3.16976627 -0.69525696 0.00614222 0.63718022 2.90423217 ## ## Residual standard error: 2.119192 ## Degrees of freedom: 400 total; 398 residualint <- intervals(vm4)int$varStruct## lower est. upper## expon 0.4562871 0.4845623 0.5128374## attr(,"label")## "Variance function:"int$sigma## lower est. upper ## 1.786737 2.119192 2.513505 ## attr(,"label")## "Residual standard error:"L’estimation « expon » de 0,4845623 dans la section « fonction de variance » est très proche de notre valeur spécifiée de 0,5. De même, l’erreur standard résiduelle estimée de 2,1191918 est proche de la « vraie » valeur de 2. Les estimations des coefficients du modèle sont également proches des valeurs que nous avons utilisées pour générer les données, mais remarquez l’incertitude dans l'(Intercept). Le test d’hypothèse dans le résumé ne peut pas exclure un intercept négatif. Encore une fois, cela n’est pas surprenant puisque y a tellement de variance non constante ainsi que le fait que nous n’avons aucune observation de x égale à 0. Puisque l’intercept est la valeur que y prend lorsque x est égale à 0, notre intercept estimé est une extrapolation à un événement que nous n’avons pas observé.
varConstPower
La fonction varConstPower nous permet de modéliser la variance comme une constante positive plus une puissance de la valeur absolue d’une covariable. C’est une mise en bouche, mais c’est fondamentalement la même chose que la fonction varPower, sauf que maintenant nous lui ajoutons une constante positive. Le code suivant simule des données pour lesquelles la fonction varConstPower serait appropriée à utiliser. Remarquez qu’il est identique au code que nous avons utilisé pour simuler les données de la section varPower, sauf que nous ajoutons 0,7 à x dans la fonction abs. Pourquoi 0,7 ? Aucune raison particulière. C’est juste une constante positive que nous avons choisie.
set.seed(4)n <- 1000x <- seq(1,10,length.out = n)y <- 1.2 + 2.1*x + rnorm(n, sd = 2*abs(0.7 + x)^1.5)d <- data.frame(y, x)ggplot(d, aes(x, y)) + geom_point()
La façon correcte de modéliser ces données est d’utiliser gls avec la fonction varConstPower. La formule unilatérale est la même que précédemment. Le résumé renvoie trois estimations pour le modèle de variance : la constante, la puissance et l’erreur standard résiduelle. Rappelons que les « vraies » valeurs sont 0,7, 1,5 et 2, respectivement.
vm5 <- gls(y ~ x, data = d, weights = varConstPower(form = ~x))summary(vm5)## Generalized least squares fit by REML## Model: y ~ x ## Data: d ## AIC BIC logLik## 9319.794 9344.323 -4654.897## ## Variance function:## Structure: Constant plus power of variance covariate## Formula: ~x ## Parameter estimates:## const power ## 0.3554921 1.3593016 ## ## Coefficients:## Value Std.Error t-value p-value## (Intercept) 2.763788 0.7607449 3.633003 3e-04## x 1.456677 0.2915621 4.996112 0e+00## ## Correlation: ## (Intr)## x -0.824## ## Standardized residuals:## Min Q1 Med Q3 Max ## -2.846999539 -0.664551194 0.006709051 0.679895874 3.287366683 ## ## Residual standard error: 2.887761 ## Degrees of freedom: 1000 total; 998 residualint <- intervals(vm5)int$varStruct## lower est. upper## const 0.02658804 0.3554921 4.753063## power 1.14258677 1.3593016 1.576016## attr(,"label")## "Variance function:"int$sigma## lower est. upper ## 1.814757 2.887761 4.595195 ## attr(,"label")## "Residual standard error:"Les intervalles capturent les « vraies » valeurs, bien que l’intervalle pour la constante soit assez large. Si nous ne connaissions pas la vraie valeur, il semblerait que la constante pourrait plausiblement être 2 ou même 4.
varComb
Enfin, la fonction varComb nous permet de modéliser des combinaisons de deux modèles de variance ou plus en multipliant ensemble les fonctions de variance correspondantes. Évidemment, cela peut accommoder des modèles de variance très complexes. Nous allons simuler quelques données de base qu’il serait approprié d’utiliser avec la fonction varComb.
Ce que nous faisons ci-dessous est de combiner deux processus de variance différents :
- L’un qui permet à l’écart-type de différer selon le sexe (« f » = \(2\), « m » = \(2 \times 2.5\))
- Un qui permet à l’écart-type d’augmenter lorsque
xaugmente, oùxest multiplié par 1,5 et exponentialisé
Pour aider à visualiser les données, nous avons limité x à une plage de 1 à 3.
set.seed(77)n <- 400x <- seq(1, 3, length.out = n)g <- sample(c("m","f"), size = n, replace = TRUE)msd <- ifelse(g=="m", 2*2.5, 2) * exp(1.5*x)y <- 1.2 + 2.1*x - 1.5*(g == "m") + 2.8*x*(g == "m") + rnorm(n, sd = msd)d <- data.frame(y, x, g)ggplot(d, aes(x, y, color = g)) + geom_point()
Le graphique montre une variance croissante à mesure que x augmente, mais aussi des différences de variance entre les sexes. La variance de y pour le groupe « m » est beaucoup plus grande que la variance de y dans le groupe « f », surtout lorsque x est supérieure à 1,5.
Pour modéliser correctement le processus de génération de données que nous avons spécifié ci-dessus et tenter de récupérer les vraies valeurs, nous utilisons la fonction varComb comme enveloppe autour de deux autres fonctions de variance : varIdent et varExp. Pourquoi ces deux fonctions ? Parce que nous avons des variances différentes entre les sexes, et nous avons une variance qui augmente exponentiellement en fonction de x.
vm6 <- gls(y ~ x*g, data = d, weights = varComb(varIdent(form = ~ 1|g), varExp(form = ~x)))summary(vm6)## Generalized least squares fit by REML## Model: y ~ x * g ## Data: d ## AIC BIC logLik## 4431.45 4459.32 -2208.725## ## Combination of variance functions: ## Structure: Different standard deviations per stratum## Formula: ~1 | g ## Parameter estimates:## f m ## 1.000000 2.437046 ## Structure: Exponential of variance covariate## Formula: ~x ## Parameter estimates:## expon ## 1.540608 ## ## Coefficients:## Value Std.Error t-value p-value## (Intercept) -5.996337 6.539552 -0.9169339 0.3597## x 8.829971 4.849272 1.8208860 0.0694## gm -17.572238 16.548540 -1.0618603 0.2889## x:gm 16.932938 12.192793 1.3887661 0.1657## ## Correlation: ## (Intr) x gm ## x -0.972 ## gm -0.395 0.384 ## x:gm 0.387 -0.398 -0.974## ## Standardized residuals:## Min Q1 Med Q3 Max ## -2.74441479 -0.67478610 -0.00262221 0.66079254 3.14975288 ## ## Residual standard error: 1.794185 ## Degrees of freedom: 400 total; 396 residualLa section de résumé contient deux sections pour modéliser la variance. La première estime le multiplicateur du groupe « m » à 2,437, ce qui est très proche de la vraie valeur de 2,5. Le paramètre exponentiel est estimé à 1,54, ce qui est extrêmement proche de la valeur réelle de 1,5 que nous avons utilisée lors de la génération des données. Enfin, l’erreur standard résiduelle est estimée à environ 1,79, proche de la vraie valeur de 2. L’appel intervals(vm6) montre des intervalles de confiance très serrés. Les préfixes A et B dans l’élément varStruct sont juste des étiquettes pour les deux modèles de variance différents.
int <- intervals(vm6)int$varStruct## lower est. upper## A.m 2.119874 2.437046 2.801673## B.expon 1.419366 1.540608 1.661850## attr(,"label")## "Variance function:"int$sigma## lower est. upper ## 1.380008 1.794185 2.332669 ## attr(,"label")## "Residual standard error:"Malheureusement, en raison de la grande variabilité exponentielle, les estimations des coefficients du modèle sont terriblement mauvaises.
À quoi bon ?
Alors pourquoi avons-nous simulé toutes ces fausses données pour ensuite tenter de récupérer les « vraies » valeurs ? A quoi cela sert-il ? Questions justes. La réponse est que cela nous aide à comprendre quelles hypothèses nous faisons lorsque nous spécifions et ajustons un modèle. Lorsque nous spécifions et ajustons le modèle suivant…
m <- lm(y ~ x1 + x2)… nous disons que nous pensons que y est approximativement égal à une somme pondérée de x1 et x2 (plus un intercept), les erreurs étant des tirages aléatoires d’une distribution normale avec une moyenne de 0 et un écart type fixe. De même, si nous spécifions et ajustons le modèle suivant…
m <- gls(y ~ x1 + x2, data = d, weights = varPower(form = ~x1))… nous disons que nous pensons que y est approximativement égal à une somme pondérée de x1 et x2 (plus une ordonnée à l’origine), les erreurs étant des tirages aléatoires d’une distribution normale avec une moyenne de 0 et un écart type qui croît comme un multiple de x1 augmenté d’une certaine puissance.
Si nous pouvons simuler des données adaptées à ces modèles, alors nous avons une meilleure compréhension des hypothèses que nous faisons lorsque nous utilisons ces modèles. J’espère que vous avez maintenant une meilleure compréhension de ce que vous pouvez faire pour modéliser la variance en utilisant le package nlme.
Pour toute question ou clarification concernant cet article, contactez le StatLab de la bibliothèque de l’UVA : [email protected]
Voir la collection complète des articles du StatLab de la bibliothèque de l’UVA.
Clay Ford
Consultant en recherche statistique
Bibliothèque de l’Université de Virginie
7 avril 2020
-
1. Le package
nlmeest peut-être mieux connu pour sa fonctionlme, qui est utilisée pour ajuster des modèles à effets mixtes (c’est-à-dire des modèles avec des effets fixes et aléatoires). Ce billet de blog démontre les fonctions de variance en utilisantgls, qui n’ajuste pas les effets aléatoires. Cependant, tout ce que nous présentons dans ce billet de blog peut être utilisé aveclme.