On charge le jeu de données
library(carData)
head(Salaries)
## rank discipline yrs.since.phd yrs.service sex salary
## 1 Prof B 19 18 Male 139750
## 2 Prof B 20 16 Male 173200
## 3 AsstProf B 4 3 Male 79750
## 4 Prof B 45 39 Male 115000
## 5 Prof B 40 41 Male 141500
## 6 AssocProf B 6 6 Male 97000
Lien entre 2 variables categorielles : chi-deux
#### Y a-t-il un lien entre le statut ('rank') et le genre ('sex') ?
table(Salaries$rank) # description de la variable 'rank'
##
## AsstProf AssocProf Prof
## 67 64 266
table(Salaries$sex) # description de la variable 'sex'
##
## Female Male
## 39 358
tableau <- table(Salaries$sex, Salaries$rank) # tableau croisé des 2 variables
tableau
##
## AsstProf AssocProf Prof
## Female 11 10 18
## Male 56 54 248
addmargins(tableau) # on fait apparaitre les effectifs marginaux
##
## AsstProf AssocProf Prof Sum
## Female 11 10 18 39
## Male 56 54 248 358
## Sum 67 64 266 397
prop.table(tableau, 1)*100 # pourcentages suivant les lignes (ex: 25.64% des femmes sont AssocProf)
##
## AsstProf AssocProf Prof
## Female 28.20513 25.64103 46.15385
## Male 15.64246 15.08380 69.27374
prop.table(tableau, 2)*100 # pourcentages suivant les colonnes (ex: 15.63% des AssocProf sont des femmes)
##
## AsstProf AssocProf Prof
## Female 16.417910 15.625000 6.766917
## Male 83.582090 84.375000 93.233083
chisq.test(tableau) # test du chi-deux pour voir si le lien est significatif
##
## Pearson's Chi-squared test
##
## data: tableau
## X-squared = 8.5259, df = 2, p-value = 0.01408
Test sur une différence de moyennes : t de Student
#### Y a-t-il une différence de salaire entre hommes et femmes ?
# une commande pour afficher les moyennes : aggregate
aggregate(Salaries$salary, by = list(Salaries$sex), FUN = mean)
## Group.1 x
## 1 Female 101002.4
## 2 Male 115090.4
# on crée un nouveau tableau avec seulement les hommes
hommes <- Salaries[Salaries$sex == "Male",]
# on crée un nouveau tableau avec seulement les femmes
femmes <- Salaries[Salaries$sex == "Female",]
# commande t de Student
t.test(hommes$salary, femmes$salary, paired=FALSE, alternative = "greater")
##
## Welch Two Sample t-test
##
## data: hommes$salary and femmes$salary
## t = 3.1615, df = 50.122, p-value = 0.001332
## alternative hypothesis: true difference in means is greater than 0
## 95 percent confidence interval:
## 6620.263 Inf
## sample estimates:
## mean of x mean of y
## 115090.4 101002.4
# paired=FALSE : on effectue un test sur deux échantillons indépendants (hommes et femmes)
# alternative = "greater" : le test est unilatéral -> l'hypothèse est que les hommes gagnent plus que les femmes
# -> l'effet du genre est significatif (p = .001) : les hommes (115090.4 $/an) gagnent plus que les femmes (101002.4 $/an)
Analyse multivariée : régression linéaire
#### Cas 1: VD = salaire, une VI = nombre d'années de service
# -> Peut-on prédire le salaire à partir du nombre d'années de service ?
library(psych) # on charge le package 'psych'
describe(Salaries$salary) # description de la variable 'salary'
## vars n mean sd median trimmed mad min max range skew
## X1 1 397 113706.5 30289.04 107300 111401.6 29355.48 57800 231545 173745 0.71
## kurtosis se
## X1 0.18 1520.16
hist(Salaries$salary,100) # histogramme
describe(Salaries$yrs.service) # description de la variable 'yrs.service'
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 397 17.61 13.01 16 16.51 14.83 0 60 60 0.65 -0.34 0.65
hist(Salaries$yrs.service,100) # histogramme
model <- lm(salary ~ yrs.service, data=Salaries) # commande de la regression
summary(model)
##
## Call:
## lm(formula = salary ~ yrs.service, data = Salaries)
##
## Residuals:
## Min 1Q Median 3Q Max
## -81933 -20511 -3776 16417 101947
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 99974.7 2416.6 41.37 < 2e-16 ***
## yrs.service 779.6 110.4 7.06 7.53e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 28580 on 395 degrees of freedom
## Multiple R-squared: 0.1121, Adjusted R-squared: 0.1098
## F-statistic: 49.85 on 1 and 395 DF, p-value: 7.529e-12
# -> l'effet du nombre d'années de service est significatif (p < 0.001) :
# une année de service en plus se traduit par 779.6 dollars par an en plus
# pour 0 année de service, le salaire annuel est de 99974.7 (Intercept)
# graphique
plot(Salaries$yrs.service, Salaries$salary, xlab="Années de service", ylab="Salaire ($/an)", xlim=c(0,60), ylim=c(55000,250000))
abline(lm(salary ~ yrs.service, data=Salaries), col="red")
#### Cas 2: VD = salaire, deux VI = nombre d'années de service et genre
# -> Peut-on prédire le salaire à partir du nombre d'années de service ET du genre ?
model1 <- lm(salary ~ sex, data=Salaries) # VI: genre
summary(model1)
##
## Call:
## lm(formula = salary ~ sex, data = Salaries)
##
## Residuals:
## Min 1Q Median 3Q Max
## -57290 -23502 -6828 19710 116455
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 101002 4809 21.001 < 2e-16 ***
## sexMale 14088 5065 2.782 0.00567 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 30030 on 395 degrees of freedom
## Multiple R-squared: 0.01921, Adjusted R-squared: 0.01673
## F-statistic: 7.738 on 1 and 395 DF, p-value: 0.005667
# -> l'effet du genre est significatif (même résultat que le t de Student). Le coefficient de régression correspond à la différence de salaire entre hommes et femmes
model2 <- lm(salary ~ yrs.service + sex, data=Salaries) # on met les 2 VI
summary(model2)
##
## Call:
## lm(formula = salary ~ yrs.service + sex, data = Salaries)
##
## Residuals:
## Min 1Q Median 3Q Max
## -81757 -20614 -3376 16779 101707
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 92356.9 4740.2 19.484 < 2e-16 ***
## yrs.service 747.6 111.4 6.711 6.74e-11 ***
## sexMale 9071.8 4861.6 1.866 0.0628 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 28490 on 394 degrees of freedom
## Multiple R-squared: 0.1198, Adjusted R-squared: 0.1154
## F-statistic: 26.82 on 2 and 394 DF, p-value: 1.201e-11
# -> l'effet du nombre d'années de service est significatif (p < 0.001) : une année de service en plus se traduit par 747.6 dollars par an en plus
# -> l'effet du genre n'est plus significatif (p = 0.06)
# -> car l'effet du genre est confondu avec l'effet du nombre d'années de service
# on peut le vérifier avec un t de Student
t.test(hommes$yrs.service, femmes$yrs.service, paired=FALSE, alternative = "greater")
##
## Welch Two Sample t-test
##
## data: hommes$yrs.service and femmes$yrs.service
## t = 4.2604, df = 58.559, p-value = 3.752e-05
## alternative hypothesis: true difference in means is greater than 0
## 95 percent confidence interval:
## 4.077543 Inf
## sample estimates:
## mean of x mean of y
## 18.27374 11.56410
# -> les hommes (18.27 ans) ont plus d'années de service que les femmes (11.56 ans)
# -> c'est cette différence qui explique la différence de salaire entre hommes et femmes
Régression linéaire : graphique
#### on effectue la regression séparément pour hommes et femmes
by(Salaries, Salaries$sex, function(x) summary(lm(salary ~ yrs.service, data=x)))
## Salaries$sex: Female
##
## Call:
## lm(formula = salary ~ yrs.service, data = x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -55205 -10417 3534 13530 54473
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 82068.5 5822.7 14.095 < 2e-16 ***
## yrs.service 1637.3 402.4 4.069 0.000238 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 21860 on 37 degrees of freedom
## Multiple R-squared: 0.3092, Adjusted R-squared: 0.2905
## F-statistic: 16.56 on 1 and 37 DF, p-value: 0.0002378
##
## ------------------------------------------------------------
## Salaries$sex: Male
##
## Call:
## lm(formula = salary ~ yrs.service, data = x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -80381 -20665 -4965 16469 102536
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 102197.1 2617.5 39.043 < 2e-16 ***
## yrs.service 705.6 116.1 6.078 3.14e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 29010 on 356 degrees of freedom
## Multiple R-squared: 0.094, Adjusted R-squared: 0.09146
## F-statistic: 36.94 on 1 and 356 DF, p-value: 3.14e-09
# -> l'effet du nombre d'années de service est plus fort chez les femmes que chez les hommes
# graphique
symb <- c(21,22) # on choisit les symboles
colors <- c("red","blue") # on choisit les couleurs
plot(Salaries$yrs.service, Salaries$salary, pch=symb[Salaries$sex], col=colors[Salaries$sex], bg=colors[Salaries$sex], cex=1.0, xlab="Années de service", ylab="Salaire ($/an)", xlim=c(0,60), ylim=c(55000,250000))
abline(lm(salary ~ yrs.service, data=Salaries[Salaries$sex == "Male",]), col="blue")
abline(lm(salary ~ yrs.service, data=Salaries[Salaries$sex == "Female",]), col="red")
legend("topleft",
legend = c("Hommes","Femmes"),
col = c("blue","red"),
pch = c(21,22),
bty = "n",
pt.cex = 1,
cex = 0.8,
text.col = "black",
horiz = F ,
inset = c(0.01, 0.01))
#### on effectue la regression séparément pour les trois statuts (Professeur, Professeur associé, Professeur assistant)
by(Salaries, Salaries$rank, function(x) summary(lm(salary ~ yrs.service, data=x)))
## Salaries$rank: AsstProf
##
## Call:
## lm(formula = salary ~ yrs.service, data = x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -20502 -6639 -1009 8372 14209
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 78926.6 1878.0 42.028 <2e-16 ***
## yrs.service 779.3 670.9 1.162 0.25
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8153 on 65 degrees of freedom
## Multiple R-squared: 0.02034, Adjusted R-squared: 0.005265
## F-statistic: 1.349 on 1 and 65 DF, p-value: 0.2496
##
## ------------------------------------------------------------
## Salaries$rank: AssocProf
##
## Call:
## lm(formula = salary ~ yrs.service, data = x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -27027 -12140 2855 9309 32178
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 98594.3 2597.7 37.95 <2e-16 ***
## yrs.service -394.7 166.5 -2.37 0.0209 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 13350 on 62 degrees of freedom
## Multiple R-squared: 0.08307, Adjusted R-squared: 0.06828
## F-statistic: 5.617 on 1 and 62 DF, p-value: 0.02092
##
## ------------------------------------------------------------
## Salaries$rank: Prof
##
## Call:
## lm(formula = salary ~ yrs.service, data = x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -65609 -21909 -3290 18035 106585
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 129494.4 3760.5 34.435 <2e-16 ***
## yrs.service -119.3 147.0 -0.812 0.418
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 27740 on 264 degrees of freedom
## Multiple R-squared: 0.002489, Adjusted R-squared: -0.001289
## F-statistic: 0.6588 on 1 and 264 DF, p-value: 0.4177
# graphique
symb <- c(21,22,23) # on choisit les symboles
colors <- c("red","blue","green") # on choisit les couleurs
plot(Salaries$yrs.service, Salaries$salary, pch=symb[Salaries$rank], col=colors[Salaries$rank], bg=colors[Salaries$rank], cex=1.0, xlab="Années de service", ylab="Salaire ($/an)", xlim=c(0,60), ylim=c(55000,250000))
abline(lm(salary ~ yrs.service, data=Salaries[Salaries$rank == "AsstProf",]), col="red")
abline(lm(salary ~ yrs.service, data=Salaries[Salaries$rank == "AssocProf",]), col="blue")
abline(lm(salary ~ yrs.service, data=Salaries[Salaries$rank == "Prof",]), col="green")
legend("topleft",
legend = c("AsstProf","AssocProf","Prof"),
col = c("red","blue","green"),
pch = c(21,22,23),
bty = "n",
pt.cex = 1,
cex = 0.8,
text.col = "black",
horiz = F ,
inset = c(0.01, 0.01))