AST 518
Md Monowar Hossain
2025-08-17
Basic R
17/08/2025
# rm(list = is())
mtcars## mpg cyl disp hp drat wt qsec vs am gear carb
## Mazda RX4 21.0 6 160.0 110 3.90 2.620 16.46 0 1 4 4
## Mazda RX4 Wag 21.0 6 160.0 110 3.90 2.875 17.02 0 1 4 4
## Datsun 710 22.8 4 108.0 93 3.85 2.320 18.61 1 1 4 1
## Hornet 4 Drive 21.4 6 258.0 110 3.08 3.215 19.44 1 0 3 1
## Hornet Sportabout 18.7 8 360.0 175 3.15 3.440 17.02 0 0 3 2
## Valiant 18.1 6 225.0 105 2.76 3.460 20.22 1 0 3 1
## Duster 360 14.3 8 360.0 245 3.21 3.570 15.84 0 0 3 4
## Merc 240D 24.4 4 146.7 62 3.69 3.190 20.00 1 0 4 2
## Merc 230 22.8 4 140.8 95 3.92 3.150 22.90 1 0 4 2
## Merc 280 19.2 6 167.6 123 3.92 3.440 18.30 1 0 4 4
## Merc 280C 17.8 6 167.6 123 3.92 3.440 18.90 1 0 4 4
## Merc 450SE 16.4 8 275.8 180 3.07 4.070 17.40 0 0 3 3
## Merc 450SL 17.3 8 275.8 180 3.07 3.730 17.60 0 0 3 3
## Merc 450SLC 15.2 8 275.8 180 3.07 3.780 18.00 0 0 3 3
## Cadillac Fleetwood 10.4 8 472.0 205 2.93 5.250 17.98 0 0 3 4
## Lincoln Continental 10.4 8 460.0 215 3.00 5.424 17.82 0 0 3 4
## Chrysler Imperial 14.7 8 440.0 230 3.23 5.345 17.42 0 0 3 4
## Fiat 128 32.4 4 78.7 66 4.08 2.200 19.47 1 1 4 1
## Honda Civic 30.4 4 75.7 52 4.93 1.615 18.52 1 1 4 2
## Toyota Corolla 33.9 4 71.1 65 4.22 1.835 19.90 1 1 4 1
## Toyota Corona 21.5 4 120.1 97 3.70 2.465 20.01 1 0 3 1
## Dodge Challenger 15.5 8 318.0 150 2.76 3.520 16.87 0 0 3 2
## AMC Javelin 15.2 8 304.0 150 3.15 3.435 17.30 0 0 3 2
## Camaro Z28 13.3 8 350.0 245 3.73 3.840 15.41 0 0 3 4
## Pontiac Firebird 19.2 8 400.0 175 3.08 3.845 17.05 0 0 3 2
## Fiat X1-9 27.3 4 79.0 66 4.08 1.935 18.90 1 1 4 1
## Porsche 914-2 26.0 4 120.3 91 4.43 2.140 16.70 0 1 5 2
## Lotus Europa 30.4 4 95.1 113 3.77 1.513 16.90 1 1 5 2
## Ford Pantera L 15.8 8 351.0 264 4.22 3.170 14.50 0 1 5 4
## Ferrari Dino 19.7 6 145.0 175 3.62 2.770 15.50 0 1 5 6
## Maserati Bora 15.0 8 301.0 335 3.54 3.570 14.60 0 1 5 8
## Volvo 142E 21.4 4 121.0 109 4.11 2.780 18.60 1 1 4 2ncol(mtcars)## [1] 11length(mtcars)## [1] 11typeof(mtcars)## [1] "list"typeof(x=5)## [1] "double"f3<-function(x,y){
r<- 3*x + x*y + log(x = x,base = exp(1))
return(r)
}
f3(x = 2, y = 2)## [1] 10.69315x<-2
typeof(x)## [1] "double"y<-1L
typeof(y)## [1] "integer"y2<-"5"
typeof(y2)## [1] "character"xv<-c(11,9,8,10,5)
xc<-c("boy","boy","girl","boy","girl")
#table(gender)
# prop.table(gender)
#genf<- as.factor(gender)
#as.numeric(gender)
#as.numeric(genf)
ind <- c(T, T, F, F, T)
xd<- data.frame(age=xv,gender=xc)
xd$gender=="boy"## [1] TRUE TRUE FALSE TRUE FALSExd[xd$gender=="boy",]## age gender
## 1 11 boy
## 2 9 boy
## 4 10 boyrm(xv)
rm(xc)
xd[xd$gender=="boy",]## age gender
## 1 11 boy
## 2 9 boy
## 4 10 boymissing value
xm<-c(1:5,NA)
is.na(xm)## [1] FALSE FALSE FALSE FALSE FALSE TRUEsum (is.na(xm))## [1] 1sum (!is.na(xm))## [1] 5na.omit(xm)## [1] 1 2 3 4 5
## attr(,"na.action")
## [1] 6
## attr(,"class")
## [1] "omit"xv<-c(11,9,8,10,5)
xc<-c("boy","boy","girl","boy",NA)
xd<- data.frame(x=xv,y=xc)
xd[!is.na(xd$y),]## x y
## 1 11 boy
## 2 9 boy
## 3 8 girl
## 4 10 boyAdvanced R
R for data Science
Problem 1
library(tidyverse)## Warning: package 'tidyverse' was built under R version 4.4.3## Warning: package 'lubridate' was built under R version 4.4.2## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.1 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(gapminder)
d <- gapminder
dim(d)## [1] 1704 6str(d)## tibble [1,704 × 6] (S3: tbl_df/tbl/data.frame)
## $ country : Factor w/ 142 levels "Afghanistan",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ continent: Factor w/ 5 levels "Africa","Americas",..: 3 3 3 3 3 3 3 3 3 3 ...
## $ year : int [1:1704] 1952 1957 1962 1967 1972 1977 1982 1987 1992 1997 ...
## $ lifeExp : num [1:1704] 28.8 30.3 32 34 36.1 ...
## $ pop : int [1:1704] 8425333 9240934 10267083 11537966 13079460 14880372 12881816 13867957 16317921 22227415 ...
## $ gdpPercap: num [1:1704] 779 821 853 836 740 ...library(palmerpenguins)
penguins## # A tibble: 344 × 8
## species island bill_length_mm bill_depth_mm flipper_length_mm body_mass_g
## <fct> <fct> <dbl> <dbl> <int> <int>
## 1 Adelie Torgersen 39.1 18.7 181 3750
## 2 Adelie Torgersen 39.5 17.4 186 3800
## 3 Adelie Torgersen 40.3 18 195 3250
## 4 Adelie Torgersen NA NA NA NA
## 5 Adelie Torgersen 36.7 19.3 193 3450
## 6 Adelie Torgersen 39.3 20.6 190 3650
## 7 Adelie Torgersen 38.9 17.8 181 3625
## 8 Adelie Torgersen 39.2 19.6 195 4675
## 9 Adelie Torgersen 34.1 18.1 193 3475
## 10 Adelie Torgersen 42 20.2 190 4250
## # ℹ 334 more rows
## # ℹ 2 more variables: sex <fct>, year <int>str(penguins)## tibble [344 × 8] (S3: tbl_df/tbl/data.frame)
## $ species : Factor w/ 3 levels "Adelie","Chinstrap",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ island : Factor w/ 3 levels "Biscoe","Dream",..: 3 3 3 3 3 3 3 3 3 3 ...
## $ bill_length_mm : num [1:344] 39.1 39.5 40.3 NA 36.7 39.3 38.9 39.2 34.1 42 ...
## $ bill_depth_mm : num [1:344] 18.7 17.4 18 NA 19.3 20.6 17.8 19.6 18.1 20.2 ...
## $ flipper_length_mm: int [1:344] 181 186 195 NA 193 190 181 195 193 190 ...
## $ body_mass_g : int [1:344] 3750 3800 3250 NA 3450 3650 3625 4675 3475 4250 ...
## $ sex : Factor w/ 2 levels "female","male": 2 1 1 NA 1 2 1 2 NA NA ...
## $ year : int [1:344] 2007 2007 2007 2007 2007 2007 2007 2007 2007 2007 ...glimpse(penguins)## Rows: 344
## Columns: 8
## $ species <fct> Adelie, Adelie, Adelie, Adelie, Adelie, Adelie, Adel…
## $ island <fct> Torgersen, Torgersen, Torgersen, Torgersen, Torgerse…
## $ bill_length_mm <dbl> 39.1, 39.5, 40.3, NA, 36.7, 39.3, 38.9, 39.2, 34.1, …
## $ bill_depth_mm <dbl> 18.7, 17.4, 18.0, NA, 19.3, 20.6, 17.8, 19.6, 18.1, …
## $ flipper_length_mm <int> 181, 186, 195, NA, 193, 190, 181, 195, 193, 190, 186…
## $ body_mass_g <int> 3750, 3800, 3250, NA, 3450, 3650, 3625, 4675, 3475, …
## $ sex <fct> male, female, female, NA, female, male, female, male…
## $ year <int> 2007, 2007, 2007, 2007, 2007, 2007, 2007, 2007, 2007…typeof(mtcars)## [1] "list"typeof(penguins)## [1] "list"is_tibble(mtcars)## [1] FALSEis_tibble(penguins)## [1] TRUEmtcars## mpg cyl disp hp drat wt qsec vs am gear carb
## Mazda RX4 21.0 6 160.0 110 3.90 2.620 16.46 0 1 4 4
## Mazda RX4 Wag 21.0 6 160.0 110 3.90 2.875 17.02 0 1 4 4
## Datsun 710 22.8 4 108.0 93 3.85 2.320 18.61 1 1 4 1
## Hornet 4 Drive 21.4 6 258.0 110 3.08 3.215 19.44 1 0 3 1
## Hornet Sportabout 18.7 8 360.0 175 3.15 3.440 17.02 0 0 3 2
## Valiant 18.1 6 225.0 105 2.76 3.460 20.22 1 0 3 1
## Duster 360 14.3 8 360.0 245 3.21 3.570 15.84 0 0 3 4
## Merc 240D 24.4 4 146.7 62 3.69 3.190 20.00 1 0 4 2
## Merc 230 22.8 4 140.8 95 3.92 3.150 22.90 1 0 4 2
## Merc 280 19.2 6 167.6 123 3.92 3.440 18.30 1 0 4 4
## Merc 280C 17.8 6 167.6 123 3.92 3.440 18.90 1 0 4 4
## Merc 450SE 16.4 8 275.8 180 3.07 4.070 17.40 0 0 3 3
## Merc 450SL 17.3 8 275.8 180 3.07 3.730 17.60 0 0 3 3
## Merc 450SLC 15.2 8 275.8 180 3.07 3.780 18.00 0 0 3 3
## Cadillac Fleetwood 10.4 8 472.0 205 2.93 5.250 17.98 0 0 3 4
## Lincoln Continental 10.4 8 460.0 215 3.00 5.424 17.82 0 0 3 4
## Chrysler Imperial 14.7 8 440.0 230 3.23 5.345 17.42 0 0 3 4
## Fiat 128 32.4 4 78.7 66 4.08 2.200 19.47 1 1 4 1
## Honda Civic 30.4 4 75.7 52 4.93 1.615 18.52 1 1 4 2
## Toyota Corolla 33.9 4 71.1 65 4.22 1.835 19.90 1 1 4 1
## Toyota Corona 21.5 4 120.1 97 3.70 2.465 20.01 1 0 3 1
## Dodge Challenger 15.5 8 318.0 150 2.76 3.520 16.87 0 0 3 2
## AMC Javelin 15.2 8 304.0 150 3.15 3.435 17.30 0 0 3 2
## Camaro Z28 13.3 8 350.0 245 3.73 3.840 15.41 0 0 3 4
## Pontiac Firebird 19.2 8 400.0 175 3.08 3.845 17.05 0 0 3 2
## Fiat X1-9 27.3 4 79.0 66 4.08 1.935 18.90 1 1 4 1
## Porsche 914-2 26.0 4 120.3 91 4.43 2.140 16.70 0 1 5 2
## Lotus Europa 30.4 4 95.1 113 3.77 1.513 16.90 1 1 5 2
## Ford Pantera L 15.8 8 351.0 264 4.22 3.170 14.50 0 1 5 4
## Ferrari Dino 19.7 6 145.0 175 3.62 2.770 15.50 0 1 5 6
## Maserati Bora 15.0 8 301.0 335 3.54 3.570 14.60 0 1 5 8
## Volvo 142E 21.4 4 121.0 109 4.11 2.780 18.60 1 1 4 2penguins## # A tibble: 344 × 8
## species island bill_length_mm bill_depth_mm flipper_length_mm body_mass_g
## <fct> <fct> <dbl> <dbl> <int> <int>
## 1 Adelie Torgersen 39.1 18.7 181 3750
## 2 Adelie Torgersen 39.5 17.4 186 3800
## 3 Adelie Torgersen 40.3 18 195 3250
## 4 Adelie Torgersen NA NA NA NA
## 5 Adelie Torgersen 36.7 19.3 193 3450
## 6 Adelie Torgersen 39.3 20.6 190 3650
## 7 Adelie Torgersen 38.9 17.8 181 3625
## 8 Adelie Torgersen 39.2 19.6 195 4675
## 9 Adelie Torgersen 34.1 18.1 193 3475
## 10 Adelie Torgersen 42 20.2 190 4250
## # ℹ 334 more rows
## # ℹ 2 more variables: sex <fct>, year <int>## tidy
filter(penguins,species=="Chinstrap")## # A tibble: 68 × 8
## species island bill_length_mm bill_depth_mm flipper_length_mm body_mass_g
## <fct> <fct> <dbl> <dbl> <int> <int>
## 1 Chinstrap Dream 46.5 17.9 192 3500
## 2 Chinstrap Dream 50 19.5 196 3900
## 3 Chinstrap Dream 51.3 19.2 193 3650
## 4 Chinstrap Dream 45.4 18.7 188 3525
## 5 Chinstrap Dream 52.7 19.8 197 3725
## 6 Chinstrap Dream 45.2 17.8 198 3950
## 7 Chinstrap Dream 46.1 18.2 178 3250
## 8 Chinstrap Dream 51.3 18.2 197 3750
## 9 Chinstrap Dream 46 18.9 195 4150
## 10 Chinstrap Dream 51.3 19.9 198 3700
## # ℹ 58 more rows
## # ℹ 2 more variables: sex <fct>, year <int>## base r
penguins[penguins$species=="Chinstrap",]## # A tibble: 68 × 8
## species island bill_length_mm bill_depth_mm flipper_length_mm body_mass_g
## <fct> <fct> <dbl> <dbl> <int> <int>
## 1 Chinstrap Dream 46.5 17.9 192 3500
## 2 Chinstrap Dream 50 19.5 196 3900
## 3 Chinstrap Dream 51.3 19.2 193 3650
## 4 Chinstrap Dream 45.4 18.7 188 3525
## 5 Chinstrap Dream 52.7 19.8 197 3725
## 6 Chinstrap Dream 45.2 17.8 198 3950
## 7 Chinstrap Dream 46.1 18.2 178 3250
## 8 Chinstrap Dream 51.3 18.2 197 3750
## 9 Chinstrap Dream 46 18.9 195 4150
## 10 Chinstrap Dream 51.3 19.9 198 3700
## # ℹ 58 more rows
## # ℹ 2 more variables: sex <fct>, year <int>penguins %>%
filter(species=="Chinstrap")## # A tibble: 68 × 8
## species island bill_length_mm bill_depth_mm flipper_length_mm body_mass_g
## <fct> <fct> <dbl> <dbl> <int> <int>
## 1 Chinstrap Dream 46.5 17.9 192 3500
## 2 Chinstrap Dream 50 19.5 196 3900
## 3 Chinstrap Dream 51.3 19.2 193 3650
## 4 Chinstrap Dream 45.4 18.7 188 3525
## 5 Chinstrap Dream 52.7 19.8 197 3725
## 6 Chinstrap Dream 45.2 17.8 198 3950
## 7 Chinstrap Dream 46.1 18.2 178 3250
## 8 Chinstrap Dream 51.3 18.2 197 3750
## 9 Chinstrap Dream 46 18.9 195 4150
## 10 Chinstrap Dream 51.3 19.9 198 3700
## # ℹ 58 more rows
## # ℹ 2 more variables: sex <fct>, year <int>filter(penguins,species=="Chinstrap",bill_length_mm>57)## # A tibble: 1 × 8
## species island bill_length_mm bill_depth_mm flipper_length_mm body_mass_g
## <fct> <fct> <dbl> <dbl> <int> <int>
## 1 Chinstrap Dream 58 17.8 181 3700
## # ℹ 2 more variables: sex <fct>, year <int>filter(penguins,
species %in% c("Chinstrap","Adelie"))## # A tibble: 220 × 8
## species island bill_length_mm bill_depth_mm flipper_length_mm body_mass_g
## <fct> <fct> <dbl> <dbl> <int> <int>
## 1 Adelie Torgersen 39.1 18.7 181 3750
## 2 Adelie Torgersen 39.5 17.4 186 3800
## 3 Adelie Torgersen 40.3 18 195 3250
## 4 Adelie Torgersen NA NA NA NA
## 5 Adelie Torgersen 36.7 19.3 193 3450
## 6 Adelie Torgersen 39.3 20.6 190 3650
## 7 Adelie Torgersen 38.9 17.8 181 3625
## 8 Adelie Torgersen 39.2 19.6 195 4675
## 9 Adelie Torgersen 34.1 18.1 193 3475
## 10 Adelie Torgersen 42 20.2 190 4250
## # ℹ 210 more rows
## # ℹ 2 more variables: sex <fct>, year <int>gapminder %>%
filter(year ==2007,lifeExp>80) %>%
select(country)## # A tibble: 13 × 1
## country
## <fct>
## 1 Australia
## 2 Canada
## 3 France
## 4 Hong Kong, China
## 5 Iceland
## 6 Israel
## 7 Italy
## 8 Japan
## 9 New Zealand
## 10 Norway
## 11 Spain
## 12 Sweden
## 13 Switzerlandgapminder %>%
filter(country %in% c("Bangladesh","India")) %>%
select(year,country,lifeExp) %>%
summarise(lifeExp)## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
## dplyr 1.1.0.
## ℹ Please use `reframe()` instead.
## ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
## always returns an ungrouped data frame and adjust accordingly.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## # A tibble: 24 × 1
## lifeExp
## <dbl>
## 1 37.5
## 2 39.3
## 3 41.2
## 4 43.5
## 5 45.3
## 6 46.9
## 7 50.0
## 8 52.8
## 9 56.0
## 10 59.4
## # ℹ 14 more rowsgapminder %>%
select()## # A tibble: 1,704 × 0penguins %>%
mutate(
mass_c= case_when(body_mass_g <= 3000 ~"small",
body_mass_g <= 4500 &
body_mass_g > 3000 ~"mediam",
body_mass_g > 4500 ~"lerge")
)## # A tibble: 344 × 9
## species island bill_length_mm bill_depth_mm flipper_length_mm body_mass_g
## <fct> <fct> <dbl> <dbl> <int> <int>
## 1 Adelie Torgersen 39.1 18.7 181 3750
## 2 Adelie Torgersen 39.5 17.4 186 3800
## 3 Adelie Torgersen 40.3 18 195 3250
## 4 Adelie Torgersen NA NA NA NA
## 5 Adelie Torgersen 36.7 19.3 193 3450
## 6 Adelie Torgersen 39.3 20.6 190 3650
## 7 Adelie Torgersen 38.9 17.8 181 3625
## 8 Adelie Torgersen 39.2 19.6 195 4675
## 9 Adelie Torgersen 34.1 18.1 193 3475
## 10 Adelie Torgersen 42 20.2 190 4250
## # ℹ 334 more rows
## # ℹ 3 more variables: sex <fct>, year <int>, mass_c <chr>penguins %>%
mutate(
mass_c= case_when(body_mass_g <= 3000 ~"small",
body_mass_g <= 4500 &
body_mass_g > 3000 ~"mediam",
body_mass_g > 4500 ~"lerge")
) %>%
count()## # A tibble: 1 × 1
## n
## <int>
## 1 344penguins %>%
group_by(sex) %>%
reframe(mean_mass=mean(body_mass_g,na.rm=T),
n=n())## # A tibble: 3 × 3
## sex mean_mass n
## <fct> <dbl> <int>
## 1 female 3862. 165
## 2 male 4546. 168
## 3 <NA> 4006. 11penguins %>%
group_by(sex) %>%
reframe(mean_mass=mean(body_mass_g,na.rm=T),
n=n()) %>%
reframe(mean(mean_mass))## # A tibble: 1 × 1
## `mean(mean_mass)`
## <dbl>
## 1 4138.penguins %>%
reframe(mean_mass=mean(body_mass_g,na.rm=T),
n=n(),.by=sex)## # A tibble: 3 × 3
## sex mean_mass n
## <fct> <dbl> <int>
## 1 male 4546. 168
## 2 female 3862. 165
## 3 <NA> 4006. 11penguins %>%
group_by(sex) %>%
reframe(mean_mass=mean(body_mass_g,na.rm=T),
n=n()) %>%
pull(mean_mass)## [1] 3862.273 4545.685 4005.556wdat<-penguins %>%
count(year,species) %>%
pivot_wider(names_from = "species",values_from = n)
wdat %>%
pivot_longer(cols = -year,names_to = "penguins",values_to = "feq")## # A tibble: 9 × 3
## year penguins feq
## <int> <chr> <int>
## 1 2007 Adelie 50
## 2 2007 Chinstrap 26
## 3 2007 Gentoo 34
## 4 2008 Adelie 50
## 5 2008 Chinstrap 18
## 6 2008 Gentoo 46
## 7 2009 Adelie 52
## 8 2009 Chinstrap 24
## 9 2009 Gentoo 44Chapter 1
Probability , Joint probability
dat_ex<-penguins %>%
count(year,species)
dat_ex %>%
mutate(joint= n/sum(n))## # A tibble: 9 × 4
## year species n joint
## <int> <fct> <int> <dbl>
## 1 2007 Adelie 50 0.145
## 2 2007 Chinstrap 26 0.0756
## 3 2007 Gentoo 34 0.0988
## 4 2008 Adelie 50 0.145
## 5 2008 Chinstrap 18 0.0523
## 6 2008 Gentoo 46 0.134
## 7 2009 Adelie 52 0.151
## 8 2009 Chinstrap 24 0.0698
## 9 2009 Gentoo 44 0.128dat_ex %>%
mutate(joint= n/sum(n)) %>%
mutate(year_1= n/sum(n),.by=year) %>% ##wrong
mutate(species_1= n/sum(n),.by=species) ##wrong## # A tibble: 9 × 6
## year species n joint year_1 species_1
## <int> <fct> <int> <dbl> <dbl> <dbl>
## 1 2007 Adelie 50 0.145 0.455 0.329
## 2 2007 Chinstrap 26 0.0756 0.236 0.382
## 3 2007 Gentoo 34 0.0988 0.309 0.274
## 4 2008 Adelie 50 0.145 0.439 0.329
## 5 2008 Chinstrap 18 0.0523 0.158 0.265
## 6 2008 Gentoo 46 0.134 0.404 0.371
## 7 2009 Adelie 52 0.151 0.433 0.342
## 8 2009 Chinstrap 24 0.0698 0.2 0.353
## 9 2009 Gentoo 44 0.128 0.367 0.355# correct
dat<-dat_ex %>%
mutate(n_s=sum(n),.by= species) %>%
mutate(n_y=sum(n),.by= year ) %>%
mutate(p_s=n_s/sum(n),p_y=n_y/sum(n),p_j=n/sum(n))
out<-dat %>%
select(species,year,p_j) %>%
pivot_wider(names_from = year,values_from = p_j)
out## # A tibble: 3 × 4
## species `2007` `2008` `2009`
## <fct> <dbl> <dbl> <dbl>
## 1 Adelie 0.145 0.145 0.151
## 2 Chinstrap 0.0756 0.0523 0.0698
## 3 Gentoo 0.0988 0.134 0.128dat %>%
filter(year==2007) %>%
select(species,p_s)## # A tibble: 3 × 2
## species p_s
## <fct> <dbl>
## 1 Adelie 0.442
## 2 Chinstrap 0.198
## 3 Gentoo 0.360dat %>%
filter(species=="Adelie") %>%
select(year,p_y)## # A tibble: 3 × 2
## year p_y
## <int> <dbl>
## 1 2007 0.320
## 2 2008 0.331
## 3 2009 0.349dat %>%
mutate(np= paste(n, "(",round(p_j,3),")") )## # A tibble: 9 × 9
## year species n n_s n_y p_s p_y p_j np
## <int> <fct> <int> <int> <int> <dbl> <dbl> <dbl> <chr>
## 1 2007 Adelie 50 152 110 0.442 0.320 0.145 50 ( 0.145 )
## 2 2007 Chinstrap 26 68 110 0.198 0.320 0.0756 26 ( 0.076 )
## 3 2007 Gentoo 34 124 110 0.360 0.320 0.0988 34 ( 0.099 )
## 4 2008 Adelie 50 152 114 0.442 0.331 0.145 50 ( 0.145 )
## 5 2008 Chinstrap 18 68 114 0.198 0.331 0.0523 18 ( 0.052 )
## 6 2008 Gentoo 46 124 114 0.360 0.331 0.134 46 ( 0.134 )
## 7 2009 Adelie 52 152 120 0.442 0.349 0.151 52 ( 0.151 )
## 8 2009 Chinstrap 24 68 120 0.198 0.349 0.0698 24 ( 0.07 )
## 9 2009 Gentoo 44 124 120 0.360 0.349 0.128 44 ( 0.128 )Treatment assignment mechanism (page 52)
\[ E[Y(0)] = -\frac{\phi(0.5)}{\Phi(0.5)} \]
\[ E[Y(1)] = \frac{2 \phi(0.5)}{1 - \Phi(0.5)} - 0.5 \]
library(tidyverse)
set.seed(1000)
n<- 1000
df<-tibble(id= 1:n,
y0= rnorm(n),
tau= -0.5+y0,
y1= tau+y0)
df## # A tibble: 1,000 × 4
## id y0 tau y1
## <int> <dbl> <dbl> <dbl>
## 1 1 -0.446 -0.946 -1.39
## 2 2 -1.21 -1.71 -2.91
## 3 3 0.0411 -0.459 -0.418
## 4 4 0.639 0.139 0.779
## 5 5 -0.787 -1.29 -2.07
## 6 6 -0.385 -0.885 -1.27
## 7 7 -0.476 -0.976 -1.45
## 8 8 0.720 0.220 0.940
## 9 9 -0.0185 -0.519 -0.537
## 10 10 -1.37 -1.87 -3.25
## # ℹ 990 more rows## true value
y0mean<- mean(df$y0-df$y1)
eya1 <- (2*dnorm(.5)/pnorm(0.5,lower.tail = F))-.5
eya0<- -dnorm(.5)/pnorm(0.5,lower.tail = T)
tau_true <- eya1 - eya0
tau_true## [1] 2.291316## prefect doctor
dat1 <- df %>%
mutate(A1 = if_else(tau > 0, 1,0),
Y = A1*y1 + (1-A1)*y0)
dat1## # A tibble: 1,000 × 6
## id y0 tau y1 A1 Y
## <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 -0.446 -0.946 -1.39 0 -0.446
## 2 2 -1.21 -1.71 -2.91 0 -1.21
## 3 3 0.0411 -0.459 -0.418 0 0.0411
## 4 4 0.639 0.139 0.779 1 0.779
## 5 5 -0.787 -1.29 -2.07 0 -0.787
## 6 6 -0.385 -0.885 -1.27 0 -0.385
## 7 7 -0.476 -0.976 -1.45 0 -0.476
## 8 8 0.720 0.220 0.940 1 0.940
## 9 9 -0.0185 -0.519 -0.537 0 -0.0185
## 10 10 -1.37 -1.87 -3.25 0 -1.37
## # ℹ 990 more rowsdat1 %>%
reframe(EY1 = mean(Y), .by = A1) %>%
mutate(tau <- EY1-EY1[1])## # A tibble: 2 × 3
## A1 EY1 `tau <- EY1 - EY1[1]`
## <dbl> <dbl> <dbl>
## 1 0 -0.501 0
## 2 1 1.70 2.20## clueless doctor
dat2<- df %>%
mutate(A2= rbinom(n = nrow(df),size = 1,prob = 0.5),Y= A2*y1+(1-A2)*y0)
dat2## # A tibble: 1,000 × 6
## id y0 tau y1 A2 Y
## <int> <dbl> <dbl> <dbl> <int> <dbl>
## 1 1 -0.446 -0.946 -1.39 1 -1.39
## 2 2 -1.21 -1.71 -2.91 0 -1.21
## 3 3 0.0411 -0.459 -0.418 0 0.0411
## 4 4 0.639 0.139 0.779 0 0.639
## 5 5 -0.787 -1.29 -2.07 1 -2.07
## 6 6 -0.385 -0.885 -1.27 1 -1.27
## 7 7 -0.476 -0.976 -1.45 1 -1.45
## 8 8 0.720 0.220 0.940 1 0.940
## 9 9 -0.0185 -0.519 -0.537 0 -0.0185
## 10 10 -1.37 -1.87 -3.25 1 -3.25
## # ℹ 990 more rowsdat2 %>%
reframe(EY1 = mean(Y), .by=A2) %>%
mutate(tau= EY1-EY1[2])## # A tibble: 2 × 3
## A2 EY1 tau
## <int> <dbl> <dbl>
## 1 1 -0.524 -0.512
## 2 0 -0.0128 0Abir code
set.seed(1000)
y0 <- rnorm(1000)
tao <- -0.5 + y0
y1 <- tao + y0
library(tidyverse)
dat <- data.frame(c(1:1000) , y0, tao, y1)
dat8 <- tibble(id = 1:1000, y0, tao, y1)
mean(y1) - mean(y0)## [1] -0.5124759ey1 <- (2*dnorm(0.5)/(1- pnorm(0.5))) - 0.5
ey0 <- -dnorm(0.5)/pnorm(0.5)
tau_true <- ey1 - ey0
tau_true## [1] 2.291316dat8 <- dat8 %>%
mutate(A1 = if_else(tao > 0, 1,0),
Y = A1*y1 + (1-A1)*y0)
dat8 %>%
reframe(EY1 = mean(Y),
.by = A1) %>%
mutate(tau <- EY1-EY1[1])## # A tibble: 2 × 3
## A1 EY1 `tau <- EY1 - EY1[1]`
## <dbl> <dbl> <dbl>
## 1 0 -0.501 0
## 2 1 1.70 2.20EY0<-mean(dat8$Y[dat8$A1 == 0])
EY1<-mean(dat8$Y[dat8$A1 == 1])
t<-EY1-EY0
t## [1] 2.200073Assignment 01
Q1
library(Matching) # contains lalonde## Warning: package 'Matching' was built under R version 4.4.3## Loading required package: MASS##
## Attaching package: 'MASS'## The following object is masked from 'package:dplyr':
##
## select## ##
## ## Matching (Version 4.10-15, Build Date: 2024-10-14)
## ## See https://www.jsekhon.com for additional documentation.
## ## Please cite software as:
## ## Jasjeet S. Sekhon. 2011. ``Multivariate and Propensity Score Matching
## ## Software with Automated Balance Optimization: The Matching package for R.''
## ## Journal of Statistical Software, 42(7): 1-52.
## ##data(lalonde)
# Look at structure
str(lalonde)## 'data.frame': 445 obs. of 12 variables:
## $ age : int 37 22 30 27 33 22 23 32 22 33 ...
## $ educ : int 11 9 12 11 8 9 12 11 16 12 ...
## $ black : int 1 0 1 1 1 1 1 1 1 0 ...
## $ hisp : int 0 1 0 0 0 0 0 0 0 0 ...
## $ married: int 1 0 0 0 0 0 0 0 0 1 ...
## $ nodegr : int 1 1 0 1 1 1 0 1 0 0 ...
## $ re74 : num 0 0 0 0 0 0 0 0 0 0 ...
## $ re75 : num 0 0 0 0 0 0 0 0 0 0 ...
## $ re78 : num 9930 3596 24910 7506 290 ...
## $ u74 : int 1 1 1 1 1 1 1 1 1 1 ...
## $ u75 : int 1 1 1 1 1 1 1 1 1 1 ...
## $ treat : int 1 1 1 1 1 1 1 1 1 1 ...# Fit linear regression of re78 on treat + covariates (associational adjustment)
fit.lm <- lm(re78 ~ factor(treat) + age + educ + factor(black) + factor(hisp )+ factor(married) + factor(nodegr) + factor(u74) + factor(u75),
data = lalonde)
summary(fit.lm)##
## Call:
## lm(formula = re78 ~ factor(treat) + age + educ + factor(black) +
## factor(hisp) + factor(married) + factor(nodegr) + factor(u74) +
## factor(u75), data = lalonde)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9703 -4478 -1537 3066 53184
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1720.83 3414.38 0.504 0.6145
## factor(treat)1 1600.87 640.02 2.501 0.0127 *
## age 57.81 45.59 1.268 0.2055
## educ 367.43 228.05 1.611 0.1079
## factor(black)1 -2103.41 1173.28 -1.793 0.0737 .
## factor(hisp)1 218.67 1559.72 0.140 0.8886
## factor(married)1 43.61 855.51 0.051 0.9594
## factor(nodegr)1 -135.06 995.50 -0.136 0.8921
## factor(u74)1 380.77 1025.37 0.371 0.7106
## factor(u75)1 -1146.33 955.70 -1.199 0.2310
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6523 on 435 degrees of freedom
## Multiple R-squared: 0.05213, Adjusted R-squared: 0.03252
## F-statistic: 2.658 on 9 and 435 DF, p-value: 0.005212Q_2a
library(tidyverse)
library(dplyr)
d<-datasets::UCBAdmissions %>% as.tibble()## Warning: `as.tibble()` was deprecated in tibble 2.0.0.
## ℹ Please use `as_tibble()` instead.
## ℹ The signature and semantics have changed, see `?as_tibble`.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.d## # A tibble: 24 × 4
## Admit Gender Dept n
## <chr> <chr> <chr> <dbl>
## 1 Admitted Male A 512
## 2 Rejected Male A 313
## 3 Admitted Female A 89
## 4 Rejected Female A 19
## 5 Admitted Male B 353
## 6 Rejected Male B 207
## 7 Admitted Female B 17
## 8 Rejected Female B 8
## 9 Admitted Male C 120
## 10 Rejected Male C 205
## # ℹ 14 more rowsd1<-d %>%
reframe(n=sum(n),.by= c(Gender,Admit)) %>%
mutate(n_g= sum(n),.by= Gender) %>%
mutate(n_a= sum(n),.by= Admit) %>%
mutate(p_g = n_g / sum(n), # Proportion by gender
p_a = n_a / sum(n), # Proportion by admission status
p_j = n / sum(n) # Joint proportion
)
d1## # A tibble: 4 × 8
## Gender Admit n n_g n_a p_g p_a p_j
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Male Admitted 1198 2691 1755 0.595 0.388 0.265
## 2 Male Rejected 1493 2691 2771 0.595 0.612 0.330
## 3 Female Admitted 557 1835 1755 0.405 0.388 0.123
## 4 Female Rejected 1278 1835 2771 0.405 0.612 0.282result<-function(data){
p_ma<-(data$p_j[1]/(data$p_j[1]+data$p_j[2]))
p_fe<- (data$p_j[3]/(data$p_j[3]+data$p_j[4]))
RD<- p_ma-p_fe
RR<- p_ma/p_fe
OR<- (p_ma/(1-p_ma))/(p_fe/(1-p_fe))
var_Or<- ((p_ma*(1-p_ma))/(data$n[2]+data$n[4]))+((p_fe*(1-p_fe))/(data$n[1]+data$n[3]))
# Output results
cat("RD:", RD, "\n")
cat("RR", RR, "\n")
cat("OR", OR, "\n")
cat("Variance-OR", var_Or, "\n")
}
result(d1)## RD: 0.1416454
## RR 1.466642
## OR 1.84108
## Variance-OR 0.0002095942Q_2b
fun_data<-function(dept){
d %>%
filter(Dept==dept) %>%
reframe(n=sum(n),.by= c(Gender,Admit)) %>%
mutate(n_g= sum(n),.by= Gender) %>%
mutate(n_a= sum(n),.by= Admit) %>%
mutate(p_g = n_g / sum(n), # Proportion by gender
p_a = n_a / sum(n), # Proportion by admission status
p_j = n / sum(n) # Joint proportion
)
}
d1<-fun_data("A")
d2<-fun_data("B")
d3<-fun_data("C")
d4<-fun_data("D")
d5<-fun_data("E")
d6<-fun_data("F")
result(d1)## RD: -0.203468
## RR 0.753095
## OR 0.349212
## Variance-OR 0.0009504239result(d2)## RD: -0.04964286
## RR 0.9269958
## OR 0.8025007
## Variance-OR 0.001671862result(d3)## RD: 0.02858996
## RR 1.08393
## OR 1.13306
## Variance-OR 0.001088301result(d4)## RD: -0.01839808
## RR 0.9473337
## OR 0.9212838
## Variance-OR 0.001268339result(d5)## RD: 0.03830116
## RR 1.160131
## OR 1.221631
## Variance-OR 0.001696714result(d6)## RD: -0.0114
## RR 0.838025
## OR 0.8278727
## Variance-OR 0.001505429Q_3a
library(tidyverse)
set.seed(1000)
n<- 1000
df<-tibble(id= 1:n,
y0= rnorm(n),
tau= -0.5+y0,
y1= tau+y0)
df## # A tibble: 1,000 × 4
## id y0 tau y1
## <int> <dbl> <dbl> <dbl>
## 1 1 -0.446 -0.946 -1.39
## 2 2 -1.21 -1.71 -2.91
## 3 3 0.0411 -0.459 -0.418
## 4 4 0.639 0.139 0.779
## 5 5 -0.787 -1.29 -2.07
## 6 6 -0.385 -0.885 -1.27
## 7 7 -0.476 -0.976 -1.45
## 8 8 0.720 0.220 0.940
## 9 9 -0.0185 -0.519 -0.537
## 10 10 -1.37 -1.87 -3.25
## # ℹ 990 more rowsdf1<-df %>%
mutate(A1=if_else(tau>0,1,0)) %>%
mutate(Y=A1*y1+(1-A1)*y0)
df1 %>%
reframe(Ey1=mean(Y), .by=A1) %>%
mutate(tau<- Ey1-Ey1[1])## # A tibble: 2 × 3
## A1 Ey1 `tau <- Ey1 - Ey1[1]`
## <dbl> <dbl> <dbl>
## 1 0 -0.501 0
## 2 1 1.70 2.20Q_3b
df2<-df1 %>%
mutate(A2=rbinom(1000,size = 1,prob = 0.5)) %>%
mutate(Y_new=A2*y1+(1-A2)*y0)
df2## # A tibble: 1,000 × 8
## id y0 tau y1 A1 Y A2 Y_new
## <int> <dbl> <dbl> <dbl> <dbl> <dbl> <int> <dbl>
## 1 1 -0.446 -0.946 -1.39 0 -0.446 1 -1.39
## 2 2 -1.21 -1.71 -2.91 0 -1.21 0 -1.21
## 3 3 0.0411 -0.459 -0.418 0 0.0411 0 0.0411
## 4 4 0.639 0.139 0.779 1 0.779 0 0.639
## 5 5 -0.787 -1.29 -2.07 0 -0.787 1 -2.07
## 6 6 -0.385 -0.885 -1.27 0 -0.385 1 -1.27
## 7 7 -0.476 -0.976 -1.45 0 -0.476 1 -1.45
## 8 8 0.720 0.220 0.940 1 0.940 1 0.940
## 9 9 -0.0185 -0.519 -0.537 0 -0.0185 0 -0.0185
## 10 10 -1.37 -1.87 -3.25 0 -1.37 1 -3.25
## # ℹ 990 more rowsdf2 %>%
reframe(Ey2=mean(Y_new), .by=A2) %>%
mutate(tau<- Ey2-Ey2[2])## # A tibble: 2 × 3
## A2 Ey2 `tau <- Ey2 - Ey2[2]`
## <int> <dbl> <dbl>
## 1 1 -0.524 -0.512
## 2 0 -0.0128 0Q4
d<-datasets::UCBAdmissions %>% as.tibble()
d## # A tibble: 24 × 4
## Admit Gender Dept n
## <chr> <chr> <chr> <dbl>
## 1 Admitted Male A 512
## 2 Rejected Male A 313
## 3 Admitted Female A 89
## 4 Rejected Female A 19
## 5 Admitted Male B 353
## 6 Rejected Male B 207
## 7 Admitted Female B 17
## 8 Rejected Female B 8
## 9 Admitted Male C 120
## 10 Rejected Male C 205
## # ℹ 14 more rowsresult1<-function(data){
data<-data %>%
reframe(n=sum(n),.by= c(Gender,Admit)) %>%
mutate(n_g= sum(n),.by= Gender) %>%
mutate(n_a= sum(n),.by= Admit) %>%
mutate(p_g = n_g / sum(n), # Proportion by gender
p_a = n_a / sum(n), # Proportion by admission status
p_j = n / sum(n) # Joint proportion
)
p_ma<-(data$p_j[1]/(data$p_j[1]+data$p_j[2]))
p_fe<- (data$p_j[3]/(data$p_j[3]+data$p_j[4]))
RD<- p_ma-p_fe
RR<- p_ma/p_fe
OR<- (p_ma/(1-p_ma))/(p_fe/(1-p_fe))
# Output results
cat("RD:", RD, "\n")
cat("RR", RR, "\n")
cat("OR", OR, "\n")
}
result1(d)## RD: 0.1416454
## RR 1.466642
## OR 1.84108## for group by department
fun_data<-function(dept){
d %>%
filter(Dept==dept) %>%
reframe(n=sum(n),.by= c(Gender,Admit)) %>%
mutate(n_g= sum(n),.by= Gender) %>%
mutate(n_a= sum(n),.by= Admit) %>%
mutate(p_g = n_g / sum(n), # Proportion by gender
p_a = n_a / sum(n), # Proportion by admission status
p_j = n / sum(n) # Joint proportion
)
}
d1<-fun_data("A")
d2<-fun_data("B")
d3<-fun_data("C")
d4<-fun_data("D")
d5<-fun_data("E")
d6<-fun_data("F")
result1(d1)## RD: -0.203468
## RR 0.753095
## OR 0.349212result1(d2)## RD: -0.04964286
## RR 0.9269958
## OR 0.8025007result1(d3)## RD: 0.02858996
## RR 1.08393
## OR 1.13306result1(d4)## RD: -0.01839808
## RR 0.9473337
## OR 0.9212838result1(d5)## RD: 0.03830116
## RR 1.160131
## OR 1.221631result1(d6)## RD: -0.0114
## RR 0.838025
## OR 0.8278727practice
d<-datasets::UCBAdmissions %>% as.tibble()
d## # A tibble: 24 × 4
## Admit Gender Dept n
## <chr> <chr> <chr> <dbl>
## 1 Admitted Male A 512
## 2 Rejected Male A 313
## 3 Admitted Female A 89
## 4 Rejected Female A 19
## 5 Admitted Male B 353
## 6 Rejected Male B 207
## 7 Admitted Female B 17
## 8 Rejected Female B 8
## 9 Admitted Male C 120
## 10 Rejected Male C 205
## # ℹ 14 more rowsd %>% count(Dept)## # A tibble: 6 × 2
## Dept n
## <chr> <int>
## 1 A 4
## 2 B 4
## 3 C 4
## 4 D 4
## 5 E 4
## 6 F 4sum(d$n)## [1] 4526d %>%
reframe(n=sum(n))## # A tibble: 1 × 1
## n
## <dbl>
## 1 4526d %>%
reframe(n=sum(n),.by= Gender)## # A tibble: 2 × 2
## Gender n
## <chr> <dbl>
## 1 Male 2691
## 2 Female 1835d %>%
reframe(n=sum(n),.by= Admit)## # A tibble: 2 × 2
## Admit n
## <chr> <dbl>
## 1 Admitted 1755
## 2 Rejected 2771d %>%
reframe(n=sum(n),.by= c(Gender,Admit))## # A tibble: 4 × 3
## Gender Admit n
## <chr> <chr> <dbl>
## 1 Male Admitted 1198
## 2 Male Rejected 1493
## 3 Female Admitted 557
## 4 Female Rejected 1278d %>%
reframe(n=sum(n),.by= c(Gender,Admit))%>%
mutate(total= sum(n),.by= Gender)## # A tibble: 4 × 4
## Gender Admit n total
## <chr> <chr> <dbl> <dbl>
## 1 Male Admitted 1198 2691
## 2 Male Rejected 1493 2691
## 3 Female Admitted 557 1835
## 4 Female Rejected 1278 1835Chapter 2
Lecture (page 57)
library(Matching)
data("lalonde")
FRT_f <- function(Y, A) {
y1 <- Y[A == 1]
y0 <- Y[A == 0]
selc <- c("statistic", "p.value")
m_d <- broom::tidy(t.test(x = y1, y = y0, var.equal = T))[selc]
m_t <- broom::tidy(t.test(x = y1, y = y0, var.equal = F))[selc]
m_W <- broom::tidy(wilcox.test(x = y1, y = y0))[selc]
#
return(tibble(method = c("D", "S", "W"), bind_rows(m_d, m_t,
m_W)))
}
out_appx <- FRT_f(A = lalonde$treat, Y = lalonde$re78)
out_appx## # A tibble: 3 × 3
## method statistic p.value
## <chr> <dbl> <dbl>
## 1 D 2.84 0.00479
## 2 S 2.67 0.00789
## 3 W 27402. 0.0109res <- list()
#
for (i in 1:100) {
a_perm <- sample(lalonde$treat, replace = F)
res[[i]] <- FRT_f(Y = lalonde$re78, A = a_perm)
}
out <- bind_rows(res) |>
left_join(out_appx |>
dplyr :: select(method, est = statistic))## Joining with `by = join_by(method)`out## # A tibble: 300 × 4
## method statistic p.value est
## <chr> <dbl> <dbl> <dbl>
## 1 D -0.651 0.515 2.84
## 2 S -0.636 0.525 2.67
## 3 W 22008. 0.121 27402.
## 4 D 0.0637 0.949 2.84
## 5 S 0.0659 0.947 2.67
## 6 W 25078. 0.436 27402.
## 7 D -0.745 0.457 2.84
## 8 S -0.765 0.445 2.67
## 9 W 23103 0.472 27402.
## 10 D 0.407 0.684 2.84
## # ℹ 290 more rows #save(out, file = here::here("data/FRT-perm.Rdata"))Lalonde data (permuntation test)
library(Matching)
library(tidyverse)
library(broom)
data("lalonde")
lalonde %>%
count(treat)## treat n
## 1 0 260
## 2 1 185## 1st way for t-test
t1 <- t.test(re78 ~ treat, var.equal = T, data = lalonde)
t1##
## Two Sample t-test
##
## data: re78 by treat
## t = -2.8353, df = 443, p-value = 0.004788
## alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
## 95 percent confidence interval:
## -3038.1113 -550.5749
## sample estimates:
## mean in group 0 mean in group 1
## 4554.802 6349.145## 2nd way for t-test
t1a <- t.test(lalonde$re78[lalonde$treat == 1],lalonde$re78[lalonde$treat == 0],
var.equal = T)
t1a##
## Two Sample t-test
##
## data: lalonde$re78[lalonde$treat == 1] and lalonde$re78[lalonde$treat == 0]
## t = 2.8353, df = 443, p-value = 0.004788
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 550.5749 3038.1113
## sample estimates:
## mean of x mean of y
## 6349.145 4554.802names(t1)## [1] "statistic" "parameter" "p.value" "conf.int" "estimate"
## [6] "null.value" "stderr" "alternative" "method" "data.name"res <- list()
n <- dim(lalonde)[1]
for (i in 1:1000) {
A <- sample(x = lalonde$treat, size = n, replace = F)
Y <- lalonde$re78
res[[i]] <- t.test(Y ~ A, var.equal = T)$statistic
}
tibble(stat = unlist(res)) %>%
reframe(p = mean(abs(stat) > abs(t1$statistic)))## # A tibble: 1 × 1
## p
## <dbl>
## 1 0.003fun_FRT <- function(Y, A, nsim = 1000) {
#
res <- list()
n <- length(A)
t1 <- t.test(Y ~ A, var.equal = T)
for (i in 1:nsim) {
A <- sample(x = A, size = n, replace = F)
res[[i]] <- t.test(Y ~ A, var.equal = T)$statistic
}
out <- tibble(stat = unlist(res))
p <- out %>%
reframe(p = mean(abs(stat) > abs(t1$statistic)))
## p <- mean(abs(out$stat) > abs(t1$stat))
##
return(list(out, p))
}
res <- fun_FRT(Y =lalonde$re78, A=lalonde$treat)
res## [[1]]
## # A tibble: 1,000 × 1
## stat
## <dbl>
## 1 -0.228
## 2 0.837
## 3 0.950
## 4 -1.35
## 5 0.334
## 6 0.332
## 7 -0.502
## 8 0.206
## 9 -0.327
## 10 1.45
## # ℹ 990 more rows
##
## [[2]]
## # A tibble: 1 × 1
## p
## <dbl>
## 1 0.004res[[1]] %>%
ggplot() +
geom_histogram(aes(stat, after_stat(density))) +
geom_density(aes(stat), col = "red", bw = .75)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## 2nd way
res[[1]] %>%
ggplot(aes(stat))+
geom_histogram(aes(stat, after_stat(density)))+
geom_density(col = "red", bw = .75)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Fisher Randomized test (Approch I)
mod <- lm(re78 ~ . - treat, data = lalonde)
mod##
## Call:
## lm(formula = re78 ~ . - treat, data = lalonde)
##
## Coefficients:
## (Intercept) age educ black hisp married
## 1.738e+03 5.681e+01 3.677e+02 -2.126e+03 1.051e+02 -6.883e+01
## nodegr re74 re75 u74 u75
## -3.762e+02 1.068e-01 3.048e-02 1.309e+03 -1.205e+03## 1st way
Yhat <- residuals(mod, newdata = lalonde)
Yhat## 1 2 3 4 5 6
## 4514.01596 -2532.72942 19078.41076 2589.39797 -3864.77974 159.15601
## 7 8 9 10 11 12
## -5433.41417 3271.35435 -4683.31230 4359.26090 4447.00818 11407.12506
## 13 14 15 16 17 18
## -3302.40888 14259.65541 145.96405 -866.08403 8553.56072 3371.57308
## 19 20 21 22 23 24
## 4405.10352 5143.55365 -4689.50914 -4168.01455 -523.49849 -2754.91986
## 25 26 27 28 29 30
## 7125.42670 4681.45175 5093.91942 4362.54475 12237.03742 -3750.61986
## 31 32 33 34 35 36
## -3980.96258 3015.58825 -3159.88535 -5262.98200 -2355.28562 13257.71765
## 37 38 39 40 41 42
## 3584.84269 546.96106 -1730.12392 -4378.63765 -4065.21273 7335.75527
## 43 44 45 46 47 48
## 4815.89184 -2676.90908 -546.53696 -4734.30535 -3245.59816 -3613.28037
## 49 50 51 52 53 54
## 9590.36269 -590.23417 -3870.51612 3012.32031 -4405.45552 -3198.06562
## 55 56 57 58 59 60
## -4427.57697 -3980.96258 -3681.98485 4215.48855 5988.14288 -5710.41133
## 61 62 63 64 65 66
## -3245.59816 -2003.38684 -337.31561 -4367.45403 -1693.68300 6456.10167
## 67 68 69 70 71 72
## -4916.75203 -5344.92321 4882.29344 1337.69927 -4181.38761 -6460.53804
## 73 74 75 76 77 78
## -2567.04446 6418.76438 2778.02557 11854.12650 13301.49655 -3840.52326
## 79 80 81 82 83 84
## 17350.30856 1824.00159 -5547.03562 -2828.25055 -3469.86795 -6028.33928
## 85 86 87 88 89 90
## -3813.70540 -5541.66999 -1936.46793 12371.05869 602.66436 -3700.08395
## 91 92 93 94 95 96
## 9434.73031 -3980.96258 -4803.13058 4663.98727 -5882.53434 -2085.67854
## 97 98 99 100 101 102
## 812.75891 -4056.66793 -4435.44837 20177.68406 -4378.63765 4773.53090
## 103 104 105 106 107 108
## -4317.90058 409.17448 2995.84376 -5090.35919 28070.96072 -1660.01037
## 109 110 111 112 113 114
## -5087.18420 -3331.45853 5756.62463 419.49414 -5270.48824 -3470.58554
## 115 116 117 118 119 120
## -5890.96386 4812.93233 16365.96357 -4827.26518 2812.97796 -5262.38774
## 121 122 123 124 125 126
## 10079.22663 4469.61234 12728.53383 -9010.71072 -2527.65134 -3573.03750
## 127 128 129 130 131 132
## -5900.80371 5746.34241 -8038.79317 -3756.23080 1713.55480 54089.78627
## 133 134 135 136 137 138
## -4578.62204 738.87565 -680.41966 4697.48486 -2669.09795 6923.47729
## 139 140 141 142 143 144
## -6888.98141 5126.54961 -105.94557 -7058.48820 -3142.80492 -5207.16737
## 145 146 147 148 149 150
## -6814.49119 -538.40853 5325.06502 -5130.56923 20291.81429 947.59296
## 151 152 153 154 155 156
## -2701.17503 -776.59325 3908.76963 2391.62670 -4968.94424 -7226.44944
## 157 158 159 160 161 162
## 34.64323 3361.30515 -4661.00640 691.65758 -682.76924 3882.98056
## 163 164 165 166 167 168
## -5581.16084 -6002.17434 4115.30294 6241.19382 -479.08036 3940.00094
## 169 170 171 172 173 174
## 8808.13077 -3854.34718 1330.49195 -5081.79839 -6042.03997 -4616.31681
## 175 176 177 178 179 180
## -6693.18873 -5285.86440 2389.32743 -421.39347 4211.28779 -3257.21656
## 181 182 183 184 185 186
## 7623.10986 26441.64522 6465.48267 -2324.40691 -3236.12462 -4321.82692
## 187 188 189 190 191 192
## 4653.71074 -3897.33399 7070.00891 5857.15494 5556.95891 4799.88386
## 193 194 195 196 197 198
## 794.73813 -522.99981 -5060.36634 -280.09715 -4238.19833 -3377.19795
## 199 200 201 202 203 204
## -6588.40916 -3416.03033 -3944.61214 3240.00235 4160.01279 -606.75058
## 205 206 207 208 209 210
## -4863.11630 13000.22342 4574.83009 -30.04540 -184.15395 12238.28728
## 211 212 213 214 215 216
## -4877.91965 -5336.44427 -31.26141 -4409.88874 656.91889 1396.49003
## 217 218 219 220 221 222
## -4557.72461 -5262.98200 -4803.13058 2931.59163 -5371.23782 -2135.89475
## 223 224 225 226 227 228
## 9492.04841 -3413.59982 -371.64240 2316.24869 -4438.62336 3005.79511
## 229 230 231 232 233 234
## -5206.17128 3054.37098 2314.82203 -4321.82692 7255.91222 -4746.31986
## 235 236 237 238 239 240
## 399.25086 -905.59355 -3022.59769 -4435.44837 -2865.57318 -4689.50914
## 241 242 243 244 245 246
## 3678.23218 -5030.37348 9880.08014 9001.06601 793.56724 258.86126
## 247 248 249 250 251 252
## -561.01910 16506.75163 -4378.63765 -2653.84804 9630.03366 -4632.69841
## 253 254 255 256 257 258
## -2128.93562 6280.54797 -3927.32684 -5490.22490 -4265.01620 -1819.90101
## 259 260 261 262 263 264
## -3883.33852 -2479.97078 7782.31076 -2139.14178 2277.80539 4447.12525
## 265 266 267 268 269 270
## -2394.22710 -5376.60345 -1371.85089 -4776.31272 -3670.09109 -8087.00190
## 271 272 273 274 275 276
## -4689.54985 -4253.00169 -5173.98778 -664.81314 1529.10932 -2477.76052
## 277 278 279 280 281 282
## -1163.93912 -1082.06101 5919.69435 1581.18072 -2323.22910 -5882.53434
## 283 284 285 286 287 288
## -3538.90058 -5547.03562 -5648.64256 -4833.12344 -4265.01620 1985.73818
## 289 290 291 292 293 294
## -4707.48748 -4689.50914 -5257.61637 -3364.18982 5264.35175 -4321.82692
## 295 296 297 298 299 300
## 716.42907 -4145.65837 -4378.63765 2193.76096 1791.51975 847.05157
## 301 302 303 304 305 306
## -4366.62314 -4961.54825 -5478.21039 -6085.15001 -338.59914 1361.87842
## 307 308 309 310 311 312
## -5282.80854 228.62967 -359.18064 -3760.06966 -5603.84635 -4321.82692
## 313 314 315 316 317 318
## 7006.10001 -6684.71215 -4366.62314 160.23750 -1734.25888 3764.84184
## 319 320 321 322 323 324
## -2146.50725 5264.90569 251.81448 4202.88521 215.29597 1191.81597
## 325 326 327 328 329 330
## 8816.10891 7207.30891 -3835.48830 -1323.25258 1197.85333 6897.64448
## 331 332 333 334 335 336
## -6325.86946 -3670.09109 -3980.96258 -1791.53258 -4462.26624 4839.25893
## 337 338 339 340 341 342
## 7010.32670 -3670.09109 5748.27184 216.65939 -3026.71695 -6474.78771
## 343 344 345 346 347 348
## 1409.37746 -3613.28037 -3705.39658 -22.58037 -4462.26624 8771.94121
## 349 350 351 352 353 354
## -1747.80254 -1280.41109 -4405.45552 2257.03199 -4462.26624 -2837.57258
## 355 356 357 358 359 360
## -729.75580 -3256.54930 748.13586 6706.29204 -3327.41081 13461.37441
## 361 362 363 364 365 366
## 573.20596 -4058.70537 5604.16083 -6806.44732 -7381.70941 -6700.51371
## 367 368 369 370 371 372
## -105.02103 -1152.81112 262.84362 187.84289 -3810.20309 -9234.53414
## 373 374 375 376 377 378
## -4513.98480 8530.30963 -4993.23561 281.79516 -5285.82875 -2326.51477
## 379 380 381 382 383 384
## -5318.41286 1061.93877 -5206.63268 -6820.17846 -5164.35197 7403.48977
## 385 386 387 388 389 390
## -5156.47440 34631.47103 4637.61917 1377.37785 3841.53606 -1186.40800
## 391 392 393 394 395 396
## -2775.95445 3503.69092 15631.95367 -6397.61236 3262.48757 -4490.57655
## 397 398 399 400 401 402
## -5241.56754 -6985.58986 38.16774 2594.70641 -379.16025 -5362.69260
## 403 404 405 406 407 408
## 8793.42898 -192.69858 -3459.20673 5261.63126 -4897.79257 -4199.68041
## 409 410 411 412 413 414
## -4367.01800 -3306.05516 12492.84451 -6906.52083 -6133.81604 6362.37062
## 415 416 417 418 419 420
## -4412.16151 -4797.07386 1817.23406 -6124.07035 -6190.24210 -5173.33694
## 421 422 423 424 425 426
## -909.34539 5934.55534 -419.95587 -7264.41145 4700.51018 -1330.37195
## 427 428 429 430 431 432
## 11471.80504 -7986.18281 -4714.85245 -1974.52058 2154.06110 13562.27414
## 433 434 435 436 437 438
## -1367.65961 -5461.20233 -4199.10263 -2524.77555 -5991.80246 -5195.39823
## 439 440 441 442 443 444
## -2492.84267 5743.43202 -7520.39769 -7146.51284 8954.21155 289.63240
## 445
## -4220.07425## 2nd and best way
estData <- broom::augment(mod, newdata = lalonde)
estData## # A tibble: 445 × 14
## age educ black hisp married nodegr re74 re75 re78 u74 u75 treat
## <int> <int> <int> <int> <int> <int> <dbl> <dbl> <dbl> <int> <int> <int>
## 1 37 11 1 0 1 1 0 0 9930. 1 1 1
## 2 22 9 0 1 0 1 0 0 3596. 1 1 1
## 3 30 12 1 0 0 0 0 0 24910. 1 1 1
## 4 27 11 1 0 0 1 0 0 7506. 1 1 1
## 5 33 8 1 0 0 1 0 0 290. 1 1 1
## 6 22 9 1 0 0 1 0 0 4056. 1 1 1
## 7 23 12 1 0 0 0 0 0 0 1 1 1
## 8 32 11 1 0 0 1 0 0 8472. 1 1 1
## 9 22 16 1 0 0 0 0 0 2164. 1 1 1
## 10 33 12 0 0 1 0 0 0 12418. 1 1 1
## # ℹ 435 more rows
## # ℹ 2 more variables: .fitted <dbl>, .resid <dbl>res1 <- fun_FRT(Y = estData$.resid, A = estData$treat)
res1## [[1]]
## # A tibble: 1,000 × 1
## stat
## <dbl>
## 1 1.24
## 2 0.0172
## 3 -1.13
## 4 0.963
## 5 -0.184
## 6 0.455
## 7 -0.827
## 8 -0.762
## 9 0.563
## 10 2.44
## # ℹ 990 more rows
##
## [[2]]
## # A tibble: 1 × 1
## p
## <dbl>
## 1 0.007res1 <- list()
set.seed(100)
for(i in 1:500) {
lalonde_n <- lalonde %>%
mutate(A = sample(treat, size = n(), replace = F))
mod_s <- lm(re78 ~ . - treat, data = lalonde_n)
res1[[i]] <- coefficients(mod_s)["A"]
}
res1## [[1]]
## A
## 66.05447
##
## [[2]]
## A
## -113.6653
##
## [[3]]
## A
## -91.76918
##
## [[4]]
## A
## -213.1234
##
## [[5]]
## A
## 140.9742
##
## [[6]]
## A
## 471.3521
##
## [[7]]
## A
## -605.5254
##
## [[8]]
## A
## -659.8794
##
## [[9]]
## A
## 690.5737
##
## [[10]]
## A
## -785.3561
##
## [[11]]
## A
## 63.14889
##
## [[12]]
## A
## 182.8178
##
## [[13]]
## A
## -281.2834
##
## [[14]]
## A
## -308.2917
##
## [[15]]
## A
## 1697.68
##
## [[16]]
## A
## 867.248
##
## [[17]]
## A
## -32.24963
##
## [[18]]
## A
## 942.6121
##
## [[19]]
## A
## -578.9598
##
## [[20]]
## A
## 446.4653
##
## [[21]]
## A
## -1158.553
##
## [[22]]
## A
## -632.3341
##
## [[23]]
## A
## 502.7908
##
## [[24]]
## A
## 141.4948
##
## [[25]]
## A
## -286.6895
##
## [[26]]
## A
## 409.3248
##
## [[27]]
## A
## 8.636311
##
## [[28]]
## A
## -721.3419
##
## [[29]]
## A
## -469.6626
##
## [[30]]
## A
## -63.14419
##
## [[31]]
## A
## -148.132
##
## [[32]]
## A
## 2054.757
##
## [[33]]
## A
## 775.0675
##
## [[34]]
## A
## -1075.856
##
## [[35]]
## A
## -220.0048
##
## [[36]]
## A
## 242.3467
##
## [[37]]
## A
## 926.9163
##
## [[38]]
## A
## 209.3311
##
## [[39]]
## A
## -633.595
##
## [[40]]
## A
## -490.8118
##
## [[41]]
## A
## -129.6891
##
## [[42]]
## A
## 263.329
##
## [[43]]
## A
## 155.1403
##
## [[44]]
## A
## -605.3427
##
## [[45]]
## A
## 97.80531
##
## [[46]]
## A
## 1122.096
##
## [[47]]
## A
## 784.3197
##
## [[48]]
## A
## -64.32241
##
## [[49]]
## A
## 1381.973
##
## [[50]]
## A
## 374.3534
##
## [[51]]
## A
## 258.5779
##
## [[52]]
## A
## 509.2898
##
## [[53]]
## A
## 680.1677
##
## [[54]]
## A
## 944.2358
##
## [[55]]
## A
## 112.1028
##
## [[56]]
## A
## 1165.506
##
## [[57]]
## A
## -661.4141
##
## [[58]]
## A
## 672.6627
##
## [[59]]
## A
## 359.6226
##
## [[60]]
## A
## -156.6569
##
## [[61]]
## A
## -553.7129
##
## [[62]]
## A
## 16.24596
##
## [[63]]
## A
## -106.0356
##
## [[64]]
## A
## 148.4261
##
## [[65]]
## A
## -46.31572
##
## [[66]]
## A
## -908.4822
##
## [[67]]
## A
## 324.5564
##
## [[68]]
## A
## -562.3346
##
## [[69]]
## A
## 1473.534
##
## [[70]]
## A
## -16.56955
##
## [[71]]
## A
## -281.8012
##
## [[72]]
## A
## 1415.02
##
## [[73]]
## A
## -1850.082
##
## [[74]]
## A
## 72.11073
##
## [[75]]
## A
## 390.2092
##
## [[76]]
## A
## 1614.478
##
## [[77]]
## A
## -268.7432
##
## [[78]]
## A
## -612.8481
##
## [[79]]
## A
## -123.4329
##
## [[80]]
## A
## 434.6072
##
## [[81]]
## A
## -641.3223
##
## [[82]]
## A
## 334.6059
##
## [[83]]
## A
## -102.0074
##
## [[84]]
## A
## 487.94
##
## [[85]]
## A
## -188.2482
##
## [[86]]
## A
## -405.0941
##
## [[87]]
## A
## -823.0153
##
## [[88]]
## A
## 10.50086
##
## [[89]]
## A
## 201.977
##
## [[90]]
## A
## 292.9196
##
## [[91]]
## A
## -77.75123
##
## [[92]]
## A
## 247.2636
##
## [[93]]
## A
## -1015.611
##
## [[94]]
## A
## -943.2803
##
## [[95]]
## A
## 28.45914
##
## [[96]]
## A
## -957.4054
##
## [[97]]
## A
## -25.3347
##
## [[98]]
## A
## 130.4186
##
## [[99]]
## A
## 121.7164
##
## [[100]]
## A
## -187.3322
##
## [[101]]
## A
## 293.4802
##
## [[102]]
## A
## -530.9342
##
## [[103]]
## A
## 287.3895
##
## [[104]]
## A
## 206.9525
##
## [[105]]
## A
## -329.4872
##
## [[106]]
## A
## 386.0107
##
## [[107]]
## A
## 637.5518
##
## [[108]]
## A
## 944.7741
##
## [[109]]
## A
## -380.1261
##
## [[110]]
## A
## 330.8401
##
## [[111]]
## A
## 172.5047
##
## [[112]]
## A
## -739.2983
##
## [[113]]
## A
## 264.5795
##
## [[114]]
## A
## -517.9332
##
## [[115]]
## A
## 736.8732
##
## [[116]]
## A
## -959.3757
##
## [[117]]
## A
## 116.9961
##
## [[118]]
## A
## -237.7446
##
## [[119]]
## A
## -6.531052
##
## [[120]]
## A
## -57.80981
##
## [[121]]
## A
## 159.7923
##
## [[122]]
## A
## -705.1892
##
## [[123]]
## A
## 300.3208
##
## [[124]]
## A
## 555.6338
##
## [[125]]
## A
## 419.6932
##
## [[126]]
## A
## -398.6607
##
## [[127]]
## A
## -579.5169
##
## [[128]]
## A
## 1987.144
##
## [[129]]
## A
## 565.9161
##
## [[130]]
## A
## 213.81
##
## [[131]]
## A
## -589.8423
##
## [[132]]
## A
## -1029.885
##
## [[133]]
## A
## -642.6676
##
## [[134]]
## A
## -612.6775
##
## [[135]]
## A
## 462.7963
##
## [[136]]
## A
## 224.2919
##
## [[137]]
## A
## -428.5049
##
## [[138]]
## A
## -287.9008
##
## [[139]]
## A
## 192.9521
##
## [[140]]
## A
## 146.507
##
## [[141]]
## A
## 119.2881
##
## [[142]]
## A
## -663.9189
##
## [[143]]
## A
## 87.06965
##
## [[144]]
## A
## -399.5818
##
## [[145]]
## A
## -99.58102
##
## [[146]]
## A
## 13.17562
##
## [[147]]
## A
## 144.1305
##
## [[148]]
## A
## -240.6245
##
## [[149]]
## A
## -843.9004
##
## [[150]]
## A
## 965.8031
##
## [[151]]
## A
## 13.31286
##
## [[152]]
## A
## 317.2717
##
## [[153]]
## A
## -155.1309
##
## [[154]]
## A
## 700.8803
##
## [[155]]
## A
## -892.4155
##
## [[156]]
## A
## -454.5984
##
## [[157]]
## A
## -140.4079
##
## [[158]]
## A
## -195.5335
##
## [[159]]
## A
## 202.2339
##
## [[160]]
## A
## 225.7401
##
## [[161]]
## A
## 39.25384
##
## [[162]]
## A
## -257.3205
##
## [[163]]
## A
## 557.8326
##
## [[164]]
## A
## 804.1909
##
## [[165]]
## A
## 252.1739
##
## [[166]]
## A
## 921.0727
##
## [[167]]
## A
## 1036.484
##
## [[168]]
## A
## 487.0155
##
## [[169]]
## A
## 674.3787
##
## [[170]]
## A
## -415.4164
##
## [[171]]
## A
## 84.36327
##
## [[172]]
## A
## -713.5484
##
## [[173]]
## A
## -330.2761
##
## [[174]]
## A
## -913.1129
##
## [[175]]
## A
## 1206.603
##
## [[176]]
## A
## -807.1853
##
## [[177]]
## A
## -58.32316
##
## [[178]]
## A
## 271.5442
##
## [[179]]
## A
## 361.367
##
## [[180]]
## A
## -738.2631
##
## [[181]]
## A
## 1247.94
##
## [[182]]
## A
## 85.59053
##
## [[183]]
## A
## 1018.023
##
## [[184]]
## A
## -96.12453
##
## [[185]]
## A
## 416.7309
##
## [[186]]
## A
## -711.209
##
## [[187]]
## A
## 29.17745
##
## [[188]]
## A
## -209.9639
##
## [[189]]
## A
## -1719.571
##
## [[190]]
## A
## 558.0898
##
## [[191]]
## A
## 19.94923
##
## [[192]]
## A
## -398.5289
##
## [[193]]
## A
## -834.6078
##
## [[194]]
## A
## 1313.749
##
## [[195]]
## A
## 271.4056
##
## [[196]]
## A
## -145.501
##
## [[197]]
## A
## 46.06664
##
## [[198]]
## A
## 1231.735
##
## [[199]]
## A
## 468.4947
##
## [[200]]
## A
## 528.1817
##
## [[201]]
## A
## -1229.222
##
## [[202]]
## A
## 56.62934
##
## [[203]]
## A
## 217.5864
##
## [[204]]
## A
## -162.1189
##
## [[205]]
## A
## 914.5572
##
## [[206]]
## A
## 181.757
##
## [[207]]
## A
## -24.50433
##
## [[208]]
## A
## 43.62679
##
## [[209]]
## A
## 720.8048
##
## [[210]]
## A
## -684.8425
##
## [[211]]
## A
## 518.7947
##
## [[212]]
## A
## -712.324
##
## [[213]]
## A
## 942.0825
##
## [[214]]
## A
## 158.0256
##
## [[215]]
## A
## -379.8432
##
## [[216]]
## A
## 842.2654
##
## [[217]]
## A
## 1287.246
##
## [[218]]
## A
## -248.7993
##
## [[219]]
## A
## -198.0328
##
## [[220]]
## A
## -685.9907
##
## [[221]]
## A
## -46.28095
##
## [[222]]
## A
## 943.0708
##
## [[223]]
## A
## -90.52343
##
## [[224]]
## A
## 134.2142
##
## [[225]]
## A
## -824.1603
##
## [[226]]
## A
## 73.05398
##
## [[227]]
## A
## 27.51404
##
## [[228]]
## A
## 912.8358
##
## [[229]]
## A
## -729.1544
##
## [[230]]
## A
## 1165.121
##
## [[231]]
## A
## -586.1886
##
## [[232]]
## A
## 228.2474
##
## [[233]]
## A
## 86.01514
##
## [[234]]
## A
## -324.1238
##
## [[235]]
## A
## 141.2199
##
## [[236]]
## A
## -590.2103
##
## [[237]]
## A
## -418.7175
##
## [[238]]
## A
## 117.8619
##
## [[239]]
## A
## 57.4512
##
## [[240]]
## A
## -453.403
##
## [[241]]
## A
## 632.974
##
## [[242]]
## A
## 784.3844
##
## [[243]]
## A
## 897.4997
##
## [[244]]
## A
## -339.5268
##
## [[245]]
## A
## -498.7917
##
## [[246]]
## A
## -95.13365
##
## [[247]]
## A
## 288.3894
##
## [[248]]
## A
## -296.5924
##
## [[249]]
## A
## 1551.619
##
## [[250]]
## A
## -980.7171
##
## [[251]]
## A
## 129.2284
##
## [[252]]
## A
## 249.9609
##
## [[253]]
## A
## -1180.206
##
## [[254]]
## A
## -697.4946
##
## [[255]]
## A
## 738.189
##
## [[256]]
## A
## -16.72748
##
## [[257]]
## A
## -12.70111
##
## [[258]]
## A
## -5.668872
##
## [[259]]
## A
## 440.1413
##
## [[260]]
## A
## -15.11824
##
## [[261]]
## A
## 264.5417
##
## [[262]]
## A
## 1414.8
##
## [[263]]
## A
## -187.4039
##
## [[264]]
## A
## 163.7146
##
## [[265]]
## A
## -202.8112
##
## [[266]]
## A
## 603.3834
##
## [[267]]
## A
## 8.265526
##
## [[268]]
## A
## 432.2541
##
## [[269]]
## A
## -84.72551
##
## [[270]]
## A
## -120.1497
##
## [[271]]
## A
## 1097.676
##
## [[272]]
## A
## -595.4002
##
## [[273]]
## A
## 489.9777
##
## [[274]]
## A
## -527.3288
##
## [[275]]
## A
## 831.3689
##
## [[276]]
## A
## 157.6372
##
## [[277]]
## A
## 455.1103
##
## [[278]]
## A
## 73.60286
##
## [[279]]
## A
## -140.5726
##
## [[280]]
## A
## 142.3358
##
## [[281]]
## A
## -145.6653
##
## [[282]]
## A
## 1106.019
##
## [[283]]
## A
## 32.56575
##
## [[284]]
## A
## -239.2024
##
## [[285]]
## A
## -1414.086
##
## [[286]]
## A
## -245.4396
##
## [[287]]
## A
## -999.5688
##
## [[288]]
## A
## 81.89557
##
## [[289]]
## A
## -199.7149
##
## [[290]]
## A
## 1275.157
##
## [[291]]
## A
## 426.3478
##
## [[292]]
## A
## -151.0131
##
## [[293]]
## A
## -328.9937
##
## [[294]]
## A
## 144.1821
##
## [[295]]
## A
## -345.3031
##
## [[296]]
## A
## 442.3834
##
## [[297]]
## A
## -1063.458
##
## [[298]]
## A
## -209.6456
##
## [[299]]
## A
## -235.6239
##
## [[300]]
## A
## -1210.055
##
## [[301]]
## A
## 553.7294
##
## [[302]]
## A
## -978.5225
##
## [[303]]
## A
## -634.8322
##
## [[304]]
## A
## 368.6036
##
## [[305]]
## A
## 650.7982
##
## [[306]]
## A
## 493.24
##
## [[307]]
## A
## -404.9121
##
## [[308]]
## A
## 819.8974
##
## [[309]]
## A
## 820.796
##
## [[310]]
## A
## 994.473
##
## [[311]]
## A
## 524.8135
##
## [[312]]
## A
## -21.9187
##
## [[313]]
## A
## 175.9918
##
## [[314]]
## A
## -1008.33
##
## [[315]]
## A
## -153.9023
##
## [[316]]
## A
## -952.5459
##
## [[317]]
## A
## -165.4528
##
## [[318]]
## A
## -417.8763
##
## [[319]]
## A
## -386.6656
##
## [[320]]
## A
## 41.34154
##
## [[321]]
## A
## -1285.265
##
## [[322]]
## A
## 269.7467
##
## [[323]]
## A
## -433.2725
##
## [[324]]
## A
## -371.9149
##
## [[325]]
## A
## 452.3918
##
## [[326]]
## A
## -511.1743
##
## [[327]]
## A
## 335.7056
##
## [[328]]
## A
## 141.305
##
## [[329]]
## A
## -825.2939
##
## [[330]]
## A
## -147.3167
##
## [[331]]
## A
## -145.6144
##
## [[332]]
## A
## -472.1275
##
## [[333]]
## A
## -300.1277
##
## [[334]]
## A
## 95.24496
##
## [[335]]
## A
## 714.7523
##
## [[336]]
## A
## -419.0102
##
## [[337]]
## A
## 792.1904
##
## [[338]]
## A
## 14.38992
##
## [[339]]
## A
## 481.1871
##
## [[340]]
## A
## 515.5395
##
## [[341]]
## A
## -652.2129
##
## [[342]]
## A
## 625.2532
##
## [[343]]
## A
## 778.7014
##
## [[344]]
## A
## 863.2858
##
## [[345]]
## A
## 861.0631
##
## [[346]]
## A
## -858.7755
##
## [[347]]
## A
## 390.5651
##
## [[348]]
## A
## 305.0473
##
## [[349]]
## A
## -449.9604
##
## [[350]]
## A
## 749.1389
##
## [[351]]
## A
## -170.1839
##
## [[352]]
## A
## -95.74228
##
## [[353]]
## A
## 71.52542
##
## [[354]]
## A
## -1020.475
##
## [[355]]
## A
## 643.4905
##
## [[356]]
## A
## 980.2366
##
## [[357]]
## A
## 787.035
##
## [[358]]
## A
## 731.663
##
## [[359]]
## A
## -353.3338
##
## [[360]]
## A
## -695.3924
##
## [[361]]
## A
## -81.63295
##
## [[362]]
## A
## 73.52
##
## [[363]]
## A
## -276.9304
##
## [[364]]
## A
## 443.5252
##
## [[365]]
## A
## 662.947
##
## [[366]]
## A
## -203.3023
##
## [[367]]
## A
## 149.7133
##
## [[368]]
## A
## 632.9308
##
## [[369]]
## A
## -114.8507
##
## [[370]]
## A
## 74.61184
##
## [[371]]
## A
## 853.7801
##
## [[372]]
## A
## 182.654
##
## [[373]]
## A
## -292.9115
##
## [[374]]
## A
## 698.9286
##
## [[375]]
## A
## 1178.12
##
## [[376]]
## A
## -266.5537
##
## [[377]]
## A
## 830.1311
##
## [[378]]
## A
## -750.274
##
## [[379]]
## A
## -651.7562
##
## [[380]]
## A
## -702.2432
##
## [[381]]
## A
## 582.8539
##
## [[382]]
## A
## -21.63603
##
## [[383]]
## A
## -469.4594
##
## [[384]]
## A
## -1037.074
##
## [[385]]
## A
## -58.97226
##
## [[386]]
## A
## 243.5284
##
## [[387]]
## A
## 164.1771
##
## [[388]]
## A
## 405.4249
##
## [[389]]
## A
## -1344.72
##
## [[390]]
## A
## -623.9568
##
## [[391]]
## A
## 350.0205
##
## [[392]]
## A
## 159.9289
##
## [[393]]
## A
## 155.1457
##
## [[394]]
## A
## -461.1445
##
## [[395]]
## A
## -85.62426
##
## [[396]]
## A
## 832.6441
##
## [[397]]
## A
## -848.382
##
## [[398]]
## A
## -134.2565
##
## [[399]]
## A
## -251.6364
##
## [[400]]
## A
## -664.0964
##
## [[401]]
## A
## -468.1463
##
## [[402]]
## A
## -49.22138
##
## [[403]]
## A
## -76.91984
##
## [[404]]
## A
## -320.3761
##
## [[405]]
## A
## 395.8724
##
## [[406]]
## A
## 109.334
##
## [[407]]
## A
## -337.5135
##
## [[408]]
## A
## 717.1054
##
## [[409]]
## A
## 702.5986
##
## [[410]]
## A
## -564.5642
##
## [[411]]
## A
## -188.9936
##
## [[412]]
## A
## -287.9507
##
## [[413]]
## A
## 1916.55
##
## [[414]]
## A
## 612.6224
##
## [[415]]
## A
## 220.7013
##
## [[416]]
## A
## 12.83211
##
## [[417]]
## A
## -729.4072
##
## [[418]]
## A
## 252.3679
##
## [[419]]
## A
## 838.2812
##
## [[420]]
## A
## 211.7105
##
## [[421]]
## A
## 717.9222
##
## [[422]]
## A
## -151.6067
##
## [[423]]
## A
## -175.5331
##
## [[424]]
## A
## 887.9908
##
## [[425]]
## A
## -136.5684
##
## [[426]]
## A
## 423.5403
##
## [[427]]
## A
## 1200.387
##
## [[428]]
## A
## -665.3687
##
## [[429]]
## A
## 328.4249
##
## [[430]]
## A
## -459.6049
##
## [[431]]
## A
## 176.7189
##
## [[432]]
## A
## -647.8855
##
## [[433]]
## A
## -317.3965
##
## [[434]]
## A
## 355.2094
##
## [[435]]
## A
## -1210.897
##
## [[436]]
## A
## -533.5167
##
## [[437]]
## A
## 57.23258
##
## [[438]]
## A
## 152.2181
##
## [[439]]
## A
## 331.2241
##
## [[440]]
## A
## 607.1667
##
## [[441]]
## A
## 102.6747
##
## [[442]]
## A
## 577.6833
##
## [[443]]
## A
## -720.8701
##
## [[444]]
## A
## 182.8281
##
## [[445]]
## A
## 52.55657
##
## [[446]]
## A
## 444.1893
##
## [[447]]
## A
## 79.18336
##
## [[448]]
## A
## -1254.47
##
## [[449]]
## A
## -269.49
##
## [[450]]
## A
## 516.5715
##
## [[451]]
## A
## -246.9689
##
## [[452]]
## A
## -23.1056
##
## [[453]]
## A
## 40.16826
##
## [[454]]
## A
## 2197.291
##
## [[455]]
## A
## 10.69798
##
## [[456]]
## A
## -1555.684
##
## [[457]]
## A
## -325.8166
##
## [[458]]
## A
## -609.9764
##
## [[459]]
## A
## 379.6217
##
## [[460]]
## A
## 877.076
##
## [[461]]
## A
## 253.9645
##
## [[462]]
## A
## -530.1104
##
## [[463]]
## A
## -860.3629
##
## [[464]]
## A
## 432.7825
##
## [[465]]
## A
## -409.7336
##
## [[466]]
## A
## 247.6428
##
## [[467]]
## A
## -337.4301
##
## [[468]]
## A
## -1238.336
##
## [[469]]
## A
## 159.4662
##
## [[470]]
## A
## -66.93999
##
## [[471]]
## A
## 670.3898
##
## [[472]]
## A
## -786.3889
##
## [[473]]
## A
## 457.482
##
## [[474]]
## A
## -212.3561
##
## [[475]]
## A
## 788.2473
##
## [[476]]
## A
## -388.0646
##
## [[477]]
## A
## -227.1111
##
## [[478]]
## A
## 276.6012
##
## [[479]]
## A
## -287.2339
##
## [[480]]
## A
## 867.4942
##
## [[481]]
## A
## 295.174
##
## [[482]]
## A
## -1481.751
##
## [[483]]
## A
## -63.01531
##
## [[484]]
## A
## 512.7197
##
## [[485]]
## A
## 1257.141
##
## [[486]]
## A
## 959.0657
##
## [[487]]
## A
## -1086.602
##
## [[488]]
## A
## 890.4503
##
## [[489]]
## A
## -42.45197
##
## [[490]]
## A
## 560.5719
##
## [[491]]
## A
## -1414.248
##
## [[492]]
## A
## 292.9131
##
## [[493]]
## A
## 1134.361
##
## [[494]]
## A
## 322.6629
##
## [[495]]
## A
## 212.5617
##
## [[496]]
## A
## 420.621
##
## [[497]]
## A
## 169.5224
##
## [[498]]
## A
## -1028.443
##
## [[499]]
## A
## -335.9968
##
## [[500]]
## A
## -197.7704hist(as.numeric(res1))
abs(as.numeric(res1))>abs(t1$statistic)## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [16] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [31] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [46] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [61] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [76] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [91] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [106] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [121] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [136] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [151] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [166] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [181] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [196] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [211] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [226] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [241] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [256] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [271] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [286] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [301] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [316] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [331] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [346] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [361] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [376] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [391] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [406] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [421] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [436] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [451] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [466] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [481] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [496] TRUE TRUE TRUE TRUE TRUEmean(abs(as.numeric(res1))>abs(t1$statistic))## [1] 1