# Data Generation

library(tidyverse)
## Warning: package 'tidyverse' was built under R version 4.1.1
## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --
## v ggplot2 3.4.2     v purrr   1.0.1
## v tibble  3.2.1     v dplyr   1.1.2
## v tidyr   1.3.0     v stringr 1.5.0
## v readr   2.1.2     v forcats 0.5.1
## Warning: package 'ggplot2' was built under R version 4.1.3
## Warning: package 'tibble' was built under R version 4.1.3
## Warning: package 'tidyr' was built under R version 4.1.3
## Warning: package 'readr' was built under R version 4.1.3
## Warning: package 'purrr' was built under R version 4.1.3
## Warning: package 'dplyr' was built under R version 4.1.3
## Warning: package 'stringr' was built under R version 4.1.3
## Warning: package 'forcats' was built under R version 4.1.1
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(broom)


city <- sample(x = c("Tel Aviv", "Jslm", "NY", "Boston", "LA"), size = 1500, replace = T)
age <- sample(x = seq(15, 105), size = 1500, replace = T)
Gender <- sample(x = c(1,2), size = 1500, replace = T)
Gender_word <- ifelse(Gender == 1, "Male", "Female")
Salary <- sample(x = seq(2500,10000), size = 1500, replace = T)
Seniority <- sample(x = seq(1,18), size = 1500, replace = T)

data <- data.frame(age, city, Gender, Gender_word, Salary, Seniority)


library(apaTables)
## Warning: package 'apaTables' was built under R version 4.1.3
# One Way ANOVA
apa.1way.table(dv = Salary, iv = Gender_word, filename = "Result_t", data = data)
## 
## 
## Descriptive statistics for Salary as a function of Gender_word.  
## 
##  Gender_word       M      SD
##       Female 6347.06 2205.34
##         Male 6303.27 2154.79
## 
## Note. M and SD represent mean and standard deviation, respectively.
## 
apa.1way.table(dv = Salary, iv = Gender_word, filename = "Result_t", data = data)
## 
## 
## Descriptive statistics for Salary as a function of Gender_word.  
## 
##  Gender_word       M      SD
##       Female 6347.06 2205.34
##         Male 6303.27 2154.79
## 
## Note. M and SD represent mean and standard deviation, respectively.
## 
m1 <- lm(Salary ~ Gender_word + Seniority, data)

# Regression
apaTables::apa.reg.table(m1, filename = "Regression_result")
## 
## 
## Regression results using Salary as the criterion
##  
## 
##        Predictor         b           b_95%_CI sr2  sr2_95%_CI             Fit
##      (Intercept) 6527.31** [6269.91, 6784.70]                                
##  Gender_wordMale    -41.95  [-262.64, 178.74] .00 [-.00, .00]                
##        Seniority    -19.04     [-40.59, 2.50] .00 [-.00, .01]                
##                                                                     R2 = .002
##                                                               95% CI[.00,.01]
##                                                                              
## 
## Note. A significant b-weight indicates the semi-partial correlation is also significant.
## b represents unstandardized regression weights. 
## sr2 represents the semi-partial correlation squared.
## Square brackets are used to enclose the lower and upper limits of a confidence interval.
## * indicates p < .05. ** indicates p < .01.
## 
av1 <- aov(Salary ~ city ,data)

av1 %>% summary()
##               Df    Sum Sq Mean Sq F value Pr(>F)
## city           4 1.521e+07 3803645     0.8  0.525
## Residuals   1495 7.104e+09 4751706
# 
apa.2way.table(city, Gender_word, Salary, show.marginal.means = T, data, 
               filename = "two-way result")
## 
## 
## Means and standard deviations for Salary as a function of a 5(city) X 2(Gender_word) design 
## 
##           Gender_word                                         
##                Female            Male         Marginal        
##      city           M      SD       M      SD        M      SD
##    Boston     5927.99 2292.16 6423.01 2231.19  6174.66 2271.70
##      Jslm     6357.33 2128.01 6332.41 2175.63  6342.99 2151.89
##        LA     6204.81 2255.00 6291.34 2223.68  6245.32 2237.05
##        NY     6839.16 2105.69 6032.18 2054.15  6411.30 2113.86
##  Tel Aviv     6432.70 2155.71 6452.95 2093.51  6442.12 2123.66
##  Marginal     6347.06 2205.34 6303.27 2154.79                 
## 
## Note. M and SD represent mean and standard deviation, respectively. 
## Marginal indicates the means and standard deviations pertaining to main effects.
apa.aov.table(lm_output = av1, filename = "aov_res", 
              table.number = 12, conf.level = 0.05)
## 
## 
## Table 12 
## 
## ANOVA results using Salary as the dependent variable
##  
## 
##    Predictor             SS   df             MS       F    p partial_eta2
##  (Intercept) 11247287059.58    1 11247287059.58 2367.00 .000             
##         city    15214579.48    4     3803644.87    0.80 .525          .00
##        Error  7103800953.43 1495     4751706.32                          
##  CI_90_partial_eta2
##                    
##          [.00, .00]
##                    
## 
## Note: Values in square brackets indicate the bounds of the 90% confidence interval for partial eta-squared
test_res <- t.test(Salary ~ Gender_word, data)

# add_r_file.docx <- read_docx()

broom::tidy(m1) %>% knitr::kable(digits = 3,align = 'c', 
                                 row.names = T, format.args = list(big.mark = ','), 
                                 caption = "Caption for this thing")
Caption for this thing
term estimate std.error statistic p.value
1 (Intercept) 6,527.305 131.218 49.744 0.000
2 Gender_wordMale -41.947 112.508 -0.373 0.709
3 Seniority -19.044 10.983 -1.734 0.083 
# table(data$city, data$Gender_word) %>% knitr::kable(align = "c") %>% 
#   body_add_table(x = "add_r_file.docx")

library(officer)
## Warning: package 'officer' was built under R version 4.1.3
t_test_result <- t.test(Salary ~ Gender_word, data = data)

# Display results in a simple format
cat("A two-sample t-test was conducted to compare payment between genders.\n")
## A two-sample t-test was conducted to compare payment between genders.
cat(sprintf("t(%.0f) = %.2f, p = %.3f, 95%% CI [%.2f, %.2f]\n",
            t_test_result$parameter,
            t_test_result$statistic,
            t_test_result$p.value,
            t_test_result$conf.int[1],
            t_test_result$conf.int[2]))
## t(1494) = 0.39, p = 0.697, 95% CI [-177.10, 264.69]
t_test_result <- t.test(Salary ~ Gender_word, data = data)

# Create a Word document
doc <- read_docx() %>%
  body_add_par("T-Test Results", style = "heading 1") %>%
  body_add_par(sprintf("A two-sample t-test was conducted to compare payment between genders.\n")) %>%
  body_add_par(sprintf("t(%.0f) = %.2f, p = %.3f, 95%% CI [%.2f, %.2f]", 
                       t_test_result$parameter, 
                       t_test_result$statistic, 
                       t_test_result$p.value, 
                       t_test_result$conf.int[1], 
                       t_test_result$conf.int[2]))

# Save the document
print(doc, target = "t_test_results.docx")


tidy(aov(Salary~Gender_word + city + city*Gender_word, data)) %>% 
  knitr::kable(caption = "2-way-anova table", row.names = T)
2-way-anova table
term df sumsq meansq statistic p.value
1 Gender_word 1 719079.2 719079.2 0.1522616 0.6964392
2 city 4 15284168.5 3821042.1 0.8090876 0.5192950
3 Gender_word:city 4 66255021.0 16563755.2 3.5072967 0.0073820
4 Residuals 1490 7036757264.2 4722655.9 NA NA 
ggplot(data, aes(y = Salary, fill = city)) + geom_boxplot(size = 0.55) + 
  facet_wrap(~Gender_word)

library(knitr)
## Warning: package 'knitr' was built under R version 4.1.2
m1 %>% tidy() %>% kable(digits = 3, row.names = T, caption = "This is the label")
This is the label
term estimate std.error statistic p.value
1 (Intercept) 6527.305 131.218 49.744 0.000
2 Gender_wordMale -41.947 112.508 -0.373 0.709
3 Seniority -19.044 10.983 -1.734 0.083 
glance(data) %>% kable(caption = "This is too", row.names = T)
## Warning: Data frame tidiers are deprecated and will be removed in an upcoming
## release of broom.
This is too
nrow ncol complete.obs na.fraction
1 1500 6 1500 0 
data.frame("Residuals"=m1$residuals, " BREAK ", 
           "Fitted Values"=m1$fitted.values, "Sequence"=seq(1,6000), "sss"=seq(1,50)+seq(2,17)) %>% head(50) %>% 
  kable(digits = 3, caption = "Values", align = "c", row.names = T)
## Warning in seq(1, 50) + seq(2, 17): longer object length is not a multiple of
## shorter object length
## Warning in data.frame(Residuals = m1$residuals, " BREAK ", `Fitted Values` =
## m1$fitted.values, : row names were found from a short variable and have been
## discarded
Values
Residuals X..BREAK.. Fitted.Values Sequence sss
1 2523.297 BREAK 6199.703 1 3
2 -892.271 BREAK 6447.271 2 5
3 -90.834 BREAK 6256.834 3 7
4 -1844.174 BREAK 6470.174 4 9
5 -1263.218 BREAK 6489.218 5 11
6 2782.729 BREAK 6447.271 6 13
7 -43.174 BREAK 6470.174 7 15
8 -1690.315 BREAK 6466.315 8 17
9 2666.991 BREAK 6333.009 9 19
10 769.001 BREAK 6393.999 10 21
11 -1801.572 BREAK 6142.572 11 23
12 1030.253 BREAK 6218.747 12 25
13 3643.393 BREAK 6222.607 13 27
14 2603.782 BREAK 6489.218 14 29
15 -3768.922 BREAK 6294.922 15 31
16 -697.660 BREAK 6180.660 16 33
17 -2121.912 BREAK 6355.912 17 19
18 -1716.227 BREAK 6428.227 18 21
19 1001.175 BREAK 6317.825 19 23
20 -1577.999 BREAK 6393.999 20 25
21 -2625.834 BREAK 6256.834 21 27
22 2652.816 BREAK 6409.184 22 29
23 1452.166 BREAK 6256.834 23 31
24 -300.965 BREAK 6313.965 24 33
25 2437.035 BREAK 6313.965 25 35
26 2229.957 BREAK 6413.043 26 37
27 2508.393 BREAK 6222.607 27 39
28 -3224.563 BREAK 6203.563 28 41
29 -1819.130 BREAK 6451.130 29 43
30 -1880.650 BREAK 6241.650 30 45
31 -2110.130 BREAK 6451.130 31 47
32 -311.184 BREAK 6409.184 32 49
33 1891.001 BREAK 6393.999 33 35
34 1278.991 BREAK 6333.009 34 37
35 1735.739 BREAK 6508.261 35 39
36 88.685 BREAK 6466.315 36 41
37 -2025.825 BREAK 6317.825 37 43
38 586.088 BREAK 6355.912 38 45
39 -955.130 BREAK 6451.130 39 47
40 804.860 BREAK 6390.140 40 49
41 2128.816 BREAK 6409.184 41 51
42 -3778.130 BREAK 6451.130 42 53
43 -1995.174 BREAK 6470.174 43 55
44 1332.262 BREAK 6279.738 44 57
45 -3862.087 BREAK 6432.087 45 59
46 -468.703 BREAK 6199.703 46 61
47 841.437 BREAK 6203.563 47 63
48 7.350 BREAK 6241.650 48 65
49 301.913 BREAK 6432.087 49 51
50 2269.393 BREAK 6222.607 50 53