Description

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))