Anova 3-czynnikowa

##Wprowadzenie: Przyjęty przeze mnie poziom istotności to alfa=0,01. Przeprowadzamy badanie wpływu zmiennych rank,discipline i sex na wielkość pensji czyli zmiennej salary. Dodatkowo potem przyjrzymy się zmiennej lata od doktoratu(yrs.since.phd) i staż pracy (yrs.service)

##Sprawdzenie czy grupy są równoliczne

table(Salaries$rank, Salaries$discipline)

##            
##               A   B
##   AsstProf   24  43
##   AssocProf  26  38
##   Prof      131 135

table(Salaries$rank, Salaries$sex)

##            
##             Female Male
##   AsstProf      11   56
##   AssocProf     10   54
##   Prof          18  248

table(Salaries$discipline, Salaries$sex)

##    
##     Female Male
##   A     18  163
##   B     21  195

Odp: Jak widać, grupy nie są równoliczne, kobiet jest dużo mniej, jedynie równolicznymi można nazwać dyscypliny A i B

Wykresy pudełkowe

pudelkowy<-ggboxplot(Salaries, x="rank",
          y="salary", color="discipline", facet.by = "sex")
pudelkowy

##Statystyki opisowe

raport <- list(
  "Pensja" =
    list(
      "Min" = ~ min(salary),
      "Max" = ~ max(salary),
      "Q1"  = ~ quantile(salary, 0.25),
      "Mediana" = ~ round(median(salary,0.5),2),
      "Q3"  = ~ quantile(salary, 0.75),
      "Średnia" = ~ round(mean(salary),2),
      "Odch. std." = ~ round(sd(salary),2)
    )  
)
tabela<-summary_table(Salaries, summaries=raport, by=c("sex", "rank", "discipline"))
kbl(tabela, digits=2, caption="Statystyki opisowe") %>%
  kable_classic(full_width=F) %>%
  kable_styling(bootstrap_options=c("striped","hover"))

Statystyki opisowe
	Female.AsstProf.A (N = 6)	Male.AsstProf.A (N = 18)	Female.AssocProf.A (N = 4)	Male.AssocProf.A (N = 22)	Female.Prof.A (N = 8)	Male.Prof.A (N = 123)	Female.AsstProf.B (N = 5)	Male.AsstProf.B (N = 38)	Female.AssocProf.B (N = 6)	Male.AssocProf.B (N = 32)	Female.Prof.B (N = 10)	Male.Prof.B (N = 125)
Min	63100	63900	62884	70000	90450	57800	74692	68404	71065	83900	105450	67559
Max	78500	85000	77500	108413	137000	205500	97032	95079	109650	126431	161101	231545
Q1	72500	72625	70696	81339	101500	102290	77000	79829	103647	95567	123298	114500
Mediana	73000	74000	74065	82350	109800	113278	80225	85408	103872	100730	128256	132825
Q3	76500	75711	75498	88174	116726	136580	92000	90736	104405	105798	143512	150376
Średnia	72933	74270	72129	85049	109632	120619	84190	84647	99436	101622	131836	133518
Odch. std.	5463	4580	6403	10612	15095	28505	9792	6900	14086	9608	17504	26514

Założenia

Obserwacje odstające

Salaries %>% 
  group_by(sex, rank, discipline) %>%
  identify_outliers(salary)

## # A tibble: 18 x 8
##    rank  discipline sex   yrs.since.phd yrs.service salary is.outlier is.extreme
##    <fct> <fct>      <fct>         <int>       <int>  <int> <lgl>      <lgl>     
##  1 Asst~ A          Fema~             7           6  63100 TRUE       FALSE     
##  2 Asso~ A          Fema~            25          22  62884 TRUE       FALSE     
##  3 Asso~ B          Fema~            14           7 109650 TRUE       TRUE      
##  4 Asso~ B          Fema~            12           9  71065 TRUE       TRUE      
##  5 Asst~ A          Male              3           1  63900 TRUE       FALSE     
##  6 Asst~ A          Male              2           0  85000 TRUE       TRUE      
##  7 Asst~ A          Male              8           4  81035 TRUE       FALSE     
##  8 Asso~ A          Male             14           8 100102 TRUE       FALSE     
##  9 Asso~ A          Male              9           7  70000 TRUE       FALSE     
## 10 Asso~ A          Male             11           1 104800 TRUE       FALSE     
## 11 Asso~ A          Male             45          39  70700 TRUE       FALSE     
## 12 Asso~ A          Male             10           1 108413 TRUE       FALSE     
## 13 Asso~ A          Male             11           8 104121 TRUE       FALSE     
## 14 Asso~ B          Male             13          11 126431 TRUE       FALSE     
## 15 Prof  A          Male             29           7 204000 TRUE       FALSE     
## 16 Prof  A          Male             42          18 194800 TRUE       FALSE     
## 17 Prof  A          Male             43          43 205500 TRUE       FALSE     
## 18 Prof  B          Male             38          38 231545 TRUE       FALSE

Widzimy 3 mocno odstające obserwacje.

Normalność

Salaries %>% 
  group_by(sex, rank, discipline) %>%
  shapiro_test(salary)

## # A tibble: 12 x 6
##    rank      discipline sex    variable statistic        p
##    <fct>     <fct>      <fct>  <chr>        <dbl>    <dbl>
##  1 AsstProf  A          Female salary       0.870 0.226   
##  2 AsstProf  B          Female salary       0.889 0.354   
##  3 AssocProf A          Female salary       0.863 0.269   
##  4 AssocProf B          Female salary       0.635 0.00117 
##  5 Prof      A          Female salary       0.934 0.549   
##  6 Prof      B          Female salary       0.974 0.923   
##  7 AsstProf  A          Male   salary       0.941 0.300   
##  8 AsstProf  B          Male   salary       0.941 0.0458  
##  9 AssocProf A          Male   salary       0.878 0.0113  
## 10 AssocProf B          Male   salary       0.967 0.416   
## 11 Prof      A          Male   salary       0.952 0.000259
## 12 Prof      B          Male   salary       0.978 0.0435

ggqqplot(Salaries, "salary") +
  facet_grid(sex~rank~discipline)

Niestety niektóre grupy nie posiadają rozkładu normalnego.

Jednorodność wariancji

wyniki<-Salaries %>%
  levene_test(salary~sex*rank*discipline)
wyniki

## # A tibble: 1 x 4
##     df1   df2 statistic        p
##   <int> <int>     <dbl>    <dbl>
## 1    11   385      9.05 2.06e-14

Wariancje nie są jednorodne

Z racji tego, że nie zostały spełnione założenia testu ANOVA(populacje nie są równoliczne, niektóre grupy nie mają rozkładu normalnego, wariancje nie są jednorodne). Możemy odrzucić obserwacje odstające i zlogarytmować naszą zmienną salary. ### Odrzucenie wartości odstających

library(stats)
library(extraoperators)

## Warning: pakiet 'extraoperators' został zbudowany w wersji R 4.1.2

obserwacje.odstajace<- Salaries %>%
  group_by(sex, rank, discipline) %>%
  identify_outliers(salary)
obserwacje.odstajace

## # A tibble: 18 x 8
##    rank  discipline sex   yrs.since.phd yrs.service salary is.outlier is.extreme
##    <fct> <fct>      <fct>         <int>       <int>  <int> <lgl>      <lgl>     
##  1 Asst~ A          Fema~             7           6  63100 TRUE       FALSE     
##  2 Asso~ A          Fema~            25          22  62884 TRUE       FALSE     
##  3 Asso~ B          Fema~            14           7 109650 TRUE       TRUE      
##  4 Asso~ B          Fema~            12           9  71065 TRUE       TRUE      
##  5 Asst~ A          Male              3           1  63900 TRUE       FALSE     
##  6 Asst~ A          Male              2           0  85000 TRUE       TRUE      
##  7 Asst~ A          Male              8           4  81035 TRUE       FALSE     
##  8 Asso~ A          Male             14           8 100102 TRUE       FALSE     
##  9 Asso~ A          Male              9           7  70000 TRUE       FALSE     
## 10 Asso~ A          Male             11           1 104800 TRUE       FALSE     
## 11 Asso~ A          Male             45          39  70700 TRUE       FALSE     
## 12 Asso~ A          Male             10           1 108413 TRUE       FALSE     
## 13 Asso~ A          Male             11           8 104121 TRUE       FALSE     
## 14 Asso~ B          Male             13          11 126431 TRUE       FALSE     
## 15 Prof  A          Male             29           7 204000 TRUE       FALSE     
## 16 Prof  A          Male             42          18 194800 TRUE       FALSE     
## 17 Prof  A          Male             43          43 205500 TRUE       FALSE     
## 18 Prof  B          Male             38          38 231545 TRUE       FALSE

Salaries$salary[Salaries$salary%in%obserwacje.odstajace$salary] <- NA
Salaries <- na.omit(Salaries)

##Zlogarytmowanie zmiennej salary

Salaries<-mutate(Salaries,nowa.kolumna=log10(salary))

####Sprawdzenie wartości odstających, normalności rozkładów i homegniczności wariancji

Salaries %>% 
  group_by(sex, rank, discipline) %>%
  identify_outliers(nowa.kolumna)

## # A tibble: 11 x 9
##    rank      discipline sex   yrs.since.phd yrs.service salary nowa.kolumna
##    <fct>     <fct>      <fct>         <int>       <int>  <int>        <dbl>
##  1 AsstProf  A          Male              8           3  69700         4.84
##  2 AsstProf  A          Male              5           3  69200         4.84
##  3 AsstProf  A          Male             11           4  78785         4.90
##  4 AssocProf A          Male             10           7  73877         4.87
##  5 AssocProf A          Male             30          23  74000         4.87
##  6 AssocProf A          Male             41          33  88600         4.95
##  7 AssocProf A          Male              8           6  88650         4.95
##  8 AssocProf A          Male             13           8  78182         4.89
##  9 AssocProf A          Male              8           5  86895         4.94
## 10 Prof      A          Male             51          51  57800         4.76
## 11 Prof      B          Male             46          45  67559         4.83
## # ... with 2 more variables: is.outlier <lgl>, is.extreme <lgl>

Salaries %>% 
  group_by(sex, rank, discipline) %>%
  shapiro_test(nowa.kolumna)

## # A tibble: 12 x 6
##    rank      discipline sex    variable     statistic      p
##    <fct>     <fct>      <fct>  <chr>            <dbl>  <dbl>
##  1 AsstProf  A          Female nowa.kolumna     0.813 0.104 
##  2 AsstProf  B          Female nowa.kolumna     0.896 0.387 
##  3 AssocProf A          Female nowa.kolumna     0.978 0.717 
##  4 AssocProf B          Female nowa.kolumna     0.916 0.517 
##  5 Prof      A          Female nowa.kolumna     0.936 0.575 
##  6 Prof      B          Female nowa.kolumna     0.976 0.939 
##  7 AsstProf  A          Male   nowa.kolumna     0.952 0.560 
##  8 AsstProf  B          Male   nowa.kolumna     0.931 0.0218
##  9 AssocProf A          Male   nowa.kolumna     0.891 0.0587
## 10 AssocProf B          Male   nowa.kolumna     0.981 0.840 
## 11 Prof      A          Male   nowa.kolumna     0.985 0.212 
## 12 Prof      B          Male   nowa.kolumna     0.986 0.236

ggqqplot(Salaries, "nowa.kolumna") +
  facet_grid(sex~rank~discipline)

wyniki2<-Salaries %>%
  levene_test(nowa.kolumna~sex*rank*discipline)
wyniki2

## # A tibble: 1 x 4
##     df1   df2 statistic        p
##   <int> <int>     <dbl>    <dbl>
## 1    11   366      8.92 4.20e-14

Na przyjętym przeze mnie poziomie istotności równym 0.01, nasze rozkłady grup można uznać za normalne(najmniejsze p jest równe 0.02). Niestety dalej posiadamy 2 obserwacje odstające oraz wariancje nie możemy uznać za jednorodne, jednakże uznałem, że przy liczbie obserwacji wynoszącej obecnie 378, test anova poradzi sobie dobrze, pomimo niespełnionych założeń. Nasze grupy nie są równoliczne, dlatego musimy użyć ANOVY typu 3.

###Anova typu 3

wyniki<-Salaries %>%
  anova_test(nowa.kolumna~sex*rank*discipline, type=c("III"))

## Coefficient covariances computed by hccm()

wyniki

## ANOVA Table (type tests)
## 
##                Effect DFn DFd      F                                p p<.05
## 1                 sex   1 366  0.270 0.603999999999999981348253186297      
## 2                rank   2 366 74.804 0.000000000000000000000000000579     *
## 3          discipline   1 366 27.566 0.000000258000000000000007365636     *
## 4            sex:rank   2 366  0.114 0.893000000000000015987211554602      
## 5      sex:discipline   1 366  0.635 0.425999999999999989785948173449      
## 6     rank:discipline   2 366  1.248 0.287999999999999978239628717347      
## 7 sex:rank:discipline   2 366  0.315 0.729999999999999982236431605997      
##        ges
## 1 0.000736
## 2 0.290000
## 3 0.070000
## 4 0.000621
## 5 0.002000
## 6 0.007000
## 7 0.002000

wyniki2<-Salaries %>%
 anova_test(nowa.kolumna~sex+rank+discipline, type=c("III"))

## Coefficient covariances computed by hccm()

wyniki2

## ANOVA Table (type tests)
## 
##       Effect DFn DFd       F
## 1        sex   1 373   0.362
## 2       rank   2 373 192.140
## 3 discipline   1 373  57.846
##                                                                p p<.05      ges
## 1 0.548000000000000042632564145606011152267456054687500000000000       0.000969
## 2 0.000000000000000000000000000000000000000000000000000000000439     * 0.507000
## 3 0.000000000000233000000000000016928299040319672030818765051663     * 0.134000

Odp: Wyniki obu analiz wariancji mówią nam, że na wysokość pensji(salary) istotny statystcznie wpływ mają dyscyplina i stopień naukowy(rank). Pomiędzy zmiennymi objaśniającymi nie wykryto żadnych istotnych statystycznie interakcji

###Czy lata od doktoratu (yrs.since.phd), staż pracy (yrs.service) mogą być istotnymi zmiennymi współ-oddziaływującymi na płacę (salary)?

Salaries$yrs.since.phd<-as.factor(Salaries$yrs.since.phd)
Salaries$yrs.service<-as.factor(Salaries$yrs.service)
wyniki4<- Salaries %>%
anova_test(nowa.kolumna~yrs.since.phd*yrs.service)

## Coefficient covariances computed by hccm()

## Note: model has aliased coefficients
##       sums of squares computed by model comparison

wyniki4

## ANOVA Table (type II tests)
## 
##                      Effect DFn DFd    F           p p<.05   ges
## 1             yrs.since.phd  49 120 3.09 0.000000291     * 0.558
## 2               yrs.service  48 120 2.71 0.000005950     * 0.520
## 3 yrs.since.phd:yrs.service 157 120 1.28 0.077000000       0.626

Salaries$yrs.since.phd<-as.numeric(Salaries$yrs.since.phd)
Salaries$yrs.service<-as.numeric(Salaries$yrs.service)
cor.test(Salaries$yrs.service,Salaries$yrs.since.phd)

## 
##  Pearson's product-moment correlation
## 
## data:  Salaries$yrs.service and Salaries$yrs.since.phd
## t = 43, df = 376, p-value <0.0000000000000002
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.893 0.927
## sample estimates:
##   cor 
## 0.912

Odp: Test pokazał, że i lata od doktoratu, i staż pracy istotnie wpływają na wysokość otrzymywanej płacy. Jak można było się domyślić, obie zmienne są współliniowe. Ich współczynnik korelacji Pearsona wyniósł aż 0.912

Czy istnieje istotna różnica w latach od doktoratu (yrs.since.phd) i stażu pracy (yrs.service) różnej rangi profesorów?

###Sprawdzamy normalność

Salaries %>%
  group_by(rank) %>%
  shapiro_test(yrs.service, yrs.since.phd) %>%
  arrange(variable)

## # A tibble: 6 x 4
##   rank      variable      statistic             p
##   <fct>     <chr>             <dbl>         <dbl>
## 1 AsstProf  yrs.service       0.922 0.000681     
## 2 AssocProf yrs.service       0.679 0.00000000136
## 3 Prof      yrs.service       0.980 0.000940     
## 4 AsstProf  yrs.since.phd     0.933 0.00204      
## 5 AssocProf yrs.since.phd     0.739 0.0000000192 
## 6 Prof      yrs.since.phd     0.970 0.0000266

Salaries %>%
  select(yrs.service, yrs.since.phd) %>%
  mshapiro_test()

## # A tibble: 1 x 2
##   statistic  p.value
##       <dbl>    <dbl>
## 1     0.877 8.10e-17

Niestety, rozkładom jest daleko od bycia normalnymi. Zarówno jednowymiarowo jak i wielowymiarowo

###Test Manova

wyniki6 <- manova(cbind(yrs.service, yrs.since.phd) ~ rank, data = Salaries)
summary(wyniki6)

##            Df Pillai approx F num Df den Df              Pr(>F)    
## rank        2  0.497       62      4    750 <0.0000000000000002 ***
## Residuals 375                                                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Odp: Tak, istnieje istotna statystycznie różnica w latach od doktoraru i stażu pracy profesorów różnej rangi, jednak należy pamiętać, że rozkłady naszych grup nie były normalne.Na dodatek, wcześniej wykazaliśmy, że zmienne yrs.service i yrs.since.phd są mocno skorelowane, co jest również przeszkodą dla użycia MANOVy.