Lai veiktu uzdevumus mājasdarbā tika veiktas dažādas datu manipulācijas, lai iegūtu noteiktas kategoriskas vērtības un datu kopas, pie kurām būtu iespējams veikt uzdevumus. Visas datu manipulācijas kodu veidā ir pieejamas zem dotā teksta, ja nu interesē:

interval_2010_2020 <- c("2010-11","2011-12","2012-13","2013-14",
                        "2014-15","2015-16","2016-17","2017-18","2018-19","2019-20")
# Vectors for each continent
north_america <- c("USA", "Canada", "Mexico", "Jamaica", "Puerto Rico", "US Virgin Islands", "Dominican Republic")
south_america <- c("Argentina", "Brazil", "Venezuela", "Uruguay", "Colombia")
europe <- c("Serbia and Montenegro", "Ukraine", "Croatia", "Lithuania", "Slovenia", "France", "Germany", 
            "England", "Turkey", "Greece", "Finland", "Ireland", "Scotland", "Poland", "Netherlands", 
            "Czech Republic", "Montenegro", "United Kingdom", "Yugoslavia", "Austria", "Switzerland", 
            "Spain", "Portugal", "Bosnia and Herzegovina", "Italy", "Denmark")
africa <- c("Nigeria", "Congo", "Democratic Republic of the Congo", "Cameroon", "Ghana", "Mali", 
            "Senegal", "Gabon", "Angola", "Egypt", "Sudan", "South Sudan", "Tunisia", "Cabo Verde")
asia_oceania <- c("China", "Japan", "South Korea", "Iran", "Georgia", "New Zealand", "Australia")


eastern_conference <- c("ATL", "BOS", "BKN", "CHA", "CHI", "CLE", "DET", "IND", "MIA", "MIL", "NYK", "ORL", "PHI", "TOR", "WAS")
western_conference <- c("DAL", "DEN", "GSW", "HOU", "LAC", "LAL", "MEM", "MIN", "NOP", "OKC", "PHX", "POR", "SAC", "SAS", "UTA")

grouped_data_conference <- data %>% 
  filter(team_abbreviation %in% c(eastern_conference, western_conference)) %>%
  mutate(conference = ifelse(team_abbreviation %in% eastern_conference, "East", "West")) %>%
  group_by(conference,season) %>%  
  summarise(
    count = n(),
    mean_height = mean(player_height, na.rm = TRUE),
    mean_weight = mean(player_weight, na.rm = TRUE),
    mean_pts = sum(pts, na.rm = TRUE)
  )
grouped_data_teams <- data %>% 
  filter(team_abbreviation %in% c(eastern_conference, western_conference)) %>%
  group_by(team_abbreviation,season) %>%  
  summarise(
    count = n(),
    mean_height = mean(player_height, na.rm = TRUE),
    mean_weight = mean(player_weight, na.rm = TRUE),
    mean_pts = sum(pts, na.rm = TRUE)
  ) %>% 
  mutate(conference = ifelse(team_abbreviation %in% eastern_conference, "East", "West"))

mean_val_conference <- grouped_data_conference %>% 
  group_by(conference) %>% 
  summarise(
    mean_mean_count = mean(count),
    mean_mean_weight = mean(mean_weight),
    mean_mean_height = mean(mean_height),
    mean_mean_pts = mean(mean_pts), 
    var_mean_height = var(mean_height),  # Variance of mean_height
    var_mean_weight = var(mean_weight),   # Variance of mean_weight (optional)
    var_mean_pts = var(mean_pts)
    
  )

east <- grouped_data_conference %>% 
  filter(conference %in% 'East')
west <- grouped_data_conference %>% 
  filter(conference %in% 'West')

mean_weight <- grouped_data_conference %>% 
  ungroup() %>% 
  select(mean_weight) %>% 
  mutate(mean_weight = as.numeric(as.character(mean_weight))) %>%
  pull(mean_weight)

east_mean_weight <- grouped_data_conference %>% 
  ungroup() %>% 
  filter(conference %in% 'East') %>% 
  select(mean_weight) %>% 
  mutate(mean_weight = as.numeric(as.character(mean_weight))) %>%
  pull(mean_weight) 

west_mean_weight <- grouped_data_conference %>% 
  ungroup() %>%
  filter(conference %in% 'West') %>% 
  select(mean_weight) %>% 
  mutate(mean_weight = as.numeric(as.character(mean_weight))) %>%
  pull(mean_weight) 


assign_continent <- function(country) {
  if (country %in% north_america) {
    return("North America")
  } else if (country %in% south_america) {
    return("South America")
  } else if (country %in% europe) {
    return("Europe")
  } else if (country %in% africa) {
    return("Africa")
  } else if (country %in% asia_oceania) {
    return("Asia/Oceania")
  } else {
    return(NA)  # Return NA if the country is not found
  }
}

# Adding a new column 'continent' to the data frame
data_continent_conference <- data %>%
  mutate(continent = sapply(country, assign_continent)) %>% 
  mutate(conference = ifelse(team_abbreviation %in% eastern_conference, "East", "West"))

player_years_summary <- data %>% 
  group_by(player_name) %>% 
  summarize(years_played = n())

data_continent <- data_continent_conference %>% 
  filter(net_rating > -15, net_rating < 15) %>% 
  # filter(continent == 'North America') %>%
  filter(conference == 'East')

1. uzdevums

(7 punkti). Izmantosim datus, kas atrasti pirmajā mājasdarbā. Izvēlēties atbilstošus datus pētāmajai problēmai (censties pārbaudīt tādas hipotēzes, kas būtu interesantas, nevis triviālas, kur, piemēram, grupu vienādību triviāli jānoraida) un veikt sekojošas hipotēžu pārbaudes:

Sadalot kādu mainīgo lielumu divās grupās, veikt divu izlašu salīdzināšanu H0 : µ1 = µ2 pret H1 : µ1 ̸= µ2. Izmantot gan t-testu, gan Vilkoksona neparametrisko rangu testu. Pārbaudīt datu normalitāti. Salīdzināt abu testu rezultātus. Aprēķināt efekta lielumu un atbilstošu jaudu, ja gribētu pētīt šādu efektu. Izdarīt secinājumus. Fiksējot jaudu 80% un mēģinot noteikt vidēju efektu, kādam būtu bijis jābūt datu apjomam šim datu piemēram?

Veicot datu manipulācijas, tika izveidotas divas skaitliskas listes ar nosaukumu east_mean_weight un west_mean_weight, kuru ieguvām grupējot visus datus pēc komandu konferencēm un veicot grupāciju pēc konferencēm un sezonām. Viena vērtība ir noteiktais vidējais spēlētāja svars konferencē uz noteikto sezonu. Ierakstu skaits abās listēs ir 27. Grafikā ir redzams density plot uz abām listēm.

#tagad testējam uz east_mean_weight un west_mean_weight

#density plots ja nu tas kko izsaka

ggplot()+
  geom_density(data=data.frame(mean_weight = east_mean_weight), aes(x=mean_weight), fill='steelblue',alpha=0.3)+
  geom_density(data=data.frame(mean_weight = west_mean_weight), aes(x=mean_weight), fill= 'tomato',alpha=0.3)


Izskatās diezgan līdzīgs, iespējams arī ir datu normalitāte. Tiks izmantots hipotēžu pārbaudē.

#datu normalitātes pārbaude, shapiro
# abi ir >pval , nevaram noraidīt, ka nav normalitāte

east_shapiro <- shapiro.test(east_mean_weight);east_shapiro
## 
##  Shapiro-Wilk normality test
## 
## data:  east_mean_weight
## W = 0.93535, p-value = 0.09357
west_shapiro <- shapiro.test(west_mean_weight);west_shapiro
## 
##  Shapiro-Wilk normality test
## 
## data:  west_mean_weight
## W = 0.9294, p-value = 0.0668


Abiem p-value > 0.05, nevaram noraidīt to, ka dati nav normāli. Tālāk tiek pildīti t-testi un vilkoksona testi.

# nenoraida datu normalitāti p-val > 0.05
variancee <- c(var(east_mean_weight),var(west_mean_weight))
# dispersija listēs ir 1.831185 1.977072

# T-testa darbošanās
# zinam ka dispersijas nav vienadas, variance is equal = false

t_tests <- t.test(east_mean_weight,west_mean_weight, var.equal = FALSE );t_tests
## 
##  Welch Two Sample t-test
## 
## data:  east_mean_weight and west_mean_weight
## t = 0.66164, df = 51.924, p-value = 0.5111
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.5051578  1.0021340
## sample estimates:
## mean of x mean of y 
##  100.5003  100.2519
# p>0.05 nenoraidam H0 hipotēzi, tad nepildam H1 hipotēzi

# Vilkoksona testa darbošanās

#te tas pats
willy_cox <- wilcox.test(east_mean_weight,west_mean_weight);willy_cox
## 
##  Wilcoxon rank sum exact test
## 
## data:  east_mean_weight and west_mean_weight
## W = 404, p-value = 0.5029
## alternative hypothesis: true location shift is not equal to 0

Abi testi mums ziņo, ka ir iespējams, ka abām listēm vidējā vērtība varētu būt vienāda, p-value > 0.05 . Vilkoksona testa pvērtība ir mazāka, tā tiecas vairāk noraidīt doto hipotēzi. Tālāk par efekta lielumu un jaudu, kādā efekts ir rēķināts.

# Efekta lieluma aprekinasana

effect <- (mean(east_mean_weight)-mean(west_mean_weight))/sd(east_mean_weight)
cat("efekts:", effect, "\n")
## efekts: 0.1836282
# Jaudas aprekinasana

power_t_test <- pwr.t.test(n = length(east_mean_weight), d = effect, sig.level = 0.05, power = NULL, type = "two.sample",alternative = "two.sided")
cat("jauda pie noteiktā efekta:", power_t_test$power, "\n")
## jauda pie noteiktā efekta: 0.1015698


Jauda ir ļoti zema, pie dotā datu skaita, būtu vai nu jāiegūst vairāk datu (sezonu skaits dotajā situācijā) vai arī jāizvēlas dati, kur efekta lielums būtu lielāks.
Tālāk par n skaitu, kur jauda un efekta lielums ir fiksēti pēc prasībām.

# datu apjoma noteikshaana

#atradām, kāds būtu vajadzīgais sample size uz to, kāds ir patreizējais efekta lielums lai būtu jēdzīgs tests.

power_t_test <- pwr.t.test(n = NULL, d = 0.5, sig.level = 0.05 , power = 0.8, type = "two.sample",alternative = "two.sided");power_t_test
## 
##      Two-sample t test power calculation 
## 
##               n = 63.76561
##               d = 0.5
##       sig.level = 0.05
##           power = 0.8
##     alternative = two.sided
## 
## NOTE: n is number in *each* group
power_wilcox <- round(power_t_test$n * 1.15,0)
cat('n skaits, ja būtu skatās pēc vilkoksona testa:', power_wilcox)
## n skaits, ja būtu skatās pēc vilkoksona testa: 73

Ja būtu mūsu datiem jauda 0.8 un efekta lielums 0.5, tad būtu vajadzīgas tikai 64 vai 73 vērtības lai būtu ticamas vērtības. 64 ir attiecīgi pēc t-testa un 73 ir attiecīgi pēc Vilkoksona testa.

Pārbaudīt datu normalitāti H0 : X ∼ N(µ, σ2) ar dažādiem testiem (izvēlēties vismaz trīs testus). Konstruēt p-vērtību līknes atkarībā no datu apjoma, izvēloties patvaļīgi datus no dotās kopas. Procedūru atkārtot, konstruējot vidējās p-vērtību līknes no kādiem 1000 atkārtojumiem. Izdarīt secinājumus. Datu normalitātes noraidīšanu/pieņemšanu pamatot ar atbilstošām aprakstošām statistikām;

Izvēlējos 4 dažādus testus, kuri tiks testēti uz mean_weight kolonnas. Saistībā ar to, ka vajag būt virs 7 vērtībām lai nemestu kļūdas testi, sāku ar 8. Rezultātā var redzēt, ka, jo lielāks n, jo lielāka iespēja, ka noraidīs H0, jeb ka dati ir normāli sadalīti. Pie n = 35 visi četri testi jau sāk noraidīt datu normalitāti.

shapiro_test <- function(data) {
  return(shapiro.test(data)$p.value)
}
lillie_test <- function(data) {
  return(lillie.test(data)$p.value)
}
cvm_test <- function(data) {
  return(cvm.test(data)$p.value)
}
ad_test <- function(data) {
  return(ad.test(data)$p.value)
}

# mean_weight #tiek izmantots dotais
n_values <- seq(8, length(mean_weight), by = 2)  # Sample sizes
repeats <- 100  # Number of repetitions

#tukši rezultāta vektori
shapiro_results <- numeric(length(n_values))
lillie_results <- numeric(length(n_values))
cvm_results <- numeric(length(n_values))
ad_results <- numeric(length(n_values))

# pie katra i (vertibu skaita, sākumā 8) veicam means aprēķinu uz sample p-vērtībām, iemetam results results vektorā
for (i in 1:length(n_values)) {
  n <- n_values[i]
  shapiro_p_values <- numeric(repeats)
  lillie_p_values <- numeric(repeats)
  cvm_p_values <- numeric(repeats)
  ad_p_values <- numeric(repeats)
  
  for (j in 1:repeats) {
    #sample
    sample_data <- sample(mean_weight, n, replace = TRUE)
    
    # testi
    shapiro_p_values[j] <- shapiro_test(sample_data)
    lillie_p_values[j] <- lillie_test(sample_data)
    cvm_p_values[j] <- cvm_test(sample_data)
    ad_p_values[j] <- ad_test(sample_data)
  }
  
  # Store average p-values for current sample size
  shapiro_results[i] <- mean(shapiro_p_values, na.rm = TRUE)
  lillie_results[i] <- mean(lillie_p_values, na.rm = TRUE)
  cvm_results[i] <- mean(cvm_p_values, na.rm = TRUE)
  ad_results[i] <- mean(ad_p_values, na.rm = TRUE)
}


results <- data.frame(
  n = n_values,
  Shapiro = shapiro_results,
  KS = lillie_results,
  CVM = cvm_results,
  AD = ad_results
)
# results_long

results_long <- melt(results, id.vars = "n", variable.name = "Test", value.name = "p_value")

# Plot average p-values for each test
ggplot(results_long, aes(x = n, y = p_value, color = Test)) +
  geom_line() +
  geom_hline(yintercept = 0.05, color = "red", linetype = "dashed") + 
  labs(title = "Vidējās p-vērtības pēc sample size, dažādiem normalitātes testiem uz mean_weights datiem",
       x = "n",
       y = "p-value",
       color = "Test")

Pārbaudīt kādu divu izvēlētu kategoriālu mainīgo lielumu neatkarību, izmantojot Hī-kvadrāta testu. Aprēķināt efekta lielumu un to interpretēt. Izdarīt atbilstošus secinājumus.

Izveidoju jaunu kolonnu, kontinents, kas ģeneralizē spēlētāja valsti pēc kontinenta ar nosaukumu ‘continent’. Izveidoju, citādi būtu jāsalīdzina liels kategorisku vērtību skaits. Tiks veikta divu ‘jaunizveidoto’ kategorisko vērtību, conference un continent salīdzināšana.

#pievienojam kontinentu uz datu kopas
data_continent_conference <- data %>%
  mutate(continent = sapply(country, assign_continent)) %>% 
  mutate(conference = ifelse(team_abbreviation %in% eastern_conference, "East", "West"))

tabel <- table(data_continent_conference$continent, data_continent_conference$conference);tabel
##                
##                 East West
##   Africa          98  102
##   Asia/Oceania    86  105
##   Europe         442  534
##   North America 5340 5703
##   South America   75   90
chisq <- chisq.test(tabel)

#lai aprēķinātu efektu
stat <- chisq.test(tabel)$statistic
nefekt <- sum(tabel) #laikam shitas ir skaits ar rows kam ir reali kontinents

efekts <- sqrt(stat/(2*nefekt))

print(efekts)
##  X-squared 
## 0.01354958

Rezultātā mēs iegūstam efekta lielumu: 0.0135, kas nozīmē ka abām grupām ir ļoti maza asociācija starp mainīgajiem. Tā ir pozitīva asociācija. Bet iespējams pozitīvā asociācija ir saistīta ar to, ka vienā konferencē ir tieši vairāk doto vērtību nekā otrā, ko var novērot no biežuma tabulas.

Salīdzināt vairākas grupas savā starpā, izmantojot ANOVA analīzi:

a. Veikt vienfaktora ANOVA analīzi. Pārbaudīt, vai parametriskie nosacījumi izpildās. Papildus pielietot Kruskala-Volisa neparametrisko testu. Veikt post-hoc analīzi, salīdzinot visas grupas pa pāriem. Noteikt un interpretēt efekta lielumu. Izdarīt secinājumus;

Izveidoju jaunu kolonnu ar nosaukumu draft_era pēc kuras grupējuma tiks izskatītas vidējās vērtības atšķirības kolonnā net_rating. Sākumā vienīgais ir svarīgi pārbaudīt, vai dati ir normāli.

#Priekš ANOVA pastamies, vai kontinenta dati ir devuši kkadas jaunas normalitātes

interval_1990_2000 <- 1990:1999
interval_2000_2010 <- 2000:2009
interval_2010_2020 <- 2010:2019


data_continent <- data_continent_conference %>% 
  filter(net_rating > -15, net_rating < 15) %>% 
  filter(continent == 'Europe') %>%
  filter(conference == 'East') %>% 
  mutate(draft_era = case_when(
    draft_year %in% interval_1990_2000 ~ "1990's",
    draft_year %in% interval_2000_2010 ~ "2000's",
    draft_year %in% interval_2010_2020 ~ "2010's",
    TRUE ~ "the rest"
  ))


ggplot(data_continent, aes(x=net_rating, color=draft_era))+
  geom_density()

qqnorm(data_continent$net_rating)
qqline(data_continent$net_rating, col ='red')

#diezgan ticams, ka ir datu normalitāte 

shapiro.test(data_continent$net_rating)
## 
##  Shapiro-Wilk normality test
## 
## data:  data_continent$net_rating
## W = 0.99504, p-value = 0.2054
#pieņemam, ka šitie dati ir normāli, par net rating, tagad skatamies citus prikolus

Dotajos grafikos un konsoles outputā varam redzēt pierādījumu, ka varam ticēt, ka dati ir normāli. Ik pēc katras grupas ir diezgan normāli, pēc qqplota varam secināt, ka ir diezgan normāli, kā arī shapiro tests nenoraida hipotēzi p-value > 0.05. Varam pildīt tālāk.

boxplot(net_rating~draft_era,data = data_continent, main= "kastu grafiki net_rating priekš dažādu ēru spēlētājiem \n no Eiropas austrumu konferencē", col=c('thistle','salmon','gold','gray'))


Te redzam kā izskatās kastu grafiki. visi ir diezgan līdzīgi, bet nosacīti var redzēt 2010 un the rest grupējumiem, ka ir mazliet zemāka mean vērtība.

#is.factor(data_continent$draft_era)

data_continent$draft_era <- factor(data_continent$draft_era)

mod<-aov(net_rating~draft_era,data = data_continent)
summary(mod)
##              Df Sum Sq Mean Sq F value Pr(>F)  
## draft_era     3    256   85.25   2.615 0.0507 .
## Residuals   411  13396   32.59                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Mēs nenoraidam to ka datu vidējās vērtības ir vienādas! Sakarā ar to, ka mums ir p > 0.05 tehniski, tad mums nav jaāveic post-hoc analīze jo nevaram noraidīt, ka means ir vienādi. Bet ja noapaļojam uz 2 vērtībām pēc komata, tad mums sanāk noraidīt hipotēzi, jo 0.5 !> 0.5. To arī skatamies tālāk.

# ja noapaļojam uz 2 vērtībām pēc komata, tad mums sanāk noraidīt hipotēzi, jo 0.5 !> 0.5.

#post hoc tests
#mums ir 4 grupas, jāveic 6 salīdzinājumi

#TukeyHSD(mod)

#nenoraida ja intervāls iekļauj 0, 
plot(TukeyHSD(mod))

#tapēc arī redzam, kapēc nenoraidas h0 hipotēze, bet ir diezgan tuvu, jo grupējumu salīdzinājuimi 2010's-2000's un the rest-2000's ir ļoti tuvu lai neiekļautu 0 ticamības intervālā.

kruskal.test(net_rating~draft_era,data = data_continent)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  net_rating by draft_era
## Kruskal-Wallis chi-squared = 8.9211, df = 3, p-value = 0.03036

Skatoties uz Tukey plot, redzam, kapēc nenoraidas h0 hipotēze, bet ir diezgan tuvu, jo grupējumu salīdzinājuimi 2010’s-2000’s un the rest-2000’s ir ļoti tuvu lai neiekļautu 0 ticamības intervālā. Ja neieskaitītu 0, tad mums arī tiktu noraidīta h0 hipotēze, ka grupu vidējie ir vienādi. Kruskala tests, kas ir neparametriskais tests, gan noraida datu vidējās vērtības vienādību, p-vērtība (0.0304) < 0.05.

b. Veikt ANOVA divfaktora analīzi atbilstoši izvēlētajiem mainīgajiem. Izmantot lekcijas slaidos ievietoto informāciju (attiecīgos r-kodus) un mēģināt interpretēt rezultātus. Lai izpildītu, papildināju data_continent ar jaunu vienkāršojumu, ka austrumu konferences komandas var grupēt pēc vēl ģeogrāfiskām divīzijām. Veidots ar mērķi, lai neizveidojas pārāk daudz grupējumu salīdzinājumu. Pētam tieši net_rating pēc grupējumiem division un draft_era.

atlantic <- c("BOS", "BKN", "NYK", "PHI", "TOR")
central <- c("CHI", "CLE", "DET", "IND", "MIL")
southeast <- c("ATL", "CHA", "MIA", "ORL", "WAS")


data_continent <- data_continent_conference %>% 
  filter(net_rating > -15, net_rating < 15) %>% 
  filter(continent == 'Europe') %>%
  filter(conference == 'East') %>% 
  mutate(draft_era = case_when(
    draft_year %in% interval_1990_2000 ~ "1990's",
    draft_year %in% interval_2000_2010 ~ "2000's",
    draft_year %in% interval_2010_2020 ~ "2010's",
    TRUE ~ "the rest"
  )) %>% 
  mutate(division = case_when(
    team_abbreviation %in% atlantic ~ "atlantic",
    team_abbreviation %in% central ~ "central",
    team_abbreviation %in% southeast ~ "southeast",
    TRUE ~ "Untitled"
  ))



ggplot(data = data_continent, aes(x = draft_era, y = net_rating, fill = division))+
  geom_boxplot()


Redzam te arī grupējumu box plotus.

#modelis bez mijiedarbības efekta
mod1 <- aov(net_rating ~ draft_era + division, data = data_continent)
summary(mod1)
##              Df Sum Sq Mean Sq F value Pr(>F)  
## draft_era     3    256   85.25   2.614 0.0508 .
## division      2     59   29.47   0.904 0.4059  
## Residuals   409  13337   32.61                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#modelis ar mijiedarbības efektu
mod2 <- aov(net_rating ~ draft_era * division, data = data_continent)
summary(mod2)
##                     Df Sum Sq Mean Sq F value  Pr(>F)   
## draft_era            3    256   85.25   2.698 0.04551 * 
## division             2     59   29.47   0.933 0.39432   
## draft_era:division   6    605  100.78   3.190 0.00453 **
## Residuals          403  12732   31.59                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# modeļu salīdzināšana
anova(mod1, mod2)
## Analysis of Variance Table
## 
## Model 1: net_rating ~ draft_era + division
## Model 2: net_rating ~ draft_era * division
##   Res.Df   RSS Df Sum of Sq      F   Pr(>F)   
## 1    409 13337                                
## 2    403 12732  6    604.69 3.1899 0.004532 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#ja p-vērtība<0.05, tad plašākais modelis labāks

Atsevišķi, abiem ir pvērtības virs 0.05, nevaram noteikt, vai abas vērtības iespaido net_rating iznākumu.
Ar mijiedarbību, ‘draft_era’ ir nozīmīgāks net_rating aprēķināšanai nekā ‘division’ (0.045 < 0.05, tikmēr 0.394 > 0.05).
Bet tāpatās sanāk, ka ir nozīmīgs mijiedarbības efekts, Pr(>F)(0.00453) < 0.01

#Aprēķinām aprakstošās statistikas pa grupām

sum <- summarySE(data_continent,
                 measurevar="net_rating",
                 groupvars=c("draft_era","division"))
sum
##    draft_era  division  N net_rating       sd        se       ci
## 1     1990's  atlantic 12 -3.3916667 5.614180 1.6206741 3.567080
## 2     1990's   central 27  0.6111111 5.853292 1.1264665 2.315485
## 3     1990's southeast  6 -1.4500000 7.212974 2.9446845 7.569552
## 4     2000's  atlantic 39 -1.9538462 4.405609 0.7054621 1.428133
## 5     2000's   central 63  0.4539683 6.005101 0.7565716 1.512365
## 6     2000's southeast 66 -0.3560606 5.480009 0.6745425 1.347154
## 7     2010's  atlantic 50 -0.1540000 5.372272 0.7597540 1.526783
## 8     2010's   central 54 -1.8222222 6.381922 0.8684696 1.741929
## 9     2010's southeast 52 -3.0750000 4.338378 0.6016247 1.207812
## 10  the rest  atlantic 18 -0.5888889 6.060663 1.4285119 3.013897
## 11  the rest   central 16 -5.0437500 6.047365 1.5118413 3.222413
## 12  the rest southeast 12 -2.7750000 7.408487 2.1386460 4.707128
#Zīmējam vidējās vērtības ar SE (standard error)

ggplot(sum, aes(x= draft_era, y=net_rating, color=division)) +
  geom_errorbar(aes(ymin=net_rating-se,
                    ymax=net_rating+se),
                width=.2, linewidth=0.7, position = position_dodge(0.3)) +
  geom_point(size=4, position = position_dodge(0.3))+
  geom_line(aes(x = draft_era, y = net_rating, group = division), position = position_dodge(0.3))

#mijiedarbības efektu raksturo taisnes
#mijiedarbības efekts ir nozīmīgs, ja taisnes nav paralēlas

#taisnes nav paralēlas, mijiedarbīgas efekts IR NOZĪMĪGS


Par nozīmīgumu mijiedarboties abām kategoriskām kolonnām ar net_rating arī pierāda dotā vizualizācija. Ja līnijas nekrustotos savā starpā, varētu minēt, ka nav nozīmīgs, bet tā kā ir krustojumi, tad mijiedarbības efekts ir nozīmīgs.

#salīdzinājumi pa pāriem
library(emmeans)
em<-emmeans(mod2,  ~ draft_era*division)
pairs(em)
##  contrast                               estimate   SE  df t.ratio p.value
##  1990's atlantic - 2000's atlantic        -1.438 1.86 403  -0.775  0.9998
##  1990's atlantic - 2010's atlantic        -3.238 1.81 403  -1.792  0.8218
##  1990's atlantic - the rest atlantic      -2.803 2.09 403  -1.338  0.9736
##  1990's atlantic - 1990's central         -4.003 1.95 403  -2.053  0.6568
##  1990's atlantic - 2000's central         -3.846 1.77 403  -2.172  0.5709
##  1990's atlantic - 2010's central         -1.569 1.79 403  -0.875  0.9993
##  1990's atlantic - the rest central        1.652 2.15 403   0.770  0.9998
##  1990's atlantic - 1990's southeast       -1.942 2.81 403  -0.691  0.9999
##  1990's atlantic - 2000's southeast       -3.036 1.76 403  -1.721  0.8578
##  1990's atlantic - 2010's southeast       -0.317 1.80 403  -0.176  1.0000
##  1990's atlantic - the rest southeast     -0.617 2.29 403  -0.269  1.0000
##  2000's atlantic - 2010's atlantic        -1.800 1.20 403  -1.499  0.9404
##  2000's atlantic - the rest atlantic      -1.365 1.60 403  -0.852  0.9995
##  2000's atlantic - 1990's central         -2.565 1.41 403  -1.823  0.8047
##  2000's atlantic - 2000's central         -2.408 1.15 403  -2.102  0.6214
##  2000's atlantic - 2010's central         -0.132 1.18 403  -0.111  1.0000
##  2000's atlantic - the rest central        3.090 1.67 403   1.852  0.7881
##  2000's atlantic - 1990's southeast       -0.504 2.46 403  -0.204  1.0000
##  2000's atlantic - 2000's southeast       -1.598 1.14 403  -1.407  0.9617
##  2000's atlantic - 2010's southeast        1.121 1.19 403   0.942  0.9986
##  2000's atlantic - the rest southeast      0.821 1.86 403   0.443  1.0000
##  2010's atlantic - the rest atlantic       0.435 1.55 403   0.281  1.0000
##  2010's atlantic - 1990's central         -0.765 1.34 403  -0.570  1.0000
##  2010's atlantic - 2000's central         -0.608 1.06 403  -0.571  1.0000
##  2010's atlantic - 2010's central          1.668 1.10 403   1.512  0.9367
##  2010's atlantic - the rest central        4.890 1.61 403   3.029  0.1042
##  2010's atlantic - 1990's southeast        1.296 2.43 403   0.534  1.0000
##  2010's atlantic - 2000's southeast        0.202 1.05 403   0.192  1.0000
##  2010's atlantic - 2010's southeast        2.921 1.11 403   2.624  0.2700
##  2010's atlantic - the rest southeast      2.621 1.81 403   1.451  0.9525
##  the rest atlantic - 1990's central       -1.200 1.71 403  -0.702  0.9999
##  the rest atlantic - 2000's central       -1.043 1.50 403  -0.694  0.9999
##  the rest atlantic - 2010's central        1.233 1.53 403   0.806  0.9997
##  the rest atlantic - the rest central      4.455 1.93 403   2.307  0.4735
##  the rest atlantic - 1990's southeast      0.861 2.65 403   0.325  1.0000
##  the rest atlantic - 2000's southeast     -0.233 1.49 403  -0.156  1.0000
##  the rest atlantic - 2010's southeast      2.486 1.54 403   1.617  0.9020
##  the rest atlantic - the rest southeast    2.186 2.09 403   1.044  0.9966
##  1990's central - 2000's central           0.157 1.29 403   0.122  1.0000
##  1990's central - 2010's central           2.433 1.32 403   1.837  0.7968
##  1990's central - the rest central         5.655 1.77 403   3.189  0.0669
##  1990's central - 1990's southeast         2.061 2.54 403   0.812  0.9997
##  1990's central - 2000's southeast         0.967 1.28 403   0.753  0.9998
##  1990's central - 2010's southeast         3.686 1.33 403   2.765  0.1994
##  1990's central - the rest southeast       3.386 1.95 403   1.736  0.8503
##  2000's central - 2010's central           2.276 1.04 403   2.184  0.5625
##  2000's central - the rest central         5.498 1.57 403   3.494  0.0261
##  2000's central - 1990's southeast         1.904 2.40 403   0.793  0.9997
##  2000's central - 2000's southeast         0.810 0.99 403   0.818  0.9996
##  2000's central - 2010's southeast         3.529 1.05 403   3.351  0.0412
##  2000's central - the rest southeast       3.229 1.77 403   1.824  0.8041
##  2010's central - the rest central         3.222 1.60 403   2.014  0.6840
##  2010's central - 1990's southeast        -0.372 2.42 403  -0.154  1.0000
##  2010's central - 2000's southeast        -1.466 1.03 403  -1.422  0.9588
##  2010's central - 2010's southeast         1.253 1.09 403   1.147  0.9923
##  2010's central - the rest southeast       0.953 1.79 403   0.531  1.0000
##  the rest central - 1990's southeast      -3.594 2.69 403  -1.336  0.9740
##  the rest central - 2000's southeast      -4.688 1.57 403  -2.993  0.1145
##  the rest central - 2010's southeast      -1.969 1.61 403  -1.225  0.9867
##  the rest central - the rest southeast    -2.269 2.15 403  -1.057  0.9962
##  1990's southeast - 2000's southeast      -1.094 2.40 403  -0.456  1.0000
##  1990's southeast - 2010's southeast       1.625 2.42 403   0.671  0.9999
##  1990's southeast - the rest southeast     1.325 2.81 403   0.471  1.0000
##  2000's southeast - 2010's southeast       2.719 1.04 403   2.609  0.2783
##  2000's southeast - the rest southeast     2.419 1.76 403   1.371  0.9683
##  2010's southeast - the rest southeast    -0.300 1.80 403  -0.167  1.0000
## 
## P value adjustment: tukey method for comparing a family of 12 estimates


Te varam redzēt grupu pāru salīdzinājumus. Mērķis tam ir atrast vietas, kur p-val <0.05 . Tas ziņotu, ka tur ir manāmas atšķirības vidējām vērtībām starp grupām. Tas ir manāms diviem pāriem: ‘2000’s central - 2010’s southeast’ un ‘2000’s central - the rest central’.

2.Uzdevums

(2 punkti). Pieņemsim, ka doti novērojumu pāri (x1, y1), . . . (xn, yn), kas iegūti no diviem nepārtrauktiem skaitliskiem gadījuma lielumiem. Izveidot programmā R funkciju, kas sadala datus kategorijās (intervālos) un izveido attiecīgu biežuma tabulu. Pielietojot hī-kvadrāta testu, šai funkcijai būtu jāizdod attiecīgo p-vērtību dažādu kategoriju izvēlei (var izveidot 2x2, 3x3 ect. tabulas un pētīt, kā mainās attiecīgā p-vērtība). Izvēlēties kādu datu piemēru no 1.mājasdarba datiem un pielietot šo testu neatkarības pārbaudei. Izdarīt secinājumus. Vai korelācijas analīze dod līdzīgu rezultātu?

Izveidoju funkciju, kas izveido biežuma tabulu un hī-kvadrātu testa rezultātus. Tika izveidota datu filtrācija no md 1 datiem, kur tiek izpētīta korelācija starp vidējo punktu un piespēļu skaitam starp top10 drafta izvēles spēlētājiem no rietumu konferences.

freq_table_chisquare <- function(x, y, bins_x, bins_y) {
  
  x_cut <- cut(x, breaks = bins_x, labels = FALSE)
  y_cut <- cut(y, breaks = bins_y, labels = FALSE)
  
  freq_table <- table(x_cut, y_cut)
  
  chi_test <- chisq.test(freq_table)
  
  return(list(freq_table = freq_table, chi_square_test = chi_test))
}

#ņemu datus no grupējuma, kur x ir pirmā raunda izvēles spēlētāju vidējais punktu skaits pa sezonām pēc austrumu konferences un y ir tas pats bet viņu vidējo piespēļu skaits

first_round_data_west <- data %>% 
  mutate(conference = ifelse(team_abbreviation %in% eastern_conference, "East", "West")) %>%
  filter(draft_number < 10) %>% 
  filter(conference %in% 'West')

x_data <- first_round_data_west$pts
y_data <- first_round_data_west$ast 

data_res_2x2 <- freq_table_chisquare(x_data,y_data,2,2)
data_res_3x3 <- freq_table_chisquare(x_data,y_data,3,3)

# Izdrukājam rezultātus uz diviem normāliem sadalījumiem
{
cat("2x2 tabula uz md1 datiem: \n")
print(data_res_2x2$freq_table)
cat("hī kvadrātu testa rezultāts")
print(data_res_2x2$chi_square_test)
}
## 2x2 tabula uz md1 datiem: 
##      y_cut
## x_cut   1   2
##     1  68   0
##     2 100  11
## hī kvadrātu testa rezultāts
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  freq_table
## X-squared = 5.5646, df = 1, p-value = 0.01833
{
cat("3x3 tabula uz md1 datiem: \n")
print(data_res_3x3$freq_table)
cat("hī kvadrātu testa rezultāts \n")
print(data_res_3x3$chi_square_test)
}
## 3x3 tabula uz md1 datiem: 
##      y_cut
## x_cut  1  2  3
##     1 41  0  0
##     2 67  5  0
##     3 34 26  6
## hī kvadrātu testa rezultāts 
## 
##  Pearson's Chi-squared test
## 
## data:  freq_table
## X-squared = 50.774, df = 4, p-value = 2.489e-10
cor_test <- cor.test(x_data, y_data);cor_test
## 
##  Pearson's product-moment correlation
## 
## data:  x_data and y_data
## t = 11.92, df = 177, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.5771153 0.7414095
## sample estimates:
##       cor 
## 0.6673025
# cik es saprotu, mēs sadalam sadalījumu intervālos (šajā gadījumā 2) un tad veicot frequency table, parādas tikai parādas 100 vērtības, kaut ir kopā 200.
# Tas notiek tapēc, jo skatās pēc x_cut intervāla un y_cut intervāla, vai iekļaujas 1:1 vai 1:2. Domu sapratu

#uz 3x3 izmet warning jo uz dažām freq vērtībām ir maza iespēja un funkcijai nepatīk

Salīdzinot gan 2x2 gan 3x3, redzam, ka abiem p-vērtība ir <0.05, kas nozīmē ka ir statistiski nozīmīga asociācija starp pts un ast vērtībām. Protams, veicot 3x3, mums p-vērtība drastiski paliek mazāka, no 0.01833 uz 2.489e-10.
Veicot korelācijas testu, redzam, ka mums vērtība iznāk 0.6673025, kas ir nozīmīga, pozitīva un liela korelācija. (diapazons ir -1:1)

3. uzdevums

(1 punkts). Pieņemsim, ka mēs vēlamies saprast, kurš tests datu normalitātei labāks: Kolmogorva- Smirnova, Šapiro iebūvētais vai Pīrsona (Hī-kvadrāta) tests. Būtu nepieciešams simulēt tā saukto empīrisko jaudu pret kādu alternatīvu un apskatīties, kuram testam jauda ātrāk aug uz 1. Alternatīvu nevarētu ņemt pārāk tālu no nulles hipotēzes (ne-pārāk atšķirīgu no normālā sadalīuma, citādi testi triviāli noraidīs), piemēram, var izvēlēties U[−1, 1] sadalījumu. Atlikt uz x-ass datu apjomu ģenerētajiem datiem un uz y-ass simulēto jaudu visiem testiem.

Tika izveidota funkcija, kur pārbauda uz noteikta n sample skaita kāda ir varbūtība, ka tiks noraidīta normalitāte uz 1000 simulācijām. Tika izmantots vienmērīgs sadalījums U[0,1]. Pēc tam tika atkārtots pie cita n sample skaita. Tādā veidā no 10 līdz 350 n_sample lielumam. Lai redzētu momentu, kur empīriskā jauda sasniedz 100. Papildus tika izveidota filtrācija kura atrod pirmo instanci, kur jauda == 100.

simulate_power_visiem <- function(n){
  n_simulations = 1000
  h0 <- 0.05
  
  rejections_Shapiro <- 0
  rejections_lillie <- 0
  rejections_Pearson <- 0
  
  for(i in 1:n_simulations){
    U_data <- runif(n,0,1) 
    
    p_shap <- shapiro.test(U_data)$p.value
    p_lillie <- lillie.test(U_data)$p.value
    p_Pearson <- pearson.test(U_data)$p.value
    
    if(p_shap < h0) rejections_Shapiro <- rejections_Shapiro +1
    if(p_lillie < h0) rejections_lillie <- rejections_lillie +1
    if(p_Pearson < h0) rejections_Pearson <- rejections_Pearson +1
  }
  
  jauda_shap <- (rejections_Shapiro/n_simulations)*100
  jauda_lillie <- (rejections_lillie/n_simulations)*100
  jauda_Pearson <- (rejections_Pearson/n_simulations)*100
  
  return( c(jauda_shap,jauda_lillie,jauda_Pearson))
}


n_skaits <- seq(10, 350, 10)
results_df <- data.frame(n=numeric(), Shapiro = numeric(), KS = numeric(), Pearson = numeric())

for (n in n_skaits){
  rezultats <- simulate_power_visiem(n=n)
  
  results_df <- rbind(results_df, data.frame(n=n, Shapiro= rezultats[1], KS = rezultats[2], Pearson = rezultats[3]))
}

#lai var izmantot us plota
results_long <- results_df %>% 
  pivot_longer(cols = c(Shapiro,KS,Pearson),
               names_to = 'test',
               values_to = 'p_value')

#lai atrastu uz x ass vietas, kur vairs nemainās jaudas vērtība uz katru testu 
result <- results_long %>%
  filter(p_value == 100) %>%
  group_by(test) %>%
  summarise(min_n = min(n)) %>%
  ungroup()


ggplot(results_long, aes(x=n, y=p_value, color = test))+
  geom_line()+
  labs(title = 'Testu jauda pie noteiktas datu kopas izmēra',
       x = 'Ierakstu skaits',
       y = 'Varbūtība, ka tiks noraidīta normalitāte, alpha = 0.05')+
 geom_vline(data = result, aes(xintercept = min_n, color = test), linetype = 'dashed')  # Add vertical lines where p_value = 100


Skatoties uz rezultātiem, varam redzēt, ka Shapiro testam jauda sasniegs maksimumu ātrāk nekā pārējiem testiem (pie n = 120), tikmēr pīrsona un KS testi nonāca pie gala jaudas samērā vienlīdzīgi, (310 un 330). Secinam, ka Shapiro tests ir jūtīgāks uz datu normalitāti nekā KS un Pīrsona testi.