Class8 : Functions
0505
In-class exercise
Q1.
Split the ChickWeight{datasets} data by individual chicks to extract separate slope estimates of regressing weight onto Time for each chick.
m0 <- with(ChickWeight, by(ChickWeight, Chick, function(x) lm(weight ~ Time, data = ChickWeight)))
sapply(m0, coef)## 18 16 15 13 9 20
## (Intercept) 27.467425 27.467425 27.467425 27.467425 27.467425 27.467425
## Time 8.803039 8.803039 8.803039 8.803039 8.803039 8.803039
## 10 8 17 19 4 6
## (Intercept) 27.467425 27.467425 27.467425 27.467425 27.467425 27.467425
## Time 8.803039 8.803039 8.803039 8.803039 8.803039 8.803039
## 11 3 1 12 2 5
## (Intercept) 27.467425 27.467425 27.467425 27.467425 27.467425 27.467425
## Time 8.803039 8.803039 8.803039 8.803039 8.803039 8.803039
## 14 7 24 30 22 23
## (Intercept) 27.467425 27.467425 27.467425 27.467425 27.467425 27.467425
## Time 8.803039 8.803039 8.803039 8.803039 8.803039 8.803039
## 27 28 26 25 29 21
## (Intercept) 27.467425 27.467425 27.467425 27.467425 27.467425 27.467425
## Time 8.803039 8.803039 8.803039 8.803039 8.803039 8.803039
## 33 37 36 31 39 38
## (Intercept) 27.467425 27.467425 27.467425 27.467425 27.467425 27.467425
## Time 8.803039 8.803039 8.803039 8.803039 8.803039 8.803039
## 32 40 34 35 44 45
## (Intercept) 27.467425 27.467425 27.467425 27.467425 27.467425 27.467425
## Time 8.803039 8.803039 8.803039 8.803039 8.803039 8.803039
## 43 41 47 49 46 50
## (Intercept) 27.467425 27.467425 27.467425 27.467425 27.467425 27.467425
## Time 8.803039 8.803039 8.803039 8.803039 8.803039 8.803039
## 42 48
## (Intercept) 27.467425 27.467425
## Time 8.803039 8.803039Q2.
Explain what does this statement do: lapply(lapply(search(), ls), length)
Q3.
The following R script uses Cushings{MASS} to demonstrates several ways to achieve the same objective in R. Explain the advantages or disadvantages of each method.
# method 2
# 先把不同Type的資料拆成list的資料，在用sapply一個個執行by column的平均
sapply(split(Cushings[,-3], Cushings$Type), function(x) apply(x, 2, mean))## a b c u
## Tetrahydrocortisone 2.966667 8.18 19.72 14.01667
## Pregnanetriol 2.440000 1.12 5.50 1.20000# method 3
# 我們把data這個data frame按照INDICES的factor拆分成若干塊小的data frames，在每塊小的data frame上運行函數FUN。
# function 先把Type去除assign給y，再把y執行by column的平均
# 接下來把每一個table轉換成list再rbind再一起
do.call("rbind", as.list(
by(Cushings, list(Cushings$Type), function(x) {
y <- subset(x, select = -Type)
apply(y, 2, mean)
}
)))## Tetrahydrocortisone Pregnanetriol
## a 2.966667 2.44
## b 8.180000 1.12
## c 19.720000 5.50
## u 14.016667 1.20# method 4
# using Tidy approach
# 資料先以Type group起來，再summarize
Cushings %>%
group_by(Type) %>%
summarize( t_m = mean(Tetrahydrocortisone), p_m = mean(Pregnanetriol))# method 5
# 使用nest＆map、map_dbl方法
Cushings %>%
nest(-Type) %>%
mutate(avg = map(data, ~ apply(., 2, mean)),
res_1 = map_dbl(avg, "Tetrahydrocortisone"),
res_2 = map_dbl(avg, "Pregnanetriol")) ## Warning: All elements of `...` must be named.
## Did you want `data = c(Tetrahydrocortisone, Pregnanetriol)`?Q4.
Go through the script in the NZ schools example and provide comments to each code chunk indicated by ‘##’. Give alternative code to perform the same calculation where appropriate.
## keep the school names with white spaces
dta <- read.csv("nzSchools.csv", as.is=2)
## 查看資料結構
str(dta)## 'data.frame': 2571 obs. of 6 variables:
## $ ID : int 1015 1052 1062 1092 1130 1018 1029 1030 1588 1154 ...
## $ Name: chr "Hora Hora School" "Morningside School" "Onerahi School" "Raurimu Avenue School" ...
## $ City: Factor w/ 541 levels "Ahaura","Ahipara",..: 533 533 533 533 533 533 533 533 533 533 ...
## $ Auth: Factor w/ 4 levels "Other","Private",..: 3 3 3 3 3 3 3 3 4 3 ...
## $ Dec : int 2 3 4 2 4 8 5 5 6 1 ...
## $ Roll: int 318 200 455 86 577 329 637 395 438 201 ...## [1] 2571 6## binning
## 每一行Roll若大於中位數，則Size為Large否則Small
dta$Size <- ifelse(dta$Roll > median(dta$Roll), "Large", "Small")
# 改寫
library(magrittr)##
## Attaching package: 'magrittr'## The following object is masked from 'package:purrr':
##
## set_names## The following object is masked from 'package:tidyr':
##
## extractdta %<>% mutate(Size = cut(Roll, c(min(dta$Roll),median(dta$Roll),max(dta$Roll)), c("Small", "Large")))
## 刪掉dta$Size
dta$Size <- NULL
## 看前六筆
head(dta)## 不知道是怎麼切的
dta$Size <- cut(dta$Roll, 3, labels=c("Small", "Mediam", "Large"))
## 看每個分類的個數
table(dta$Size)##
## Small Mediam Large
## 2555 15 1# 改寫
do.call(rbind, with(dta, tapply(Roll, list(Size), function(x) c(min = min(x), max = max(x), len = length(x))))) ## min max len
## Small 5 1821 2555
## Mediam 1888 3022 15
## Large 5546 5546 1## sorting
## Roll以大到小排列
dta$RollOrd <- order(dta$Roll, decreasing=T)
## 以Roll的排序看前六筆資料
head(dta[dta$RollOrd, ])##
## Other Private State State Integrated
## 1 99 2144 327##
## Other Private State State Integrated
## 1 99 2144 327## [1] "table"## Dec
## Auth 1 2 3 4 5 6 7 8 9 10
## Other 1 0 0 0 0 0 0 0 0 0
## Private 0 0 2 6 2 2 6 11 12 38
## State 259 230 208 219 214 215 188 200 205 205
## State Integrated 12 22 35 28 38 34 45 45 37 31## [1] 295.4737## [1] 308.798## dta$Auth: Other
## [1] 51
## --------------------------------------------------------
## dta$Auth: Private
## [1] 308.798
## --------------------------------------------------------
## dta$Auth: State
## [1] 300.6301
## --------------------------------------------------------
## dta$Auth: State Integrated
## [1] 258.3792## : Other
## [1] 51
## --------------------------------------------------------
## : Private
## [1] 308.798
## --------------------------------------------------------
## : State
## [1] 300.6301
## --------------------------------------------------------
## : State Integrated
## [1] 258.3792## Other Private State State Integrated
## 51.0000 308.7980 300.6301 258.3792## splitting
## 以Auth分組把Roll分開
rollsByAuth <- split(dta$Roll, dta$Auth)
str(rollsByAuth);class(rollsByAuth)## List of 4
## $ Other : int 51
## $ Private : int [1:99] 255 39 154 73 83 25 95 85 94 729 ...
## $ State : int [1:2144] 318 200 455 86 577 329 637 395 201 267 ...
## $ State Integrated: int [1:327] 438 26 191 560 151 114 126 171 211 57 ...## [1] "list"## $Other
## [1] 51
##
## $Private
## [1] 308.798
##
## $State
## [1] 300.6301
##
## $`State Integrated`
## [1] 258.3792## Other Private State State Integrated
## 51.0000 308.7980 300.6301 258.3792## [1] FALSE FALSE FALSE FALSE FALSE TRUE## : Other
## : FALSE
## [1] 51
## --------------------------------------------------------
## : Private
## : FALSE
## [1] 151.4
## --------------------------------------------------------
## : State
## : FALSE
## [1] 261.7487
## --------------------------------------------------------
## : State Integrated
## : FALSE
## [1] 183.237
## --------------------------------------------------------
## : Other
## : TRUE
## [1] NA
## --------------------------------------------------------
## : Private
## : TRUE
## [1] 402.5362
## --------------------------------------------------------
## : State
## : TRUE
## [1] 338.8243
## --------------------------------------------------------
## : State Integrated
## : TRUE
## [1] 311.2135## : Other
## : FALSE
## [1] 51
## --------------------------------------------------------
## : Private
## : FALSE
## [1] 151.4
## --------------------------------------------------------
## : State
## : FALSE
## [1] 261.7487
## --------------------------------------------------------
## : State Integrated
## : FALSE
## [1] 183.237
## --------------------------------------------------------
## : Other
## : TRUE
## [1] NA
## --------------------------------------------------------
## : Private
## : TRUE
## [1] 402.5362
## --------------------------------------------------------
## : State
## : TRUE
## [1] 338.8243
## --------------------------------------------------------
## : State Integrated
## : TRUE
## [1] 311.2135## FALSE TRUE
## Other 51.0000 NA
## Private 151.4000 402.5362
## State 261.7487 338.8243
## State Integrated 183.2370 311.2135## $Other.FALSE
## [1] 51
##
## $Private.FALSE
## [1] 151.4
##
## $State.FALSE
## [1] 261.7487
##
## $`State Integrated.FALSE`
## [1] 183.237
##
## $Other.TRUE
## [1] NaN
##
## $Private.TRUE
## [1] 402.5362
##
## $State.TRUE
## [1] 338.8243
##
## $`State Integrated.TRUE`
## [1] 311.2135## Other.FALSE Private.FALSE State.FALSE
## 51.0000 151.4000 261.7487
## State Integrated.FALSE Other.TRUE Private.TRUE
## 183.2370 NaN 402.5362
## State.TRUE State Integrated.TRUE
## 338.8243 311.2135## : Other
## [1] 51 51
## --------------------------------------------------------
## : Private
## [1] 7 1663
## --------------------------------------------------------
## : State
## [1] 5 5546
## --------------------------------------------------------
## : State Integrated
## [1] 18 1475## $Other
## [1] 51 51
##
## $Private
## [1] 7 1663
##
## $State
## [1] 5 5546
##
## $`State Integrated`
## [1] 18 1475## $Other
## [1] 51 51
##
## $Private
## [1] 7 1663
##
## $State
## [1] 5 5546
##
## $`State Integrated`
## [1] 18 1475## Other Private State State Integrated
## [1,] 51 7 5 18
## [2,] 51 1663 5546 1475Q5.
Go through the script in the NCEA 2007 example and provide comments to each code chunk indicated by ‘##’. Give alternative code to perform the same calculation where appropriate.
## [1] 88 4## 'data.frame': 88 obs. of 4 variables:
## $ Name : chr "Al-Madinah School" "Alfriston College" "Ambury Park Centre for Riding Therapy" "Aorere College" ...
## $ Level1: num 61.5 53.9 33.3 39.5 71.2 22.1 50.8 57.3 89.3 59.8 ...
## $ Level2: num 75 44.1 20 50.2 78.9 30.8 34.8 49.8 89.7 65.7 ...
## $ Level3: num 0 0 0 30.6 55.5 26.3 48.9 44.6 88.6 50.4 ...## Level1 Level2 Level3
## 62.26705 61.06818 47.97614## $Level1
## [1] 62.26705
##
## $Level2
## [1] 61.06818
##
## $Level3
## [1] 47.97614## Level1 Level2 Level3
## 62.26705 61.06818 47.97614## Level1 Level2 Level3
## 62.26705 61.06818 47.97614## Level1 Level2 Level3
## [1,] 2.8 0.0 0.0
## [2,] 97.4 95.7 95.7## $Level1
## [1] 2.8 97.4
##
## $Level2
## [1] 0.0 95.7
##
## $Level3
## [1] 0.0 95.7## Level1 Level2 Level3
## [1,] 2.8 0.0 0.0
## [2,] 97.4 95.7 95.7Exercises
Q1.
Use the data in the high schools example to solve the following problems:
- test if any pairs of the five variables: read, write, math, science, and socst, are different in means.
summ <- data.frame(mean=apply(dta[,7:11],2,mean,na.rm=T),sd=apply(dta[,7:11],2,sd,na.rm=T),len=rep(dim(dta)[1],5))
summ$se <- summ$sd/sqrt(summ$len)
summ$subject <- row.names(summ)
summsumm %>% ggplot() +
aes(reorder(subject, mean), mean) +
geom_errorbar(aes(ymin=mean-se,
ymax=mean+se),
width=.2, size=.3) +
geom_point(size=rel(3))
t <- dta %>% dplyr::select(read, write, math, science, socst) %>% gather(., key = key, value = score)
# anova
anova(lm(score ~ key, data=t))##
## Pairwise comparisons using t tests with pooled SD
##
## data: t$score and t$key
##
## math read science socst
## read 0.68 - - -
## science 0.47 0.76 - -
## socst 0.81 0.86 0.63 -
## write 0.90 0.58 0.39 0.71
##
## P value adjustment method: none## read write math science socst
## 0.8676815 0.3775863 1.0000000 0.3775863 0.8676815## read write math science socst
## [1,] "read" "write" "math" "science" "socst"
## [2,] "socst" "science" "math" "write" "read"皆無顯著不同
- test if the 4 different ethnic groups have the same mean scores for each of the 5 variables (individually): read, write, math, science, and socst.
dta %>%
dplyr::select(race, read, write, math, science, socst) %>%
gather(-race, key = subject, value = score) %>%
group_by(race, subject) %>%
summarize(mean = mean(score, na.rm = T),
se = sd(score, na.rm = T)/sqrt(n())) %>%
ggplot() +
facet_wrap(. ~ race, nrow = 1) +
aes(reorder(subject, mean), mean) +
geom_errorbar(aes(ymin=mean-se,
ymax=mean+se),
width=.2, size=.3) +
geom_point(size=rel(2),
aes(color = subject)) +
theme(legend.position = 'none',
axis.text.x = element_text(angle = 45))
t <- dta %>% dplyr::select(race, read, write, math, science, socst) %>% gather(-race, key = key, value = score)
# method 1pairwise.t.test
by(t, t$race, function(x) pairwise.t.test(x$score, x$key, p.adjust.method="none"))## t$race: african-amer
##
## Pairwise comparisons using t tests with pooled SD
##
## data: x$score and x$key
##
## math read science socst
## read 0.986 - - -
## science 0.128 0.124 - -
## socst 0.335 0.344 0.014 -
## write 0.604 0.616 0.043 0.655
##
## P value adjustment method: none
## --------------------------------------------------------
## t$race: asian
##
## Pairwise comparisons using t tests with pooled SD
##
## data: x$score and x$key
##
## math read science socst
## read 0.170 - - -
## science 0.137 0.907 - -
## socst 0.110 0.815 0.907 -
## write 0.851 0.120 0.096 0.075
##
## P value adjustment method: none
## --------------------------------------------------------
## t$race: hispanic
##
## Pairwise comparisons using t tests with pooled SD
##
## data: x$score and x$key
##
## math read science socst
## read 0.76 - - -
## science 0.63 0.86 - -
## socst 0.88 0.65 0.53 -
## write 0.70 0.93 0.92 0.59
##
## P value adjustment method: none
## --------------------------------------------------------
## t$race: white
##
## Pairwise comparisons using t tests with pooled SD
##
## data: x$score and x$key
##
## math read science socst
## read 0.97 - - -
## science 0.88 0.85 - -
## socst 0.80 0.83 0.69 -
## write 0.94 0.91 0.94 0.75
##
## P value adjustment method: nonelapply(split(dta[,7:11], list(dta$race)), function(x) {
t <- x %>% gather(., key = key, value = score)
pairwise.t.test(t$score, t$key, p.adjust.method="none")$p.value
})## $`african-amer`
## math read science socst
## read 0.9857150 NA NA NA
## science 0.1283322 0.1240346 NA NA
## socst 0.3348230 0.3438066 0.01448692 NA
## write 0.6038592 0.6163771 0.04333649 0.6546046
##
## $asian
## math read science socst
## read 0.1702313 NA NA NA
## science 0.1374818 0.9065999 NA NA
## socst 0.1099473 0.8145132 0.90659992 NA
## write 0.8511080 0.1203698 0.09570242 0.07536262
##
## $hispanic
## math read science socst
## read 0.7604310 NA NA NA
## science 0.6296594 0.8565836 NA NA
## socst 0.8788077 0.6474875 0.5268951 NA
## write 0.6968597 0.9324883 0.9227852 0.5879346
##
## $white
## math read science socst
## read 0.9664872 NA NA NA
## science 0.8792554 0.8464059 NA NA
## socst 0.8009781 0.8336136 0.6872001 NA
## write 0.9425822 0.9092070 0.9360251 0.7458623# t.test
lapply(split(dta[,7:11], list(dta$race)), function(x) {
mapply(function(a,b) t.test(a,b)$p.value, x[,1:5], x[,5:1])
})## $`african-amer`
## read write math science socst
## 0.36781536 0.06378806 1.00000000 0.06378806 0.36781536
##
## $asian
## read write math science socst
## 0.81045778 0.09453115 1.00000000 0.09453115 0.81045778
##
## $hispanic
## read write math science socst
## 0.6914577 0.9159023 1.0000000 0.9159023 0.6914577
##
## $white
## read write math science socst
## 0.8456548 0.9316598 1.0000000 0.9316598 0.8456548## Df Sum Sq Mean Sq F value Pr(>F)
## race 3 8595 2865.0 31.672 <2e-16 ***
## key 4 98 24.4 0.270 0.898
## race:key 12 1018 84.8 0.938 0.508
## Residuals 975 88197 90.5
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 5 observations deleted due to missingness## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = score ~ race * key, data = t)
##
## $race
## diff lwr upr p adj
## asian-african-amer 7.15959596 3.043308 11.275884 0.0000502
## hispanic-african-amer 0.14828962 -3.181230 3.477809 0.9994614
## white-african-amer 7.18800190 4.564812 9.811192 0.0000000
## hispanic-asian -7.01130634 -11.002145 -3.020468 0.0000408
## white-asian 0.02840594 -3.395367 3.452179 0.9999965
## white-hispanic 7.03971228 4.618124 9.461300 0.0000000
##
## $key
## diff lwr upr p adj
## read-math -0.4150000 -3.014248 2.184248 0.9924683
## science-math -0.7598824 -3.375739 1.855974 0.9324144
## socst-math -0.2400000 -2.839248 2.359248 0.9991065
## write-math 0.1300000 -2.469248 2.729248 0.9999214
## science-read -0.3448824 -2.960739 2.270974 0.9963984
## socst-read 0.1750000 -2.424248 2.774248 0.9997439
## write-read 0.5450000 -2.054248 3.144248 0.9789740
## socst-science 0.5198824 -2.095974 3.135739 0.9827672
## write-science 0.8898824 -1.725974 3.505739 0.8853917
## write-socst 0.3700000 -2.229248 2.969248 0.9951543
##
## $`race:key`
## diff lwr
## asian:math-african-amer:math 10.52272727 -2.16478804
## hispanic:math-african-amer:math 0.66666667 -9.56656468
## white:math-african-amer:math 7.22241379 -0.83972084
## african-amer:read-african-amer:math 0.05000000 -10.63826156
## asian:read-african-amer:math 5.15909091 -7.52842441
## hispanic:read-african-amer:math -0.08333333 -10.31656468
## white:read-african-amer:math 7.17413793 -0.88799670
## african-amer:science-african-amer:math -4.32894737 -15.15693069
## asian:science-african-amer:math 4.70454545 -7.98296986
## hispanic:science-african-amer:math -0.53260870 -10.86646421
## white:science-african-amer:math 7.39788732 -0.67456354
## african-amer:socst-african-amer:math 2.70000000 -7.98826156
## asian:socst-african-amer:math 4.25000000 -8.43751532
## hispanic:socst-african-amer:math 1.04166667 -9.19156468
## white:socst-african-amer:math 6.93275862 -1.12937601
## african-amer:write-african-amer:math 1.45000000 -9.23826156
## asian:write-african-amer:math 11.25000000 -1.43751532
## hispanic:write-african-amer:math -0.29166667 -10.52489802
## white:write-african-amer:math 7.30517241 -0.75696222
## hispanic:math-asian:math -9.85606061 -22.16269027
## white:math-asian:math -3.30031348 -13.87065522
## african-amer:read-asian:math -10.47272727 -23.16024259
## asian:read-asian:math -5.36363636 -19.77568531
## hispanic:read-asian:math -10.60606061 -22.91269027
## white:read-asian:math -3.34858934 -13.91893109
## african-amer:science-asian:math -14.85167464 -27.65711617
## asian:science-asian:math -5.81818182 -20.23023076
## hispanic:science-asian:math -11.05533597 -23.44576290
## white:science-asian:math -3.12483995 -13.70305212
## african-amer:socst-asian:math -7.82272727 -20.51024259
## asian:socst-asian:math -6.27272727 -20.68477622
## hispanic:socst-asian:math -9.48106061 -21.78769027
## white:socst-asian:math -3.58996865 -14.16031040
## african-amer:write-asian:math -9.07272727 -21.76024259
## asian:write-asian:math 0.72727273 -13.68477622
## hispanic:write-asian:math -10.81439394 -23.12102360
## white:write-asian:math -3.21755486 -13.78789660
## white:math-hispanic:math 6.55574713 -0.89261535
## african-amer:read-hispanic:math -0.61666667 -10.84989802
## asian:read-hispanic:math 4.49242424 -7.81420542
## hispanic:read-hispanic:math -0.75000000 -10.50700326
## white:read-hispanic:math 6.50747126 -0.94089121
## african-amer:science-hispanic:math -4.99561404 -15.37469451
## asian:science-hispanic:math 4.03787879 -8.26875087
## hispanic:science-hispanic:math -1.19927536 -11.06176280
## white:science-hispanic:math 6.73122066 -0.72830692
## african-amer:socst-hispanic:math 2.03333333 -8.19989802
## asian:socst-hispanic:math 3.58333333 -8.72329633
## hispanic:socst-hispanic:math 0.37500000 -9.38200326
## white:socst-hispanic:math 6.26609195 -1.18227052
## african-amer:write-hispanic:math 0.78333333 -9.44989802
## asian:write-hispanic:math 10.58333333 -1.72329633
## hispanic:write-hispanic:math -0.95833333 -10.71533659
## white:write-hispanic:math 6.63850575 -0.80985673
## african-amer:read-white:math -7.17241379 -15.23454842
## asian:read-white:math -2.06332288 -12.63366463
## hispanic:read-white:math -7.30574713 -14.75410960
## white:read-white:math -0.04827586 -4.01779655
## african-amer:science-white:math -11.55136116 -19.79783329
## asian:science-white:math -2.51786834 -13.08821008
## hispanic:science-white:math -7.75502249 -15.34103919
## white:science-white:math 0.17547353 -3.81495786
## african-amer:socst-white:math -4.52241379 -12.58454842
## asian:socst-white:math -2.97241379 -13.54275554
## hispanic:socst-white:math -6.18074713 -13.62910960
## white:socst-white:math -0.28965517 -4.25917586
## african-amer:write-white:math -5.77241379 -13.83454842
## asian:write-white:math 4.02758621 -6.54275554
## hispanic:write-white:math -7.51408046 -14.96244293
## white:write-white:math 0.08275862 -3.88676207
## asian:read-african-amer:read 5.10909091 -7.57842441
## hispanic:read-african-amer:read -0.13333333 -10.36656468
## white:read-african-amer:read 7.12413793 -0.93799670
## african-amer:science-african-amer:read -4.37894737 -15.20693069
## asian:science-african-amer:read 4.65454545 -8.03296986
## hispanic:science-african-amer:read -0.58260870 -10.91646421
## white:science-african-amer:read 7.34788732 -0.72456354
## african-amer:socst-african-amer:read 2.65000000 -8.03826156
## asian:socst-african-amer:read 4.20000000 -8.48751532
## hispanic:socst-african-amer:read 0.99166667 -9.24156468
## white:socst-african-amer:read 6.88275862 -1.17937601
## african-amer:write-african-amer:read 1.40000000 -9.28826156
## asian:write-african-amer:read 11.20000000 -1.48751532
## hispanic:write-african-amer:read -0.34166667 -10.57489802
## white:write-african-amer:read 7.25517241 -0.80696222
## hispanic:read-asian:read -5.24242424 -17.54905391
## white:read-asian:read 2.01504702 -8.55529472
## african-amer:science-asian:read -9.48803828 -22.29347980
## asian:science-asian:read -0.45454545 -14.86659440
## hispanic:science-asian:read -5.69169960 -18.08212654
## white:science-asian:read 2.23879641 -8.33941576
## african-amer:socst-asian:read -2.45909091 -15.14660623
## asian:socst-asian:read -0.90909091 -15.32113985
## hispanic:socst-asian:read -4.11742424 -16.42405391
## white:socst-asian:read 1.77366771 -8.79667403
## african-amer:write-asian:read -3.70909091 -16.39660623
## asian:write-asian:read 6.09090909 -8.32113985
## hispanic:write-asian:read -5.45075758 -17.75738724
## white:write-asian:read 2.14608150 -8.42426024
## white:read-hispanic:read 7.25747126 -0.19089121
## african-amer:science-hispanic:read -4.24561404 -14.62469451
## asian:science-hispanic:read 4.78787879 -7.51875087
## hispanic:science-hispanic:read -0.44927536 -10.31176280
## white:science-hispanic:read 7.48122066 0.02169308
## african-amer:socst-hispanic:read 2.78333333 -7.44989802
## asian:socst-hispanic:read 4.33333333 -7.97329633
## hispanic:socst-hispanic:read 1.12500000 -8.63200326
## white:socst-hispanic:read 7.01609195 -0.43227052
## african-amer:write-hispanic:read 1.53333333 -8.69989802
## asian:write-hispanic:read 11.33333333 -0.97329633
## hispanic:write-hispanic:read -0.20833333 -9.96533659
## white:write-hispanic:read 7.38850575 -0.05985673
## african-amer:science-white:read -11.50308530 -19.74955743
## asian:science-white:read -2.46959248 -13.03993422
## hispanic:science-white:read -7.70674663 -15.29276333
## white:science-white:read 0.22374939 -3.76668200
## african-amer:socst-white:read -4.47413793 -12.53627256
## asian:socst-white:read -2.92413793 -13.49447968
## hispanic:socst-white:read -6.13247126 -13.58083374
## white:socst-white:read -0.24137931 -4.21090000
## african-amer:write-white:read -5.72413793 -13.78627256
## asian:write-white:read 4.07586207 -6.49447968
## hispanic:write-white:read -7.46580460 -14.91416707
## white:write-white:read 0.13103448 -3.83848621
## asian:science-african-amer:science 9.03349282 -3.77194870
## hispanic:science-african-amer:science 3.79633867 -6.68196546
## white:science-african-amer:science 11.72683469 3.47027665
## african-amer:socst-african-amer:science 7.02894737 -3.79903596
## asian:socst-african-amer:science 8.57894737 -4.22649416
## hispanic:socst-african-amer:science 5.37061404 -5.00846644
## white:socst-african-amer:science 11.26170599 3.01523386
## african-amer:write-african-amer:science 5.77894737 -5.04903596
## asian:write-african-amer:science 15.57894737 2.77350584
## hispanic:write-african-amer:science 4.03728070 -6.34179978
## white:write-african-amer:science 11.63411978 3.38764765
## hispanic:science-asian:science -5.23715415 -17.62758108
## white:science-asian:science 2.69334187 -7.88487030
## african-amer:socst-asian:science -2.00454545 -14.69206077
## asian:socst-asian:science -0.45454545 -14.86659440
## hispanic:socst-asian:science -3.66287879 -15.96950845
## white:socst-asian:science 2.22821317 -8.34212858
## african-amer:write-asian:science -3.25454545 -15.94206077
## asian:write-asian:science 6.54545455 -7.86659440
## hispanic:write-asian:science -4.99621212 -17.30284178
## white:write-asian:science 2.60062696 -7.96971479
## white:science-hispanic:science 7.93049602 0.33351651
## african-amer:socst-hispanic:science 3.23260870 -7.10124682
## asian:socst-hispanic:science 4.78260870 -7.60781824
## hispanic:socst-hispanic:science 1.57427536 -8.28821208
## white:socst-hispanic:science 7.46536732 -0.12064939
## african-amer:write-hispanic:science 1.98260870 -8.35124682
## asian:write-hispanic:science 11.78260870 -0.60781824
## hispanic:write-hispanic:science 0.24094203 -9.62154541
## white:write-hispanic:science 7.83778111 0.25176440
## african-amer:socst-white:science -4.69788732 -12.77033819
## asian:socst-white:science -3.14788732 -13.72609950
## hispanic:socst-white:science -6.35622066 -13.81574824
## white:socst-white:science -0.46512870 -4.45556009
## african-amer:write-white:science -5.94788732 -14.02033819
## asian:write-white:science 3.85211268 -6.72609950
## hispanic:write-white:science -7.68955399 -15.14908157
## white:write-white:science -0.09271491 -4.08314630
## asian:socst-african-amer:socst 1.55000000 -11.13751532
## hispanic:socst-african-amer:socst -1.65833333 -11.89156468
## white:socst-african-amer:socst 4.23275862 -3.82937601
## african-amer:write-african-amer:socst -1.25000000 -11.93826156
## asian:write-african-amer:socst 8.55000000 -4.13751532
## hispanic:write-african-amer:socst -2.99166667 -13.22489802
## white:write-african-amer:socst 4.60517241 -3.45696222
## hispanic:socst-asian:socst -3.20833333 -15.51496300
## white:socst-asian:socst 2.68275862 -7.88758312
## african-amer:write-asian:socst -2.80000000 -15.48751532
## asian:write-asian:socst 7.00000000 -7.41204894
## hispanic:write-asian:socst -4.54166667 -16.84829633
## white:write-asian:socst 3.05517241 -7.51516933
## white:socst-hispanic:socst 5.89109195 -1.55727052
## african-amer:write-hispanic:socst 0.40833333 -9.82489802
## asian:write-hispanic:socst 10.20833333 -2.09829633
## hispanic:write-hispanic:socst -1.33333333 -11.09033659
## white:write-hispanic:socst 6.26350575 -1.18485673
## african-amer:write-white:socst -5.48275862 -13.54489325
## asian:write-white:socst 4.31724138 -6.25310037
## hispanic:write-white:socst -7.22442529 -14.67278776
## white:write-white:socst 0.37241379 -3.59710690
## asian:write-african-amer:write 9.80000000 -2.88751532
## hispanic:write-african-amer:write -1.74166667 -11.97489802
## white:write-african-amer:write 5.85517241 -2.20696222
## hispanic:write-asian:write -11.54166667 -23.84829633
## white:write-asian:write -3.94482759 -14.51516933
## white:write-hispanic:write 7.59683908 0.14847661
## upr p adj
## asian:math-african-amer:math 23.21024259 0.2631891
## hispanic:math-african-amer:math 10.89989802 1.0000000
## white:math-african-amer:math 15.28454842 0.1483404
## african-amer:read-african-amer:math 10.73826156 1.0000000
## asian:read-african-amer:math 17.84660623 0.9968257
## hispanic:read-african-amer:math 10.14989802 1.0000000
## white:read-african-amer:math 15.23627256 0.1568362
## african-amer:science-african-amer:math 6.49903596 0.9974395
## asian:science-african-amer:math 17.39206077 0.9990490
## hispanic:science-african-amer:math 9.80124682 1.0000000
## white:science-african-amer:math 15.47033819 0.1217669
## african-amer:socst-african-amer:math 13.38826156 0.9999971
## asian:socst-african-amer:math 16.93751532 0.9997706
## hispanic:socst-african-amer:math 11.27489802 1.0000000
## white:socst-african-amer:math 14.99489325 0.2047564
## african-amer:write-african-amer:math 12.13826156 1.0000000
## asian:write-african-amer:math 23.93751532 0.1614594
## hispanic:write-african-amer:math 9.94156468 1.0000000
## white:write-african-amer:math 15.36730704 0.1345907
## hispanic:math-asian:math 2.45056906 0.3261825
## white:math-asian:math 7.27002827 0.9999185
## african-amer:read-asian:math 2.21478804 0.2714264
## asian:read-asian:math 9.04841258 0.9990007
## hispanic:read-asian:math 1.70056906 0.2014395
## white:read-asian:math 7.22175240 0.9998988
## african-amer:science-asian:math -2.04623312 0.0063886
## asian:science-asian:math 8.59386712 0.9971022
## hispanic:science-asian:math 1.33509096 0.1533994
## white:science-asian:math 7.45337222 0.9999649
## african-amer:socst-asian:math 4.86478804 0.8054915
## asian:socst-asian:math 8.13932167 0.9927227
## hispanic:socst-asian:math 2.82556906 0.4014787
## white:socst-asian:math 6.98037309 0.9997203
## african-amer:write-asian:math 3.61478804 0.5518237
## asian:write-asian:math 15.13932167 1.0000000
## hispanic:write-asian:math 1.49223572 0.1735494
## white:write-asian:math 7.35278689 0.9999444
## white:math-hispanic:math 14.00410960 0.1713528
## african-amer:read-hispanic:math 9.61656468 1.0000000
## asian:read-hispanic:math 16.79905391 0.9992315
## hispanic:read-hispanic:math 9.00700326 1.0000000
## white:read-hispanic:math 13.95583374 0.1815882
## african-amer:science-hispanic:math 5.38346644 0.9775931
## asian:science-hispanic:math 16.34450845 0.9998302
## hispanic:science-hispanic:math 8.66321208 1.0000000
## white:science-hispanic:math 14.19074824 0.1394872
## african-amer:socst-hispanic:math 12.26656468 0.9999999
## asian:socst-hispanic:math 15.88996300 0.9999719
## hispanic:socst-hispanic:math 10.13200326 1.0000000
## white:socst-hispanic:math 13.71445443 0.2392844
## african-amer:write-hispanic:math 11.01656468 1.0000000
## asian:write-hispanic:math 22.88996300 0.2046615
## hispanic:write-hispanic:math 8.79866992 1.0000000
## white:write-hispanic:math 14.08686822 0.1547971
## african-amer:read-white:math 0.88972084 0.1571462
## asian:read-white:math 8.50701886 1.0000000
## hispanic:read-white:math 0.14261535 0.0620609
## white:read-white:math 3.92124483 1.0000000
## african-amer:science-white:math -3.30488903 0.0001368
## asian:science-white:math 8.05247341 0.9999989
## hispanic:science-white:math -0.16900578 0.0385485
## white:science-white:math 4.16590492 1.0000000
## african-amer:socst-white:math 3.53972084 0.9049180
## asian:socst-white:math 7.59792795 0.9999837
## hispanic:socst-white:math 1.26761535 0.2622867
## white:socst-white:math 3.67986552 1.0000000
## african-amer:write-white:math 2.28972084 0.5493137
## asian:write-white:math 14.59792795 0.9986307
## hispanic:write-white:math -0.06571799 0.0451556
## white:write-white:math 4.05227931 1.0000000
## asian:read-african-amer:read 17.79660623 0.9971933
## hispanic:read-african-amer:read 10.09989802 1.0000000
## white:read-african-amer:read 15.18627256 0.1660118
## african-amer:science-african-amer:read 6.44903596 0.9970373
## asian:science-african-amer:read 17.34206077 0.9991775
## hispanic:science-african-amer:read 9.75124682 1.0000000
## white:science-african-amer:read 15.42033819 0.1293385
## african-amer:socst-african-amer:read 13.33826156 0.9999978
## asian:socst-african-amer:read 16.88751532 0.9998068
## hispanic:socst-african-amer:read 11.22489802 1.0000000
## white:socst-african-amer:read 14.94489325 0.2158400
## african-amer:write-african-amer:read 12.08826156 1.0000000
## asian:write-african-amer:read 23.88751532 0.1673702
## hispanic:write-african-amer:read 9.89156468 1.0000000
## white:write-african-amer:read 15.31730704 0.1427760
## hispanic:read-asian:read 7.06420542 0.9943695
## white:read-asian:read 12.58538877 1.0000000
## african-amer:science-asian:read 3.31740325 0.4801739
## asian:science-asian:read 13.95750349 1.0000000
## hispanic:science-asian:read 6.69872733 0.9864870
## white:science-asian:read 12.81700859 0.9999998
## african-amer:socst-asian:read 10.22842441 1.0000000
## asian:socst-asian:read 13.50295803 1.0000000
## hispanic:socst-asian:read 8.18920542 0.9997746
## white:socst-asian:read 12.34400946 1.0000000
## african-amer:write-asian:read 8.97842441 0.9999701
## asian:write-asian:read 20.50295803 0.9948838
## hispanic:write-asian:read 6.85587209 0.9910725
## white:write-asian:read 12.71642325 0.9999999
## white:read-hispanic:read 14.70583374 0.0666644
## african-amer:science-hispanic:read 6.13346644 0.9965795
## asian:science-hispanic:read 17.09450845 0.9981988
## hispanic:science-hispanic:read 9.41321208 1.0000000
## white:science-hispanic:read 14.94074824 0.0483562
## african-amer:socst-hispanic:read 13.01656468 0.9999903
## asian:socst-hispanic:read 16.63996300 0.9995335
## hispanic:socst-hispanic:read 10.88200326 1.0000000
## white:socst-hispanic:read 14.46445443 0.0941640
## african-amer:write-hispanic:read 11.76656468 1.0000000
## asian:write-hispanic:read 23.63996300 0.1165135
## hispanic:write-hispanic:read 9.54866992 1.0000000
## white:write-hispanic:read 14.83686822 0.0547934
## african-amer:science-white:read -3.25661317 0.0001517
## asian:science-white:read 8.10074927 0.9999992
## hispanic:science-white:read -0.12072992 0.0415606
## white:science-white:read 4.21418078 1.0000000
## african-amer:socst-white:read 3.58799670 0.9131063
## asian:socst-white:read 7.64620381 0.9999874
## hispanic:socst-white:read 1.31589121 0.2758883
## white:socst-white:read 3.72814138 1.0000000
## african-amer:write-white:read 2.33799670 0.5660173
## asian:write-white:read 14.64620381 0.9983974
## hispanic:write-white:read -0.01744212 0.0486727
## white:write-white:read 4.10055517 1.0000000
## asian:science-african-amer:science 21.83893435 0.5787297
## hispanic:science-african-amer:science 14.27464281 0.9993068
## white:science-african-amer:science 19.98339274 0.0000968
## african-amer:socst-african-amer:science 17.85693069 0.7294724
## asian:socst-african-amer:science 21.38438889 0.6757990
## hispanic:socst-african-amer:science 15.74969451 0.9537082
## white:socst-african-amer:science 19.50817812 0.0002518
## african-amer:write-african-amer:science 16.60693069 0.9382984
## asian:write-african-amer:science 28.38438889 0.0027826
## hispanic:write-african-amer:science 14.41636118 0.9982028
## white:write-african-amer:science 19.88059191 0.0001146
## hispanic:science-asian:science 7.15327278 0.9948762
## white:science-asian:science 13.27155404 0.9999967
## african-amer:socst-asian:science 10.68296986 1.0000000
## asian:socst-asian:science 13.95750349 1.0000000
## hispanic:socst-asian:science 8.64375087 0.9999606
## white:socst-asian:science 12.79855491 0.9999999
## african-amer:write-asian:science 9.43296986 0.9999962
## asian:write-asian:science 20.95750349 0.9881051
## hispanic:write-asian:science 7.31041754 0.9968890
## white:write-asian:science 13.17096870 0.9999981
## white:science-hispanic:science 15.52747553 0.0296894
## african-amer:socst-hispanic:science 13.56646421 0.9999162
## asian:socst-hispanic:science 17.17303563 0.9983755
## hispanic:socst-hispanic:science 11.43676280 1.0000000
## white:socst-hispanic:science 15.05138402 0.0598600
## african-amer:write-hispanic:science 12.31646421 1.0000000
## asian:write-hispanic:science 24.17303563 0.0857498
## hispanic:write-hispanic:science 10.10342947 1.0000000
## white:write-hispanic:science 15.42379782 0.0338258
## african-amer:socst-white:science 3.37456354 0.8722678
## asian:socst-white:science 7.43032485 0.9999607
## hispanic:socst-white:science 1.10330692 0.2188010
## white:socst-white:science 3.52530269 1.0000000
## african-amer:write-white:science 2.12456354 0.4915050
## asian:write-white:science 14.43032485 0.9992566
## hispanic:write-white:science -0.23002641 0.0347985
## white:write-white:science 3.89771648 1.0000000
## asian:socst-african-amer:socst 14.23751532 1.0000000
## hispanic:socst-african-amer:socst 8.57489802 1.0000000
## white:socst-african-amer:socst 12.29489325 0.9469295
## african-amer:write-african-amer:socst 9.43826156 1.0000000
## asian:write-african-amer:socst 21.23751532 0.6652797
## hispanic:write-african-amer:socst 7.24156468 0.9999701
## white:write-african-amer:socst 12.66730704 0.8897359
## hispanic:socst-asian:socst 9.09829633 0.9999951
## white:socst-asian:socst 13.25310037 0.9999968
## african-amer:write-asian:socst 9.88751532 0.9999997
## asian:write-asian:socst 21.41204894 0.9753607
## hispanic:write-asian:socst 7.76496300 0.9991084
## white:write-asian:socst 13.62551416 0.9999749
## white:socst-hispanic:socst 13.33945443 0.3499432
## african-amer:write-hispanic:socst 10.64156468 1.0000000
## asian:write-hispanic:socst 22.51496300 0.2629362
## hispanic:write-hispanic:socst 8.42366992 1.0000000
## white:write-hispanic:socst 13.71186822 0.2399616
## african-amer:write-white:socst 2.57937601 0.6486508
## asian:write-white:socst 14.88758312 0.9966444
## hispanic:write-white:socst 0.22393719 0.0699780
## white:write-white:socst 4.34193448 1.0000000
## asian:write-african-amer:write 22.48751532 0.3962937
## hispanic:write-african-amer:write 8.49156468 1.0000000
## white:write-african-amer:write 13.91730704 0.5207219
## hispanic:write-asian:write 0.76496300 0.0982451
## white:write-asian:write 6.62551416 0.9989626
## white:write-hispanic:write 15.04520155 0.0396354african-amer science:socst science:write顯著
- Perform all pairwise simple regressions for these variables: read, write, math, science, and socst.
## read write math science socst
## (Intercept) 18.4161841 21.1309429 5.627211e-14 24.3565993 21.1259457
## x 0.6507527 0.5843927 1.000000e+00 0.5455811 0.5935322z <- combn(c("read", "write", "math", "science", "socst"), 2, simplify = FALSE)
z1 <- matrix(unlist(z), ncol = 2, byrow = TRUE)[,1];z1## [1] "read" "read" "read" "read" "write" "write" "write"
## [8] "math" "math" "science"## [1] "write" "math" "science" "socst" "math" "science" "socst"
## [8] "science" "socst" "socst"## read read.1 read.2 read.3 write
## (Intercept) 23.9594444 21.0381578 20.6602830 18.4161841 20.4377529
## x 0.5517051 0.6051473 0.5998076 0.6507527 0.6102747
## write.1 write.2 math math.1 science
## (Intercept) 21.1309429 16.2535476 17.49497 19.5572360 26.1686196
## x 0.5843927 0.6850109 0.65494 0.6239484 0.5034689# method 3
results <- lapply(seq_along(z), function (n) {
df <- dta[,colnames(dta) %in% unlist(z[n])]
result <- lm(df[,1] ~ df[,2])
return(result)})
names(results) <- paste(z1, z2, sep = " ~ ")
results## $`read ~ write`
##
## Call:
## lm(formula = df[, 1] ~ df[, 2])
##
## Coefficients:
## (Intercept) df[, 2]
## 18.1622 0.6455
##
##
## $`read ~ math`
##
## Call:
## lm(formula = df[, 1] ~ df[, 2])
##
## Coefficients:
## (Intercept) df[, 2]
## 14.0725 0.7248
##
##
## $`read ~ science`
##
## Call:
## lm(formula = df[, 1] ~ df[, 2])
##
## Coefficients:
## (Intercept) df[, 2]
## 18.157 0.654
##
##
## $`read ~ socst`
##
## Call:
## lm(formula = df[, 1] ~ df[, 2])
##
## Coefficients:
## (Intercept) df[, 2]
## 21.1259 0.5935
##
##
## $`write ~ math`
##
## Call:
## lm(formula = df[, 1] ~ df[, 2])
##
## Coefficients:
## (Intercept) df[, 2]
## 19.8872 0.6247
##
##
## $`write ~ science`
##
## Call:
## lm(formula = df[, 1] ~ df[, 2])
##
## Coefficients:
## (Intercept) df[, 2]
## 24.3566 0.5456
##
##
## $`write ~ socst`
##
## Call:
## lm(formula = df[, 1] ~ df[, 2])
##
## Coefficients:
## (Intercept) df[, 2]
## 24.792 0.534
##
##
## $`math ~ science`
##
## Call:
## lm(formula = df[, 1] ~ df[, 2])
##
## Coefficients:
## (Intercept) df[, 2]
## 21.3917 0.6003
##
##
## $`math ~ socst`
##
## Call:
## lm(formula = df[, 1] ~ df[, 2])
##
## Coefficients:
## (Intercept) df[, 2]
## 27.7456 0.4751
##
##
## $`science ~ socst`
##
## Call:
## lm(formula = df[, 1] ~ df[, 2])
##
## Coefficients:
## (Intercept) df[, 2]
## 30.1197 0.4167Q2.
The formula P = L (r/(1-(1+r)^(-M)) describes the payment you have to make per month for M number of months if you take out a loan of L amount today at a monthly interest rate of r. Compute how much you will have to pay per month for 10, 15, 20, 25, or 30 years if you borrow NT$5,000,000, 10,000,000, or 15,000,000 from a bank that charges you 2%, 5%, or 7% for the monthly interest rate.
payPerMonth <- function(loan,monthlyInterestRate,year){
M <- 12 * year
p <- loan*(monthlyInterestRate/(1-(1+monthlyInterestRate)^(-M)))
# cat(" loan:",loan,"\n",
# "monthly interest rate:", monthlyInterestRate,"\n",
# "year:",year,"\n",
# "you have to pay $",p,"per month.\n\n")
return(p)
}
Loan <- c(5000000, 10000000, 15000000)
MonthlyInterestRate <- c(0.02, 0.05, 0.07)
Year <- c(10, 15, 20, 25, 30)
all <- data.frame(Loan = rep(Loan, each = 15),
Rate = rep(MonthlyInterestRate, each = 5, 3),
Year = rep(Year, 9))
all$payPerMonth <- payPerMonth(all$Loan, all$Rate, all$Year)
allQ3.
The following R script is an attempt to demonstrate the correspondence between parameter estimations by the least square method and the maximum likelihood method for the case of simple linear regression with a constant normal error term. 1. Construct a function from the script so that any deviance value for pairs of parameter estimates can be found. 2. Generalize the function further so that it will work with any data sets that can be modeled by a simple linear regression with a constant normal error term.
m0 <- lm(weight ~ height, data=women)
# summary(m0)
param <- c(coef(m0)[1], coef(m0)[2])
a <- param[1]
b <- param[2]
yhat <- a + b*women$height
e <- summary(m0)$sigma
lkhd <- dnorm(women$weight, mean=yhat, sd=e)
dvnc <- -2 * sum(log(lkhd))
ci_a <- coef(m0)[1] + unlist(summary(m0))$coefficients3*c(-2,2)
ci_b <- coef(m0)[2] + unlist(summary(m0))$coefficients4*c(-2,2)
bb <- expand.grid(a=seq(ci_a[1], ci_a[2], len=50),
b=seq(ci_b[1], ci_b[2], len=50))
## Not working yet
dvnc_fun <- function(x){
m0 <- lm(weight ~ height, data=women)
a <- x[1]
b <- x[2]
yhat <- a + b*women$height
e <- summary(m0)$sigma
lkhd <- dnorm(women$weight, mean=yhat, sd=e)
dvnc <- -2 * sum(log(lkhd))
return(dvnc)
}
bb$d <- apply(bb, 1, dvnc_fun)
head(bb)lmWithCIs <- function(m0){
ci_a <- coef(m0)[1] + unlist(summary(m0))$coefficients3*c(-2,2)
ci_b <- coef(m0)[2] + unlist(summary(m0))$coefficients4*c(-2,2)
print(paste0('95% CI intercept: ',min(ci_a),' ~ ',max(ci_a)))
print(paste0('95% CI slope : ',min(ci_b),' ~ ',max(ci_b)))
}
lmWithCIs(m0)## [1] "95% CI intercept: -99.3905547644725 ~ -75.6427785688604"
## [1] "95% CI slope : 3.26772700906759 ~ 3.6322729909324"Q4.
Modify this R script to create a function to compute the c-statistic illustrated with the data set in the article: Tryon, W.W. (1984). A simplified time-series analysis for evaluating treatment interventions. Journal of Applied Behavioral Analysis, 34(4), 230-233.
## 'data.frame': 42 obs. of 1 variable:
## $ nc: int 28 46 39 45 24 20 35 37 36 40 ...#
# plot figure 1
#
plot(1:42, dta[,1], xlab="Observations", ylab="Number of Children")
lines(1:42, dta[,1])
abline(v=10, lty=2)
abline(v=32, lty=2)
segments(1, mean(dta[1:10,1]),10, mean(dta[1:10,1]),col="red")
segments(11, mean(dta[11:32,1]),32, mean(dta[11:32,1]),col="red")
segments(33, mean(dta[33:42,1]),42, mean(dta[33:42,1]),col="red")
#
# calculate c-stat for first baseline phase
#
cden <- 1-(sum(diff(dta[1:10,1])^2)/(2*(10-1)*var(dta[1:10,1])))
sc <- sqrt((10-2)/((10-1)*(10+1)))
pval <- 1-pnorm(cden/sc)
pval## [1] 0.2866238#
# calculate c-stat for first baseline plus group tokens
#
n <- 32
cden <- 1-(sum(diff(dta[1:n,1])^2)/(2*(n-1)*var(dta[1:n,1])))
sc <- sqrt((n-2)/((n-1)*(n+1)))
pval <- 1-pnorm(cden/sc)
list(z=cden/sc,pvalue=pval)## $z
## [1] 3.879054
##
## $pvalue
## [1] 5.243167e-05cStat <- function(x){
n <- length(x)
cden <- 1-(sum(diff(x)^2)/(2*(10-1)*var(x)))
sc <- sqrt((n-2)/((n-1)*(n+1)))
pval <- 1-pnorm(cden/sc)
return(pval)
}
cStat(dta[1:10,1])## [1] 0.2866238GcStat <- function(x){
n <- length(x)
cden <- 1-(sum(diff(dta[1:n,1])^2)/(2*(n-1)*var(dta[1:n,1])))
sc <- sqrt((n-2)/((n-1)*(n+1)))
pval <- 1-pnorm(cden/sc)
return(list(z=cden/sc,pvalue=pval))
}
GcStat(dta[1:32,1])## $z
## [1] 3.879054
##
## $pvalue
## [1] 5.243167e-05Q5.
Plot the likelihood functoion to estimate the probability of graduate admission by gender, respectively, for Dept A in UCBAdmissions{datasets}. Construct approximate 95%-CI for each gender. Do they overlap?
dta <- UCBAdmissions
n <- 100
n_M <- sum(dta[,1,1])
n_F <- sum(dta[,2,1])
p_M <- dta[1,1,1] / n_M
p_F <- dta[1,2,1] / n_F
y_M <- rbinom(n, 1, p_M)
y_F <- rbinom(n, 1, p_F)
theta <- seq(0.01, 0.99, by=.01)
llkhd_M <- sum(y_M)*log(theta)+(n-sum(y_M))*log(1-theta)
llkhd_F <- sum(y_F)*log(theta)+(n-sum(y_F))*log(1-theta)
aggr <- data.frame(theta = theta, Male = llkhd_M, Female = llkhd_F)ggplot(aes(x=theta,y=Male), data = aggr) +
geom_line(aes(x=theta,y=Male), color = 'orange') +
geom_line(aes(x=theta,y=Female), color = 'purple') +
labs(y='Likehood',title='Grid Research') +
geom_vline(xintercept = mean(y_M), color = 'orange', linetype = 2) +
geom_point(aes(mean(y_M),-200), color = 'orange', size = 3) +
geom_errorbar(aes(
xmin=mean(y_M) - 2*sqrt(mean(y_M)*(1-mean(y_M)))/sqrt(n_M),
xmax=mean(y_M) + 2*sqrt(mean(y_M)*(1-mean(y_M)))/sqrt(n_M),
y=-200), width=.2, size=.4, color = 'orange') +
geom_vline(xintercept = mean(y_F), color = 'purple', linetype = 2) +
geom_point(aes(mean(y_F),-200), color = 'purple', size = 3) +
geom_errorbar(aes(
xmin=mean(y_F) - 2*sqrt(mean(y_F)*(1-mean(y_F)))/sqrt(n_F),
xmax=mean(y_F) + 2*sqrt(mean(y_F)*(1-mean(y_F)))/sqrt(n_F),
y=-200), width=.2, size=.4, color = 'purple')
Q6.
(Bonus) The data set contains inter-response times (in milliseconds) in the resting activity of a single neuron recorded from the spinal cord of a cat. Write a function to fit an exponential distribution to the data. More specifically, estimate the rate parameter of the exponential distribution using the maximum likelihood method. Source: McGill, W.J. (1963). Luce, Bush, & Galanter, eds. Handbook of Mathematical Psychology.
log_like_lambda <- function(lambda, x){
length(x) * log(lambda) - lambda*length(x)*mean(x)
}
optimize(f = log_like_lambda, dta$RT, interval = c(1,15), maximum = TRUE)## $maximum
## [1] 1.000043
##
## $objective
## [1] -15988.67