RegLinSimple
E.Ing NEKUI TIEFANG BIEDIANT
2023-08-25
library(knitr)
library(rmarkdown)
library(markdown) X = c(20,24,28,22,32,28,32,36,41,41)
Y = c(16,18,23,24,28,29,26,31,32,34)
df=data.frame(X,Y)
dfplot(X,Y, col="blue", xlab = "Quantité d'engrais en kg", ylab = "rendement de mais en quintal", main = " nuage des points (Xi, Yi) ")
model =lm(Y~X)
coef(model)## (Intercept) X
## 4.3927701 0.7140536summary(model)##
## Call:
## lm(formula = Y ~ X)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.5301 -1.5983 -0.4557 0.8654 4.6137
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.3928 3.9718 1.106 0.300874
## X 0.7141 0.1273 5.609 0.000505 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.825 on 8 degrees of freedom
## Multiple R-squared: 0.7973, Adjusted R-squared: 0.7719
## F-statistic: 31.46 on 1 and 8 DF, p-value: 0.0005049plot(X,Y, col="blue", xlab = "Quantité d'engrais en kg", ylab = "rendement de mais en quintal", main = " nuage des points (Xi, Yi) ")
abline(model, col = "red")
l’equation de la droite ajustée est: \[Ŷ_{i}=4.3927701 +0.7140536 *X_{i}\]
Ŷ=fitted(model)
resid=residuals.lm(model)
df2=data.frame(X,Y,Ŷ,resid)
df2validité du model ce model est valide si les conditions suivantes sont respectées: 1-linéarité 2-indépendance des residus (test de Durbin-watson pour la correlation des residus) 3-Homoscedasticité 4-normalité des residus (shapiro-wilk) 5-absence de valeurs aberrantes (toute les valeurs doivent êtres représentes)
5-1) linearité
library(car)
scatterplot(Y~X,data = df, regLine=list(method=lm, lwd=1, col="red") )
Un lissage permet de se donner une idée précise de la relation
existante entre deux variables et d’estimer la linéarité de cette
relation ! existance des valeurs peux ou pas expliqué par le
medel
ou encore
library(lmtest)
raintest(model)##
## Rainbow test
##
## data: model
## Rain = 2.6149, df1 = 5, df2 = 3, p-value = 0.2293p-value sup a 5% linearité accepté
5-2)indépendance des residus
library(performance)
library(lmtest)
check_autocorrelation(model) ## OK: Residuals appear to be independent and not autocorrelated (p = 0.334).dwtest(model)##
## Durbin-Watson test
##
## data: model
## DW = 1.6837, p-value = 0.1808
## alternative hypothesis: true autocorrelation is greater than 0il faut une p-value supérieure à 0,05 pour avoir indépendance.
conclusion: residus indépendance !!!
Attention : ce test n’évalue que l’autocorrélation entre résidus et le résidus qui suit directement (n+1) et non ceux séparés d’un intervalle (résidus+2) qui pourraient être autocorrélés dans un phénomène cyclique.
*On peut aussi visualiser l’indépendance avec ce graphique suivant :
library(stats)
acf(residuals(model), main="Regression Y = f(X)")
L’interprétation de ce graphique se fait de la manière suivante : -Le premier bâtonnet est très élevé, c’est l’auto-corrélation des résidus avec eux-même ! -Le deuxième bâtonnet indique l’auto-corrélation entre les résidus et les résidus n+1 : il y a auto-corrélation dès que le bâtonnet (lag) dépasse les pointillés. -Le troisième bâtonnet entre les résidus n et les résidus n+2… etc.
5-3)homoscedasticité
library(car)
plot(model, which = 3)
On cherche ici une courbe rouge plane. L’homogénéité est à rejeter si celle-ci n’est pas horizontale.
*homoscedas rejeté
library("lmtest") ;
bptest(model) # librairie lmtest##
## studentized Breusch-Pagan test
##
## data: model
## BP = 3.6436, df = 1, p-value = 0.05628ncvTest(model) # librairie car## Non-constant Variance Score Test
## Variance formula: ~ fitted.values
## Chisquare = 2.224827, Df = 1, p = 0.13581L’homogénéité est rejeté si la p-value est inférieure à 0,05
difficile de conclure !!!gqtest(model)##
## Goldfeld-Quandt test
##
## data: model
## GQ = 0.39808, df1 = 3, df2 = 3, p-value = 0.7654
## alternative hypothesis: variance increases from segment 1 to 2Vérification par le test de Goldfeld-Quandt (homogénéité : p-value > 0,05)
library(performance)
check_heteroscedasticity(model)## OK: Error variance appears to be homoscedastic (p = 0.136).Homoscedasticité des residus verifié !!!
5-4)normalité
hist(residuals.lm(model),col="yellow",freq=F)
densite <- density(residuals.lm(model)) # estimer la densité que représente ces différentes valeurs
lines(densite) # Superposer une ligne de densité à l'histogramme
autre methodes recommandé
shapiro.test(residuals.lm(model))##
## Shapiro-Wilk normality test
##
## data: residuals.lm(model)
## W = 0.9315, p-value = 0.4628check_normality(model)## OK: residuals appear as normally distributed (p = 0.409).Normalité des residus Validé…!!!
5-5) absence de valeurs aberrantes (toute les valeurs doivent êtres présentes)
plot(model, which = 1) 
plot(model, which = 2) 
plot(model, which = 3)
plot(model, which = 4) 
plot(model, which = 5) # ce qui m'interesse !!!!
plot(model, which = 6) 
le lien de la documentation de tout les cas: https://data.library.virginia.edu/diagnostic-plots/
aucune valeur en dehors de la cook’s distance Donc pas de valeur extrême
influenceIndexPlot(model) 
pas reelement de données extrême essayons kamem d’eliminer quelque valeur qui s’y raproche et comparons les deux modèles
model2 <- lm(Y~X, data=df[-c(2,4),])
compareCoefs(model ,model2) ## Calls:
## 1: lm(formula = Y ~ X)
## 2: lm(formula = Y ~ X, data = df[-c(2, 4), ])
##
## Model 1 Model 2
## (Intercept) 4.39 3.36
## SE 3.97 4.28
##
## X 0.714 0.745
## SE 0.127 0.130
## boxplot(df2$resid)
boxplot(X)
boxplot(Y)
limites de turkey bien aux limites de notre cadrant Donc pas de valeur
extrême quelque soit la variable
CONCLUSION NOTRE MODELE \[Ŷ_{i}=4.3927701 +0.7140536 *X_{i}\] PEUT ETRE CONSIDERE COMME VALIDE PUISQUE TOUTES LES CONDITIONS SONT VALIDES
library(KefiR)
valreg## function (reg, verbose = TRUE, nvar = 5, boot = TRUE, alpha = 0.05,
## conf.level = 0.95, plot = TRUE, data = c(), raintest_alpha = 0.05,
## dwtest_alpha = 0.03, shapiro_alpha = 0.05, bptest_alpha = 0.05)
## {
## if (plot == T) {
## boot <- TRUE
## }
## error <- TRUE
## nvar1 <- length(coef(reg))
## if (nvar1 > nvar) {
## if (verbose == TRUE) {
## cat("More than ", nvar, " variables\n")
## }
## error = FALSE
## }
## else if (length(reg$fitted.values) < (nvar + 1)) {
## if (verbose == TRUE) {
## cat("Not enough values in the subset\n")
## }
## error = FALSE
## }
## else if ((length(summary(reg)[[4]][, 4])) != (length(coef(reg)))) {
## error = FALSE
## }
## else {
## if (verbose == TRUE) {
## cat("01- Analysis of the p-values of the model and its coefficients.\n")
## }
## pval_mdl <- pf(summary(reg)$fstatistic[1], summary(reg)$fstatistic[2],
## summary(reg)$fstatistic[3], lower.tail = FALSE)
## if (pval_mdl > alpha) {
## if (verbose == TRUE) {
## cat("\tWarning!\n\tBad significance of the model. p-value:",
## pval_mdl, "\n")
## }
## error = FALSE
## }
## else {
## if (verbose == TRUE) {
## cat("\tGood significance of the model. p-value:",
## pval_mdl, "\n")
## }
## }
## pval_coeff <- summary(reg)[[4]][, 4]
## if (max(pval_coeff) > alpha) {
## if (verbose == TRUE) {
## cat("\tWarning!\n\tBad significance of the coefficients. max(p.value) :",
## max(pval_coeff), "\n")
## }
## error = FALSE
## }
## else {
## if (verbose == TRUE) {
## cat("\tGood significance of the coefficients. max(pval_coeff) :",
## max(pval_coeff), "\n")
## }
## }
## if (verbose == TRUE) {
## cat("02- Analysis of the adequacy of model (Equivalence between the global model and the model established on the best points.).\n")
## }
## pvalt <- raintest(reg, order.by = "mahalanobis")$p.value
## if (pvalt < raintest_alpha) {
## if (verbose == TRUE) {
## cat("\tWarning!\n\tRainbow test ordered by mahalanobis (raintest()) - Bad adequacy. p.value : ",
## pvalt, "\n")
## }
## error = FALSE
## }
## else {
## if (verbose == TRUE) {
## cat("\tRainbow test (raintest()) - Good adequacy. p.value : ",
## pvalt, "\n")
## }
## }
## if (verbose == TRUE) {
## cat("03- Analysis of independence of the residuals.\n")
## }
## if (length(data) > 0) {
## ypred <- predict(reg, data)
## }
## else {
## ypred <- predict(reg, reg$model)
## }
## pvalt <- try(dwtest(reg, ypred)$p.value)
## if (pvalt < dwtest_alpha) {
## if (verbose == TRUE) {
## cat("\tWarning!\n\tDurbin-Watson test (dwtest()) - Bad independence of the residuals. p.value : ",
## pvalt, "\n")
## }
## error = FALSE
## }
## else {
## if (verbose == TRUE) {
## cat("\tDurbin-Watson test (dwtest()) - Good independence of the residuals. p.value : ",
## pvalt, "\n")
## }
## }
## if (verbose == TRUE) {
## cat("04- Analysis of distribution of residuals.\n")
## }
## pvalt <- shapiro.test(residuals(reg))$p.value
## if (pvalt < shapiro_alpha) {
## if (verbose == TRUE) {
## cat("\tWarning!\n\tShapiro-Wilk test (shapiro.test()) - Non-normal distribution of residuals. p.value : ",
## pvalt, "\n")
## }
## error = FALSE
## }
## else {
## if (verbose == TRUE) {
## cat("\tShapiro-Wilk test (shapiro.test()) - Normal distribution of residuals. p.value : ",
## pvalt, "\n")
## }
## }
## if (verbose == TRUE) {
## cat("05- Analysis of variance of residuals.\n")
## }
## if (length(coef(reg)) >= 2) {
## pvalt <- bptest(reg)$p.value
## if (pvalt < bptest_alpha) {
## if (verbose == TRUE) {
## cat("\tWarning!\n\tBreush-Pagan test (bptest()) - Non-constant variance of the residuals. p.value : ",
## pvalt, "\n")
## }
## error = FALSE
## }
## else {
## if (verbose == TRUE) {
## cat("\tBreush-Pagan test (bptest()) - Constant variance of the residuals. p.value : ",
## pvalt, "\n")
## }
## }
## }
## cooksd <- cooks.distance(reg)
## if (verbose == TRUE) {
## cat("06- Analysis of leverage effect.\n")
## }
## if (max(cooksd, na.rm = TRUE) > 1) {
## if (verbose == TRUE) {
## cat("\tWarning!\n\tCook's distance (cooks.distance()) - Leverage effect. max(cooks.distance())",
## max(cooksd, na.rm = TRUE), "\n")
## }
## error = FALSE
## }
## else {
## if (verbose == TRUE) {
## cat("\tCook's distance (cooks.distance()) - No leverage effect. max(cooks.distance())",
## max(cooksd, na.rm = TRUE), "\n")
## }
## }
## if (boot == TRUE) {
## if (verbose == TRUE) {
## cat("07- Analysis of solidity of model by boostrap.\n")
## }
## if (length(data) > 0) {
## bootres <- bootreg(reg, verbose = FALSE, alpha = alpha,
## conf.level = conf.level, plot = plot, data = data)
## }
## else {
## bootres <- bootreg(reg, verbose = FALSE, alpha = alpha,
## conf.level = conf.level, plot = plot)
## }
## if (bootres == FALSE) {
## if (verbose == TRUE) {
## cat("\tWarning!\n\tBootstrap (bootreg()) - Fragility of the model in boostrap. Please, use bootreg()\n")
## }
## error = FALSE
## }
## else {
## if (verbose == TRUE) {
## cat("\tBootstrap (bootreg()) - Solidity of the model in boostrap.\n")
## }
## }
## }
## }
## return(error)
## }
## <bytecode: 0x562eea413f30>
## <environment: namespace:KefiR>valreg(model, plot = TRUE )## 01- Analysis of the p-values of the model and its coefficients.
## Good significance of the model. p-value: 0.0005048663
## Warning!
## Bad significance of the coefficients. max(p.value) : 0.3008742
## 02- Analysis of the adequacy of model (Equivalence between the global model and the model established on the best points.).
## Rainbow test (raintest()) - Good adequacy. p.value : 0.5017592
## 03- Analysis of independence of the residuals.
## Durbin-Watson test (dwtest()) - Good independence of the residuals. p.value : 0.7984325
## 04- Analysis of distribution of residuals.
## Shapiro-Wilk test (shapiro.test()) - Normal distribution of residuals. p.value : 0.4628321
## 05- Analysis of variance of residuals.
## Breush-Pagan test (bptest()) - Constant variance of the residuals. p.value : 0.05628425
## 06- Analysis of leverage effect.
## Cook's distance (cooks.distance()) - No leverage effect. max(cooks.distance()) 0.4045383
## 07- Analysis of solidity of model by boostrap.
## Warning!
## Bootstrap (bootreg()) - Fragility of the model in boostrap. Please, use bootreg()## [1] FALSEbootreg(model)
_____FIN____