using_dummy_variables
This is an old revision of the document!
Table of Contents
Categorical variables
2 groups
rscr01
rm(list=ls())
library(dplyr)
library(ggplot2)
ss <- function(x) {
sum((x-mean(x))^2)
}
sp <- function(x, y) {
sum((x-mean(x))*(y-mean(y)))
}
# 1. Generate numerical variable
set.seed(12)
n <- 30
n.1 <- round(rnorm(n, mean = 50, sd = 20),0) # Continuous variable
n.2 = round(rnorm(n, 55, 20),0)
n.3 = round(rnorm(n, 60, 22),0)
sc <- data.frame(n.1,n.2,n.3)
sc.s <- stack(sc)
sc.s <- sc.s[-2]
head(sc.s)
c.1 <- (rep("lo", n))
c.2 <- (rep("mid", n))
c.3 <- rep("hi", n)
ca <- data.frame(c.1,c.2,c.3)
head(ca)
ca.s <- stack(ca)
ca.s <- ca.s[-2]
head(ca.s)
df <- data.frame(sc.s, ca.s)
colnames(df) <-c("score", "group")
df$group <- factor(df$group)
str(df)
head(df)
amean <- aggregate(df[, 1], list(df$group), mean)
colnames(amean)=c("group", "mean")
amean
boxplot(df$score~df$group)
# 세 그룹 간의 차이를 보려면
# aov 테스트를 (F-test) 하면
# 된다.
fm.df <- aov(score~group, data=df)
summary(fm.df)
# 이것을 regression을 해도 된다
m.df <- lm(score~group, data=df)
summary(m.df)
anova(m.df)
summary(fm.df)
boxplot(df$score~df$group)
amean
# 2. Manually create dummy variables (k - 1 = 3 variables needed)
df$lo <- ifelse(df$group == "lo", 1, 0)
df$mid <- ifelse(df$group == "mid", 1, 0)
df$hi <- ifelse(df$group == "hi", 1, 0)
head(df)
m.tmp <- lm(score~lo+mid, data=df)
summary(m.tmp)
summary(m.df)
m.tmp2 <- lm(score~mid+hi, data=df)
summary(m.tmp2)
b.hi <- m.tmp$coefficients[1]
b.lo <- m.tmp$coefficients[2]
b.mid <- m.tmp$coefficients[3]
b.lo
b.mid
b.hi
mean.hi <- b.hi
mean.mid <- mean.hi+b.mid
mean.lo <- mean.hi+b.lo
data.frame(mean.hi, mean.lo, mean.mid)
amean
am.df <- aov(score~group, df)
summary(am.df)
posthoc.am.df <- TukeyHSD(am.df)
posthoc.am.df
# t-test의 의미 = baseline과 (절편) 차이가 있는가 테스트
# basedline의 t-test는 baseline이 0과 다른가?
# 다른 조건과 차이를 테스트하는 것은 아님 (주의)
summary(m.df)
# 위에서 grouplo의 significance는 lo 그룹과 hi 그룹이
# (절편, baseline) 차이가 있는가를 테스트한 것
# 마찬가지로 groupmid의 t-test는 baseline인 grouphi와
# 차이가 있는지 본 것으로 not significant
# 여기서 보지 못한 것은 mid와 lo그룹 간의 차이가 있는지를
# 보는 것이다
# 이를 위해서는 아래를 보면 된다
summary(m.tmp2)
# 위의 결과는 baseline인 lo가 mid와는 차이가 없는 것을
# 나타낸다. lo와 hi는 존재
# 위는 두개의 regression을 해서 그룹 간 차이를
# post hoc test처럼 한 것이고 이를 한 번에 처리하는
# 방법은
# install.packages("emmeans")
library(emmeans)
# Calculate estimated marginal means for each group
gmeans <- emmeans(m.df, specs = ~ group)
# Perform pairwise comparisons using Tukey's adjustment (Default)
pairs(gmeans, adjust = "tukey")
rout01
> rm(list=ls())
>
> library(dplyr)
> library(ggplot2)
>
> ss <- function(x) {
+ sum((x-mean(x))^2)
+ }
> sp <- function(x, y) {
+ sum((x-mean(x))*(y-mean(y)))
+ }
>
> # 1. Generate numerical variable
> set.seed(12)
> n <- 30
> n.1 <- round(rnorm(n, mean = 50, sd = 20),0) # Continuous variable
> n.2 = round(rnorm(n, 55, 20),0)
> n.3 = round(rnorm(n, 60, 22),0)
> sc <- data.frame(n.1,n.2,n.3)
> sc.s <- stack(sc)
> sc.s <- sc.s[-2]
> head(sc.s)
values
1 20
2 82
3 31
4 32
5 10
6 45
>
> c.1 <- (rep("lo", n))
> c.2 <- (rep("mid", n))
> c.3 <- rep("hi", n)
> ca <- data.frame(c.1,c.2,c.3)
> head(ca)
c.1 c.2 c.3
1 lo mid hi
2 lo mid hi
3 lo mid hi
4 lo mid hi
5 lo mid hi
6 lo mid hi
> ca.s <- stack(ca)
> ca.s <- ca.s[-2]
> head(ca.s)
values
1 lo
2 lo
3 lo
4 lo
5 lo
6 lo
>
> df <- data.frame(sc.s, ca.s)
> colnames(df) <-c("score", "group")
> df$group <- factor(df$group)
> str(df)
'data.frame': 90 obs. of 2 variables:
$ score: num 20 82 31 32 10 45 44 37 48 59 ...
$ group: Factor w/ 3 levels "hi","lo","mid": 2 2 2 2 2 2 2 2 2 2 ...
>
> head(df)
score group
1 20 lo
2 82 lo
3 31 lo
4 32 lo
5 10 lo
6 45 lo
> amean <- aggregate(df[, 1], list(df$group), mean)
> colnames(amean)=c("group", "mean")
> amean
group mean
1 hi 62.33333
2 lo 46.96667
3 mid 54.10000
> boxplot(df$score~df$group)
>
> # 세 그룹 간의 차이를 보려면
> # aov 테스트를 (F-test) 하면
> # 된다.
> fm.df <- aov(score~group, data=df)
> summary(fm.df)
Df Sum Sq Mean Sq F value Pr(>F)
group 2 3548 1774.0 5.278 0.00686 **
Residuals 87 29242 336.1
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>
> # 이것을 regression을 해도 된다
> m.df <- lm(score~group, data=df)
> summary(m.df)
Call:
lm(formula = score ~ group, data = df)
Residuals:
Min 1Q Median 3Q Max
-49.333 -11.242 -1.533 12.000 43.033
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 62.333 3.347 18.622 < 2e-16 ***
grouplo -15.367 4.734 -3.246 0.00166 **
groupmid -8.233 4.734 -1.739 0.08552 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 18.33 on 87 degrees of freedom
Multiple R-squared: 0.1082, Adjusted R-squared: 0.0877
F-statistic: 5.278 on 2 and 87 DF, p-value: 0.006863
> anova(m.df)
Analysis of Variance Table
Response: score
Df Sum Sq Mean Sq F value Pr(>F)
group 2 3548.1 1774.03 5.278 0.006863 **
Residuals 87 29242.3 336.12
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> summary(fm.df)
Df Sum Sq Mean Sq F value Pr(>F)
group 2 3548 1774.0 5.278 0.00686 **
Residuals 87 29242 336.1
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>
> boxplot(df$score~df$group)
> amean
group mean
1 hi 62.33333
2 lo 46.96667
3 mid 54.10000
>
> # 2. Manually create dummy variables (k - 1 = 3 variables needed)
> df$lo <- ifelse(df$group == "lo", 1, 0)
> df$mid <- ifelse(df$group == "mid", 1, 0)
> df$hi <- ifelse(df$group == "hi", 1, 0)
>
> head(df)
score group lo mid hi
1 20 lo 1 0 0
2 82 lo 1 0 0
3 31 lo 1 0 0
4 32 lo 1 0 0
5 10 lo 1 0 0
6 45 lo 1 0 0
> m.tmp <- lm(score~lo+mid, data=df)
> summary(m.tmp)
Call:
lm(formula = score ~ lo + mid, data = df)
Residuals:
Min 1Q Median 3Q Max
-49.333 -11.242 -1.533 12.000 43.033
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 62.333 3.347 18.622 < 2e-16 ***
lo -15.367 4.734 -3.246 0.00166 **
mid -8.233 4.734 -1.739 0.08552 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 18.33 on 87 degrees of freedom
Multiple R-squared: 0.1082, Adjusted R-squared: 0.0877
F-statistic: 5.278 on 2 and 87 DF, p-value: 0.006863
> summary(m.df)
Call:
lm(formula = score ~ group, data = df)
Residuals:
Min 1Q Median 3Q Max
-49.333 -11.242 -1.533 12.000 43.033
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 62.333 3.347 18.622 < 2e-16 ***
grouplo -15.367 4.734 -3.246 0.00166 **
groupmid -8.233 4.734 -1.739 0.08552 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 18.33 on 87 degrees of freedom
Multiple R-squared: 0.1082, Adjusted R-squared: 0.0877
F-statistic: 5.278 on 2 and 87 DF, p-value: 0.006863
> m.tmp2 <- lm(score~mid+hi, data=df)
> summary(m.tmp2)
Call:
lm(formula = score ~ mid + hi, data = df)
Residuals:
Min 1Q Median 3Q Max
-49.333 -11.242 -1.533 12.000 43.033
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 46.967 3.347 14.031 < 2e-16 ***
mid 7.133 4.734 1.507 0.13545
hi 15.367 4.734 3.246 0.00166 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 18.33 on 87 degrees of freedom
Multiple R-squared: 0.1082, Adjusted R-squared: 0.0877
F-statistic: 5.278 on 2 and 87 DF, p-value: 0.006863
>
> b.hi <- m.tmp$coefficients[1]
> b.lo <- m.tmp$coefficients[2]
> b.mid <- m.tmp$coefficients[3]
> b.lo
lo
-15.36667
> b.mid
mid
-8.233333
> b.hi
(Intercept)
62.33333
>
> mean.hi <- b.hi
> mean.mid <- mean.hi+b.mid
> mean.lo <- mean.hi+b.lo
> data.frame(mean.hi, mean.lo, mean.mid)
mean.hi mean.lo mean.mid
(Intercept) 62.33333 46.96667 54.1
> amean
group mean
1 hi 62.33333
2 lo 46.96667
3 mid 54.10000
>
> am.df <- aov(score~group, df)
> summary(am.df)
Df Sum Sq Mean Sq F value Pr(>F)
group 2 3548 1774.0 5.278 0.00686 **
Residuals 87 29242 336.1
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> posthoc.am.df <- TukeyHSD(am.df)
> posthoc.am.df
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = score ~ group, data = df)
$group
diff lwr upr p adj
lo-hi -15.366667 -26.654078 -4.079255 0.0046872
mid-hi -8.233333 -19.520745 3.054078 0.1965275
mid-lo 7.133333 -4.154078 18.420745 0.2926865
>
> # t-test의 의미 = baseline과 (절편) 차이가 있는가 테스트
> # basedline의 t-test는 baseline이 0과 다른가?
> # 다른 조건과 차이를 테스트하는 것은 아님 (주의)
> summary(m.df)
Call:
lm(formula = score ~ group, data = df)
Residuals:
Min 1Q Median 3Q Max
-49.333 -11.242 -1.533 12.000 43.033
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 62.333 3.347 18.622 < 2e-16 ***
grouplo -15.367 4.734 -3.246 0.00166 **
groupmid -8.233 4.734 -1.739 0.08552 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 18.33 on 87 degrees of freedom
Multiple R-squared: 0.1082, Adjusted R-squared: 0.0877
F-statistic: 5.278 on 2 and 87 DF, p-value: 0.006863
> # 위에서 grouplo의 significance는 lo 그룹과 hi 그룹이
> # (절편, baseline) 차이가 있는가를 테스트한 것
> # 마찬가지로 groupmid의 t-test는 baseline인 grouphi와
> # 차이가 있는지 본 것으로 not significant
> # 여기서 보지 못한 것은 mid와 lo그룹 간의 차이가 있는지를
> # 보는 것이다
> # 이를 위해서는 아래를 보면 된다
> summary(m.tmp2)
Call:
lm(formula = score ~ mid + hi, data = df)
Residuals:
Min 1Q Median 3Q Max
-49.333 -11.242 -1.533 12.000 43.033
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 46.967 3.347 14.031 < 2e-16 ***
mid 7.133 4.734 1.507 0.13545
hi 15.367 4.734 3.246 0.00166 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 18.33 on 87 degrees of freedom
Multiple R-squared: 0.1082, Adjusted R-squared: 0.0877
F-statistic: 5.278 on 2 and 87 DF, p-value: 0.006863
> # 위의 결과는 baseline인 lo가 mid와는 차이가 없는 것을
> # 나타낸다. lo와 hi는 존재
> # 위는 두개의 regression을 해서 그룹 간 차이를
> # post hoc test처럼 한 것이고 이를 한 번에 처리하는
> # 방법은
>
> # install.packages("emmeans")
> library(emmeans)
> # Calculate estimated marginal means for each group
> gmeans <- emmeans(m.df, specs = ~ group)
> # Perform pairwise comparisons using Tukey's adjustment (Default)
> pairs(gmeans, adjust = "tukey")
contrast estimate SE df t.ratio p.value
hi - lo 15.37 4.73 87 3.246 0.0047
hi - mid 8.23 4.73 87 1.739 0.1965
lo - mid -7.13 4.73 87 -1.507 0.2927
P value adjustment: tukey method for comparing a family of 3 estimates
2 group: e.g. 2
data: for spss
elemapi2.sav
elemapi2_categories.sps
data: for r
elemapi2.csv
rscrp.02
####################
# 다른 예
####################
rm(list=ls())
# library(dplyr)
library(ggplot2)
ss <- function(x) {
sum((x-mean(x))^2)
}
sp <- function(x, y) {
sum((x-mean(x))*(y-mean(y)))
}
datavar <- read.csv("http://commres.net/_media/r/elemapi2.csv")
head(datavar)
rout.02
> ####################
> # 다른 예
> ####################
> rm(list=ls())
>
> # library(dplyr)
> library(ggplot2)
>
> ss <- function(x) {
+ sum((x-mean(x))^2)
+ }
> sp <- function(x, y) {
+ sum((x-mean(x))*(y-mean(y)))
+ }
>
> datavar <- read.csv("http://commres.net/_media/r/elemapi2.csv")
> head(datavar)
snum dnum api00 api99 growth meals ell yr_rnd mobility acs_k3 acs_46 not_hsg hsg some_col
1 906 41 693 600 93 67 9 0 11 16 22 0 0 0
2 889 41 570 501 69 92 21 0 33 15 32 0 0 0
3 887 41 546 472 74 97 29 0 36 17 25 0 0 0
4 876 41 571 487 84 90 27 0 27 20 30 36 45 9
5 888 41 478 425 53 89 30 0 44 18 31 50 50 0
6 4284 98 858 844 14 10 3 0 10 20 33 1 8 24
col_grad grad_sch avg_ed full emer enroll mealcat collcat
1 0 0 NA 76 24 247 2 1
2 0 0 NA 79 19 463 3 1
3 0 0 NA 68 29 395 3 1
4 9 0 1.91 87 11 418 3 1
5 0 0 1.50 87 13 520 3 1
6 36 31 3.89 100 0 343 1 2
>
datavar data set의 변인 설명
Variable Labels Variable Position Label snum 1 school number dnum 2 district number api00 3 api 2000 api99 4 api 1999 growth 5 growth 1999 to 2000 meals 6 pct free meals ell 7 english language learners yr_rnd 8 year round school mobility 9 pct 1st year in school acs_k3 10 avg class size k-3 acs_46 11 avg class size 4-6 not_hsg 12 parent not hsg hsg 13 parent hsg some_col 14 parent some college col_grad 15 parent college grad grad_sch 16 parent grad school avg_ed 17 avg parent ed full 18 pct full credential emer 19 pct emer credential enroll 20 number of students mealcat 21 Percentage free meals in 3 categories collcat 22 <none> Variables in the working file yr_rnd: 0 = 방학있음 1 = 방학없음 mealcat: 1 = 0-46% free meals 2 = 47-80 3 = 81-100
위의 변인들 중에서 “무방학학교”가 성적에 어떤 영향을 미칠 것인가를 알아 보기 위해서 regression 테스트를 시행하였다. 아래는 그 결과이다.
rscrp.03
# data 정리
df <- datavar[,c(1,3,6,7,8,9,18,21)]
head(df)
str(df)
df$yr_rnd <- factor(df$yr_rnd, levels=c(0,1), labels=c("break", "nobr"))
df$mealcat <- factor(df$mealcat, levels=c(1:3), labels=c("to46", "to80", "to100"))
str(df)
# regression
m.br <- lm(api00~yr_rnd, data=df)
summary(m.br)
# making dummy variables
df$nobr <- ifelse(df$yr_rnd == "nobr", 1, 0)
df$br<- ifelse(df$yr_rnd == "break", 1, 0)
head(df)
m.br2 <- lm(api00~nobr,data=df)
summary(m.br2)
m.br3 <- lm(api00~br, data=df)
summary(m.br3)
# baseline 변화
df$yr_rnd <- relevel(factor(df$yr_rnd), ref = "nobr")
m.br4 <- lm(api00 ~ yr_rnd, data = df)
summary(m.br4)
# 원래로 복원
df$yr_rnd <- relevel(factor(df$yr_rnd), ref = "break")
head(df)
tmp <- aggregate(df[, 2], list(df$yr_rnd), mean)
mean.br <- tmp[[2]][1]
mean.nobr <- tmp[[2]][2]
data.frame(mean.br, mean.nobr)
boxplot(df$api00~df$yr_rnd)
t.out <- t.test(df$api00~df$yr_rnd, var.eq=T)
t.out
t.out$estimate
# str(t.out)
t.value <- t.out$statistic
f.out <- aov(df$api00~df$yr_rnd, data=df)
summary(f.out)
f.value <- summary(f.out)[[1]][1, "F value"]
summary(m.br)
f.value.lm <- summary(m.br)$fstatistic[1]
t.value.lm <- summary(m.br)$coefficients["yr_rndnobr", "t value"]
data.frame(t.value, t.value^2, f.value, f.value.lm, t.value.lm)
plot(df$yr_rnd,df$api00)
df %>%
ggplot(aes(x = yr_rnd, y = api00, color = factor(yr_rnd))) +
geom_point(size = 3) +
geom_smooth(method = "lm", se = FALSE) +
scale_color_brewer(palette = "Set1", name = "breaks") +
labs(title = "Regression: with yr_rnd break",
x = " (yr_rnd)",
y = "api00")
summary(m.br)
rout.03
> # data 정리
> df <- datavar[,c(1,3,6,7,8,9,18,21)]
> head(df)
snum api00 meals ell yr_rnd mobility full mealcat
1 906 693 67 9 0 11 76 2
2 889 570 92 21 0 33 79 3
3 887 546 97 29 0 36 68 3
4 876 571 90 27 0 27 87 3
5 888 478 89 30 0 44 87 3
6 4284 858 10 3 0 10 100 1
> str(df)
'data.frame': 400 obs. of 8 variables:
$ snum : int 906 889 887 876 888 4284 4271 2910 2899 2887 ...
$ api00 : int 693 570 546 571 478 858 918 831 860 737 ...
$ meals : int 67 92 97 90 89 10 5 2 5 29 ...
$ ell : int 9 21 29 27 30 3 2 3 6 15 ...
$ yr_rnd : int 0 0 0 0 0 0 0 0 0 0 ...
$ mobility: int 11 33 36 27 44 10 16 44 10 17 ...
$ full : int 76 79 68 87 87 100 100 96 100 96 ...
$ mealcat : int 2 3 3 3 3 1 1 1 1 1 ...
> df$yr_rnd <- factor(df$yr_rnd, levels=c(0,1), labels=c("break", "nobr"))
> df$mealcat <- factor(df$mealcat, levels=c(1:3), labels=c("to46", "to80", "to100"))
> str(df)
'data.frame': 400 obs. of 8 variables:
$ snum : int 906 889 887 876 888 4284 4271 2910 2899 2887 ...
$ api00 : int 693 570 546 571 478 858 918 831 860 737 ...
$ meals : int 67 92 97 90 89 10 5 2 5 29 ...
$ ell : int 9 21 29 27 30 3 2 3 6 15 ...
$ yr_rnd : Factor w/ 2 levels "break","nobr": 1 1 1 1 1 1 1 1 1 1 ...
$ mobility: int 11 33 36 27 44 10 16 44 10 17 ...
$ full : int 76 79 68 87 87 100 100 96 100 96 ...
$ mealcat : Factor w/ 3 levels "to46","to80",..: 2 3 3 3 3 1 1 1 1 1 ...
>
> # regression
> m.br <- lm(api00~yr_rnd, data=df)
> summary(m.br)
Call:
lm(formula = api00 ~ yr_rnd, data = df)
Residuals:
Min 1Q Median 3Q Max
-273.539 -95.662 0.967 103.341 297.967
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 684.54 7.14 95.88 <2e-16 ***
yr_rndnobr -160.51 14.89 -10.78 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 125.3 on 398 degrees of freedom
Multiple R-squared: 0.226, Adjusted R-squared: 0.2241
F-statistic: 116.2 on 1 and 398 DF, p-value: < 2.2e-16
>
> # making dummy variables
> df$nobr <- ifelse(df$yr_rnd == "nobr", 1, 0)
> df$br<- ifelse(df$yr_rnd == "break", 1, 0)
> head(df)
snum api00 meals ell yr_rnd mobility full mealcat nobr br
1 906 693 67 9 break 11 76 to80 0 1
2 889 570 92 21 break 33 79 to100 0 1
3 887 546 97 29 break 36 68 to100 0 1
4 876 571 90 27 break 27 87 to100 0 1
5 888 478 89 30 break 44 87 to100 0 1
6 4284 858 10 3 break 10 100 to46 0 1
>
> m.br2 <- lm(api00~nobr,data=df)
> summary(m.br2)
Call:
lm(formula = api00 ~ nobr, data = df)
Residuals:
Min 1Q Median 3Q Max
-273.539 -95.662 0.967 103.341 297.967
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 684.54 7.14 95.88 <2e-16 ***
nobr -160.51 14.89 -10.78 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 125.3 on 398 degrees of freedom
Multiple R-squared: 0.226, Adjusted R-squared: 0.2241
F-statistic: 116.2 on 1 and 398 DF, p-value: < 2.2e-16
>
> m.br3 <- lm(api00~br, data=df)
> summary(m.br3)
Call:
lm(formula = api00 ~ br, data = df)
Residuals:
Min 1Q Median 3Q Max
-273.539 -95.662 0.967 103.341 297.967
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 524.03 13.06 40.11 <2e-16 ***
br 160.51 14.89 10.78 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 125.3 on 398 degrees of freedom
Multiple R-squared: 0.226, Adjusted R-squared: 0.2241
F-statistic: 116.2 on 1 and 398 DF, p-value: < 2.2e-16
>
> # baseline 변화
> df$yr_rnd <- relevel(factor(df$yr_rnd), ref = "nobr")
> m.br4 <- lm(api00 ~ yr_rnd, data = df)
> summary(m.br4)
Call:
lm(formula = api00 ~ yr_rnd, data = df)
Residuals:
Min 1Q Median 3Q Max
-273.539 -95.662 0.967 103.341 297.967
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 524.03 13.06 40.11 <2e-16 ***
yr_rndbreak 160.51 14.89 10.78 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 125.3 on 398 degrees of freedom
Multiple R-squared: 0.226, Adjusted R-squared: 0.2241
F-statistic: 116.2 on 1 and 398 DF, p-value: < 2.2e-16
> # 원래로 복원
> df$yr_rnd <- relevel(factor(df$yr_rnd), ref = "break")
>
> head(df)
snum api00 meals ell yr_rnd mobility full mealcat nobr br
1 906 693 67 9 break 11 76 to80 0 1
2 889 570 92 21 break 33 79 to100 0 1
3 887 546 97 29 break 36 68 to100 0 1
4 876 571 90 27 break 27 87 to100 0 1
5 888 478 89 30 break 44 87 to100 0 1
6 4284 858 10 3 break 10 100 to46 0 1
> tmp <- aggregate(df[, 2], list(df$yr_rnd), mean)
> mean.br <- tmp[[2]][1]
> mean.nobr <- tmp[[2]][2]
> data.frame(mean.br, mean.nobr)
mean.br mean.nobr
1 684.539 524.0326
>
> boxplot(df$api00~df$yr_rnd)
>
> t.out <- t.test(df$api00~df$yr_rnd, var.eq=T)
> t.out
Two Sample t-test
data: df$api00 by df$yr_rnd
t = 10.782, df = 398, p-value < 2.2e-16
alternative hypothesis: true difference in means between group break and group nobr is not equal to 0
95 percent confidence interval:
131.2390 189.7737
sample estimates:
mean in group break mean in group nobr
684.5390 524.0326
> t.out$estimate
mean in group break mean in group nobr
684.5390 524.0326
> # str(t.out)
> t.value <- t.out$statistic
>
> f.out <- aov(df$api00~df$yr_rnd, data=df)
> summary(f.out)
Df Sum Sq Mean Sq F value Pr(>F)
df$yr_rnd 1 1825001 1825001 116.2 <2e-16 ***
Residuals 398 6248671 15700
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> f.value <- summary(f.out)[[1]][1, "F value"]
>
> summary(m.br)
Call:
lm(formula = api00 ~ yr_rnd, data = df)
Residuals:
Min 1Q Median 3Q Max
-273.539 -95.662 0.967 103.341 297.967
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 684.54 7.14 95.88 <2e-16 ***
yr_rndnobr -160.51 14.89 -10.78 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 125.3 on 398 degrees of freedom
Multiple R-squared: 0.226, Adjusted R-squared: 0.2241
F-statistic: 116.2 on 1 and 398 DF, p-value: < 2.2e-16
> f.value.lm <- summary(m.br)$fstatistic[1]
> t.value.lm <- summary(m.br)$coefficients["yr_rndnobr", "t value"]
>
> data.frame(t.value, t.value^2, f.value, f.value.lm, t.value.lm)
t.value t.value.2 f.value f.value.lm t.value.lm
t 10.7815 116.2407 116.2407 116.2407 -10.7815
>
> plot(df$yr_rnd,df$api00)
>
> df %>%
+ ggplot(aes(x = yr_rnd, y = api00, color = factor(yr_rnd))) +
+ geom_point(size = 3) +
+ geom_smooth(method = "lm", se = FALSE) +
+ scale_color_brewer(palette = "Set1", name = "breaks") +
+ labs(title = "Regression: with yr_rnd break",
+ x = " (yr_rnd)",
+ y = "api00")
`geom_smooth()` using formula = 'y ~ x'
>
> summary(m.br)
Call:
lm(formula = api00 ~ yr_rnd, data = df)
Residuals:
Min 1Q Median 3Q Max
-273.539 -95.662 0.967 103.341 297.967
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 684.54 7.14 95.88 <2e-16 ***
yr_rndnobr -160.51 14.89 -10.78 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 125.3 on 398 degrees of freedom
Multiple R-squared: 0.226, Adjusted R-squared: 0.2241
F-statistic: 116.2 on 1 and 398 DF, p-value: < 2.2e-16
>
위의 아웃풋을 살펴 보면, 학생들의 성적이 가지는 총 변량의 (Sum of Square Total) 약 22.6% 를 방학이 있고 없고로 구분되는 yr_rnd 변인이 설명을 하고 있으며, 이는 통계적으로 유의미한 것이다 (F(1, 398) = 116.241, p < .001). 위의 regression output은 yr_rnd 변인 중에서 방학이 있는 특성이 baseline이 되어 있으며, 이를 염두에 두고 regression식을 적어 보면 아래와 같다.
$\hat{\text{api00}} = 684.54 - 160.51 * \text{yr_rndnobr} $
이 때, $\text{yr_rndnobr} $ 은 no break인 경우를 1로 보고 대입된 dummy variable이다. 즉,
- X: 0 = break
- X: 1 = no break
이므로 x=0 일때를 대입해 보면, 즉, 방학이있는학교의 경우는 684.54 의 추정치를 엊을 수 있다. 이 점수는 방학이 있는 학교의 api점수 평균이 된다. 반면에 x = 1 일때를 대입해 보면 즉, 방학이없는학교의 경우에는 684.54 - 160.51 = 524.03 의 추정치를 엊을 수 있다. b coefficient가 이 역할 (차이를 나타내는 역할을 하는데)에 대한 유의성에 대한 판단은 t-test로 하는데, 이 t값은 -10.782며 이는 F값인 116.241의 제곱근이다 (즉, $t^2 = F$ ). 사실, 이 상황은 정확히 t-test를 해야할 상황이므로 (두 그룹에 대한 성적평균의 차이), t-test를 해야 하지만 이와 같이 regression을 하여도 동일한 결과를 보게된다 (같은 의미에서 F-test를 했어도 마찬가지). 뒷 부분의 t-test와 F-test는 이를 확인한 것.
위의 그래프는 두 그룹의 (break, no break 학교그룹) api 점수 분포를 나타내 주는 것이다. 그리고 직선은 두 그룹 평균을 연결한 선으로 $\hat{Y} = 684.539 - 160.506 X$ 으로 표현된다.
이와 같이 종류변인(category, nominal)을 가지고서도 regression 테스트를 할 수 있으며, 사실 이는 t-test나 F-test와 다르지 않다. 위에서 주의해야 할 점은 두 변인의 종류를 손으로 coding할 때, 1과 2가 아닌, 0과 1로 하였다는 점이다. 이렇게 하는 이유는 해석하기에 편하기 때문이며, 이것이 보통의 방법이다. 그러나, 1과 2로 coding 데이터를 이용해도 크게 다른지 않은 결과를 구하게 된다. 다른 점이라면, 절편에 해당되는 상수값이 다르게 되며, coefficient값은 위의 분석과 동일한 값을 갖게 된다. 그리고 일반적으로 r에서는 dummy variables을 얻기 위한 coding을 하지 않아도 r이 스스로 구하여 계산을 한다. 즉, dummy variable 을 구하지 않은 채로 lm(api00 ~ yr_rnd, data=df) 만으로도 계산이 된다.
Regression with a continuous and a categorical IV
rs.04
#######################################
#######################################
# with a continuous and a category IV
#######################################
#######################################
m.mealsyr_rnd <- lm(api00~meals+yr_rnd, data=df)
summary(m.mealsyr_rnd)
m.mealsxyr_rnd <- lm(api00~meals*yr_rnd, data=df)
summary(m.mealsxyr_rnd)
# 중요: 위에서 yr_rnd의 설명력이 t test를 보면
# 사라짐. 왜?
# interaction을 추가하면서 yr_rnd의 설명력은
# (방학이 있고 없음의 설명력은) meals의 percentage가
# 0 이 될때의 것으로 설정된다. 0이 되는 데이터가
# 극단치이고 이 경우에 nobr인 학교도 없으므로 계산된
# se값이 극대화되어 t값이 작아지게 된다. 이를 바로
# 잡으려면 meals의 0이 중간에 가도록 한 후에
# regression을 하면 meals가 0 일때가 극단치가 되지
# 않으므로 결과가 옳게 나온다
df %>%
ggplot(aes(x = meals, y = api00, color = factor(yr_rnd))) +
geom_point(size = 3) +
geom_smooth(method = "lm", se = FALSE) +
scale_color_brewer(palette = "Set1", name = "yr_rnd") +
labs(title = "Multiple Regression: Interaction by yr_rnd break",
x = " (free meals percentage)",
y = "api00")
# 1. 평균을 빼준 값을 새로운 변인으로 저장
df$meals_centered <-
df$meals - mean(df$meals)
# 2. Run the model again using the centered variable
m_centered <- lm(api00 ~ meals_centered * yr_rnd,
data = df)
# 3. Check the new summary
summary(m_centered)
df %>%
ggplot(aes(x = meals_centered, y = api00, color = factor(yr_rnd))) +
geom_point(size = 3) +
geom_smooth(method = "lm", se = FALSE) +
scale_color_brewer(palette = "Set1", name = "yr_rnd") +
labs(title = "Multiple Regression: Interaction by yr_rnd break",
x = " (free meals percentage CENTERED)",
y = "api00")
# Install the package if you do not have it
# install.packages("interactions")
# Plot the interaction
library(interactions)
interact_plot(m.mealsxyr_rnd, pred = meals, modx = yr_rnd)
m.ellyr_rnd <- lm(api00~ell+yr_rnd, data=df)
summary(m.ellyr_rnd)
m.ellxyr_rnd <- lm(api00~ell*yr_rnd, data=df)
summary(m.ellxyr_rnd)
df %>%
ggplot(aes(x = ell, y = api00, color = factor(yr_rnd))) +
geom_point(size = 3) +
geom_smooth(method = "lm", se = FALSE) +
scale_color_brewer(palette = "Set1", name = "Cylinders") +
labs(title = "Multiple Regression: Interaction by yr_rnd break",
x = " (ell, english language learner)",
y = "api00")
# ell mean centred 이번 경우는 큰 차이를 보이지 않는다
df$ell_centered <- df$ell - mean(df$ell)
m.ell.centered <- lm(api00 ~ ell_centered * yr_rnd, data = df)
summary(m.ell.centered)
coefs <- summary(m.ellxyr_rnd)$coefficients
coefs
api.br <- coefs[1]
api.nobr <- coefs[1]+coefs[3]
slope.red <- coefs[2]
slope.blue <- coefs[2]+coefs[4]
data.frame(api.br, api.nobr, slope.red, slope.blue)
# br과 nobr에 따라서 ell의 slope가 달라지는것이
# interaction effects
ro.04
> #######################################
> #######################################
> # with a continuous and a category IV
> #######################################
> #######################################
> m.mealsyr_rnd <- lm(api00~meals+yr_rnd, data=df)
> summary(m.mealsyr_rnd)
Call:
lm(formula = api00 ~ meals + yr_rnd, data = df)
Residuals:
Min 1Q Median 3Q Max
-174.882 -41.542 -1.044 39.454 176.007
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 885.6192 6.4712 136.855 < 2e-16 ***
meals -3.7921 0.1036 -36.595 < 2e-16 ***
yr_rndnobr -40.3291 7.8479 -5.139 4.35e-07 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 59.99 on 397 degrees of freedom
Multiple R-squared: 0.823, Adjusted R-squared: 0.8221
F-statistic: 923.2 on 2 and 397 DF, p-value: < 2.2e-16
>
> m.mealsxyr_rnd <- lm(api00~meals*yr_rnd, data=df)
> summary(m.mealsxyr_rnd)
Call:
lm(formula = api00 ~ meals * yr_rnd, data = df)
Residuals:
Min 1Q Median 3Q Max
-175.78 -41.96 -0.76 39.06 175.49
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 885.0046 6.7647 130.827 <2e-16 ***
meals -3.7805 0.1100 -34.354 <2e-16 ***
yr_rndnobr -31.8737 27.9104 -1.142 0.254
meals:yr_rndnobr -0.1041 0.3299 -0.316 0.752
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 60.06 on 396 degrees of freedom
Multiple R-squared: 0.8231, Adjusted R-squared: 0.8217
F-statistic: 614.1 on 3 and 396 DF, p-value: < 2.2e-16
>
> # 중요: 위에서 yr_rnd의 설명력이 t test를 보면
> # 사라짐. 왜?
> # interaction을 추가하면서 yr_rnd의 설명력은
> # (방학이 있고 없음의 설명력은) meals의 percentage가
> # 0 이 될때의 것으로 설정된다. 0이 되는 데이터가
> # 극단치이고 이 경우에 nobr인 학교도 없으므로 계산된
> # se값이 극대화되어 t값이 작아지게 된다. 이를 바로
> # 잡으려면 meals의 0이 중간에 가도록 한 후에
> # regression을 하면 meals가 0 일때가 극단치가 되지
> # 않으므로 결과가 옳게 나온다
>
> df %>%
+ ggplot(aes(x = meals, y = api00, color = factor(yr_rnd))) +
+ geom_point(size = 3) +
+ geom_smooth(method = "lm", se = FALSE) +
+ scale_color_brewer(palette = "Set1", name = "yr_rnd") +
+ labs(title = "Multiple Regression: Interaction by yr_rnd break",
+ x = " (free meals percentage)",
+ y = "api00")
`geom_smooth()` using formula = 'y ~ x'
>
>
> # 1. 평균을 빼준 값을 새로운 변인으로 저장
> df$meals_centered <-
+ df$meals - mean(df$meals)
>
> # 2. Run the model again using the centered variable
> m_centered <- lm(api00 ~ meals_centered * yr_rnd,
+ data = df)
>
> # 3. Check the new summary
> summary(m_centered)
Call:
lm(formula = api00 ~ meals_centered * yr_rnd, data = df)
Residuals:
Min 1Q Median 3Q Max
-175.78 -41.96 -0.76 39.06 175.49
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 656.9827 3.5150 186.910 < 2e-16 ***
meals_centered -3.7805 0.1100 -34.354 < 2e-16 ***
yr_rndnobr -38.1550 10.4473 -3.652 0.000295 ***
meals_centered:yr_rndnobr -0.1041 0.3299 -0.316 0.752386
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 60.06 on 396 degrees of freedom
Multiple R-squared: 0.8231, Adjusted R-squared: 0.8217
F-statistic: 614.1 on 3 and 396 DF, p-value: < 2.2e-16
>
> df %>%
+ ggplot(aes(x = meals_centered, y = api00, color = factor(yr_rnd))) +
+ geom_point(size = 3) +
+ geom_smooth(method = "lm", se = FALSE) +
+ scale_color_brewer(palette = "Set1", name = "yr_rnd") +
+ labs(title = "Multiple Regression: Interaction by yr_rnd break",
+ x = " (free meals percentage CENTERED)",
+ y = "api00")
`geom_smooth()` using formula = 'y ~ x'
>
>
> # Install the package if you do not have it
> # install.packages("interactions")
>
> # Plot the interaction
> library(interactions)
> interact_plot(m.mealsxyr_rnd, pred = meals, modx = yr_rnd)
>
>
>
> m.ellyr_rnd <- lm(api00~ell+yr_rnd, data=df)
> summary(m.ellyr_rnd)
Call:
lm(formula = api00 ~ ell + yr_rnd, data = df)
Residuals:
Min 1Q Median 3Q Max
-265.843 -55.349 0.334 70.673 195.778
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 784.3978 7.2878 107.63 < 2e-16 ***
ell -4.0427 0.2094 -19.31 < 2e-16 ***
yr_rndnobr -41.8420 12.3441 -3.39 0.00077 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 90.1 on 397 degrees of freedom
Multiple R-squared: 0.6008, Adjusted R-squared: 0.5988
F-statistic: 298.8 on 2 and 397 DF, p-value: < 2.2e-16
> m.ellxyr_rnd <- lm(api00~ell*yr_rnd, data=df)
> summary(m.ellxyr_rnd)
Call:
lm(formula = api00 ~ ell * yr_rnd, data = df)
Residuals:
Min 1Q Median 3Q Max
-270.514 -57.800 5.994 65.845 206.275
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 794.2576 7.8422 101.280 < 2e-16 ***
ell -4.4418 0.2420 -18.354 < 2e-16 ***
yr_rndnobr -110.5779 24.7935 -4.460 1.07e-05 ***
ell:yr_rndnobr 1.4884 0.4673 3.185 0.00156 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 89.08 on 396 degrees of freedom
Multiple R-squared: 0.6108, Adjusted R-squared: 0.6078
F-statistic: 207.1 on 3 and 396 DF, p-value: < 2.2e-16
> df %>%
+ ggplot(aes(x = ell, y = api00, color = factor(yr_rnd))) +
+ geom_point(size = 3) +
+ geom_smooth(method = "lm", se = FALSE) +
+ scale_color_brewer(palette = "Set1", name = "Cylinders") +
+ labs(title = "Multiple Regression: Interaction by yr_rnd break",
+ x = " (ell, english language learner)",
+ y = "api00")
`geom_smooth()` using formula = 'y ~ x'
>
>
> # ell mean centred 이번 경우는 큰 차이를 보이지 않는다
> df$ell_centered <- df$ell - mean(df$ell)
> m.ell.centered <- lm(api00 ~ ell_centered * yr_rnd, data = df)
> summary(m.ell.centered)
Call:
lm(formula = api00 ~ ell_centered * yr_rnd, data = df)
Residuals:
Min 1Q Median 3Q Max
-270.514 -57.800 5.994 65.845 206.275
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 654.5514 5.3323 122.752 < 2e-16 ***
ell_centered -4.4418 0.2420 -18.354 < 2e-16 ***
yr_rndnobr -63.7652 14.0117 -4.551 7.12e-06 ***
ell_centered:yr_rndnobr 1.4884 0.4673 3.185 0.00156 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 89.08 on 396 degrees of freedom
Multiple R-squared: 0.6108, Adjusted R-squared: 0.6078
F-statistic: 207.1 on 3 and 396 DF, p-value: < 2.2e-16
>
>
> coefs <- summary(m.ellxyr_rnd)$coefficients
> coefs
Estimate Std. Error t value Pr(>|t|)
(Intercept) 794.257647 7.8421894 101.280090 3.235236e-285
ell -4.441819 0.2420092 -18.353922 6.903303e-55
yr_rndnobr -110.577884 24.7935271 -4.459950 1.069371e-05
ell:yr_rndnobr 1.488362 0.4673174 3.184906 1.562618e-03
> api.br <- coefs[1]
> api.nobr <- coefs[1]+coefs[3]
> slope.red <- coefs[2]
> slope.blue <- coefs[2]+coefs[4]
> data.frame(api.br, api.nobr, slope.red, slope.blue)
api.br api.nobr slope.red slope.blue
1 794.2576 683.6798 -4.441819 -2.953456
> # br과 nobr에 따라서 ell의 slope가 달라지는것이
> # interaction effects
>
3 or more groups
만약에 ANOVA 테스트에서와 같이 종류가 3개 이상인 변인은 어떻게 처리해야 할까? 아래는 이를 regression으로 테스트 한 결과이다.
> mod2 <- lm(api00 ~ factor(mealcat), data=datavar)
> mod2
Call:
lm(formula = api00 ~ factor(mealcat), data = datavar)
Coefficients:
(Intercept) factor(mealcat)2 factor(mealcat)3
805.7 -166.3 -301.3
> summary(mod2)
Call:
lm(formula = api00 ~ factor(mealcat), data = datavar)
Residuals:
Min 1Q Median 3Q Max
-253.394 -47.883 0.282 52.282 185.620
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 805.718 6.169 130.60 <2e-16 ***
factor(mealcat)2 -166.324 8.708 -19.10 <2e-16 ***
factor(mealcat)3 -301.338 8.629 -34.92 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 70.61 on 397 degrees of freedom
Multiple R-squared: 0.7548, Adjusted R-squared: 0.7536
F-statistic: 611.1 on 2 and 397 DF, p-value: < 2.2e-16
>
아까와 같은 두 집단 간의 비교는 가능하였지만, 세 집단인 이 경우에 regression 테스트는 x를 연속변인으로 하는 선형관계를 구하게 된다. 즉, 두집단으로 이루어진 변인의 경우, 한 집단을 0으로 보았을 때 다른 집단은 자동으로 1이 되고, 이 둘을 비교하는 것이었는데, 3집단인 경우, 어느 한 집단을 0으로 놓고 본다고 비교할 집단이 두개나 되므로 비교가 어려워진다. 이는 현실에 맞지 않으므로 대개의 경우에는 집단 수에 해당하는 변인을 가외로 (가변인 혹은 dummy variable) 만든 후, 이를 가지고 regression을 하게 된다.
SPSS의 경우에는 아래와 같이 recode작업을 할 수 있다.
compute mealcat1 = 0. if mealcat = 1 mealcat1 = 1. compute mealcat2 = 0. if mealcat = 2 mealcat2 = 1. compute mealcat3 = 0. if mealcat = 3 mealcat3 = 1. execute.
위는 해당 카데고리를 1로 만들고, 나머지를 0으로 만들어서 2분화 하는 작업이다. 이렇게 하면 3개의 새로운 변인이 만들어지게 되는데 (mealcat1, mealcat2, mealcat3), 이 세개의 변인 중에서 2개만을 취해서 regression 테스트를 한다. SPSS의 경우에는
regression /dependent api00 /method = enter mealcat2 mealcat3. 주의. 세개의 변인을 모두 넣지 않는다.
Model Summary Model R R Square Adjusted R Square Std. Error of the Estimate 1 .869a .755 .754 70.612 a. Predictors: (Constant), mealcat3, mealcat2 ANOVA(b) Model Sum of Squares df Mean Square F Sig. 1 Regression 6094197.670 2 3047098.835 611.121 .000a Residual 1979474.328 397 4986.081 Total 8073671.997 399 a. Predictors: (Constant), mealcat3, mealcat2 b. Dependent Variable: api 2000 Coefficients(a) Unstandardized Coefficients Standardized Coefficients Model B Std. Error Beta t Sig. 1 (Constant) 805.718 6.169 130.599 .000 mealcat2 -166.324 8.708 -.550 -19.099 .000 mealcat3 -301.338 8.629 -1.007 -34.922 .000 a. Dependent Variable: api 2000
| Coefficients(a) | ||||||
| Unstandardized Coefficients | Standardized Coefficients | |||||
| Model | B | Std. Error | Beta | t | Sig. | |
| 1 | (Constant) | 805.718 | 6.169 | 130.599 | .000 | |
| mealcat2 | -166.324 | 8.708 | -.550 | -19.099 | .000 | |
| mealcat3 | -301.338 | 8.629 | -1.007 | -34.922 | .000 | |
| a. Dependent Variable: api 2000 | ||||||
위에서
$\hat{Y} = 805.718 - 166.324 \; \text{mealcat2} - 301.338 \; \text{mealcat3} $
이에 대한 해석은
- mealcat2와 mealcat3이 0일 때 (즉 mealcat1 변인의 상태일 때), 805.718
- mealcat3이 0일 때, 805.718-166.324 = 639.39 의 상황
- mealcat2가 0일 때, 805.718-301.338 = 504.38 의 상황이다.
- 그리고, 이 값들은 바로 각 그룹의 평균값이 된다.
Report api 2000 Percentage free meals in 3 categories Mean N Std. Deviation 0-46% free meals 805.72 131 65.669 47-80% free meals 639.39 132 82.135 81-100% free meals 504.38 137 62.727 Total 647.62 400 142.249
위에서, mealcat1대신에 mealcat3 그룹을 빼고 사용했어도, 결과를 해석하는데는 지장이 없다.
마지막으로, 위의 테스트는 이전에 언급되었던 FactorialAnova와 동일한 것이다.
glm api00 by mealcat /print=parameter.
| Between-Subjects Factors | |||
| Value Label | N | ||
| Percentage free meals in 3 categories | 1 | 0-46% free meals | 131 |
| 2 | 47-80% free meals | 132 | |
| 3 | 81-100% free meals | 137 | |
| Tests of Between-Subjects Effects | |||||
| Dependent Variable:api 2000 | |||||
| Source | Type III Sum of Squares | df | Mean Square | F | Sig. |
| Corrected Model | 6.094E6 | 2 | 3047098.835 | 611.121 | .000 |
| Intercept | 1.688E8 | 1 | 1.688E8 | 33863.695 | .000 |
| mealcat | 6094197.670 | 2 | 3047098.835 | 611.121 | .000 |
| Error | 1979474.328 | 397 | 4986.081 | ||
| Total | 1.758E8 | 400 | |||
| Corrected Total | 8073671.997 | 399 | |||
| a. R Squared = .755 (Adjusted R Squared = .754) | |||||
| Parameter Estimates | |||||||
| Dependent Variable:api 2000 | |||||||
| 95% Confidence Interval | |||||||
| Parameter | B | Std. Error | t | Sig. | Lower Bound | Upper Bound | |
| Intercept | 504.380 | 6.033 | 83.606 | .000 | 492.519 | 516.240 | |
| [mealcat=1] | 301.338 | 8.629 | 34.922 | .000 | 284.374 | 318.302 | |
| [mealcat=2] | 135.014 | 8.612 | 15.677 | .000 | 118.083 | 151.945 | |
| [mealcat=3] | 0a | . | . | . | . | . | |
| a. This parameter is set to zero because it is redundant. | |||||||
혹은 Oneway ANOVA
ONEWAY api00 BY mealcat /STATISTICS DESCRIPTIVES EFFECTS HOMOGENEITY /PLOT MEANS /MISSING ANALYSIS /POSTHOC=TUKEY SCHEFFE ALPHA(0.05).
| Sum of Squares | df | Mean Square | F | Sig. | |
| Between Groups | 6094197.67 | 2 | 3047098.835 | 611.120953 | .000 |
| Within Groups | 1979474.328 | 397 | 4986.08143 | ||
| Total | 8073671.998 | 399 |
2 variables, categorical
위에서 사용된 2 개의 독립변인을 모두 넣어서 regression을 할 수도 있다. 위에서 언급한 경로를 따른다면, 이는 FactorialAnova의 한 종류일 것이다.
regression /dep api00 /method = enter yr_rnd mealcat1 mealcat2.
| Model Summary | ||||
| Model | R | R Square | Adjusted R Square | Std. Error of the Estimate |
| 1 | .876a | .767 | .765 | 68.893 |
| a. Predictors: (Constant), mealcat2, year round school, mealcat1 | ||||
| ANOVA(b) | ||||||
| Model | Sum of Squares | df | Mean Square | F | Sig. | |
| 1 | Regression | 6194144.303 | 3 | 2064714.768 | 435.017 | .000a |
| Residual | 1879527.694 | 396 | 4746.282 | |||
| Total | 8073671.997 | 399 | ||||
| a. Predictors: (Constant), mealcat2, year round school, mealcat1 | ||||||
| b. Dependent Variable: api 2000 | ||||||
| Coefficients(a) | ||||||
|---|---|---|---|---|---|---|
| Unstandardized Coefficients | Standardized Coefficients | |||||
| Model | B | Std. Error | Beta | t | Sig. | |
| 1 | (Constant) | 526.330 | 7.585 | 69.395 | .000 | |
| year round school | -42.960 | 9.362 | -.127 | -4.589 | .000 | |
| mealcat1 | 281.683 | 9.446 | .930 | 29.821 | .000 | |
| mealcat2 | 117.946 | 9.189 | .390 | 12.836 | .000 | |
| a. Dependent Variable: api 2000 | ||||||
똑같은 분석이지만 뒤의 두 변인의 효과를 따로 보기 위해서 뽑은 결과이다.
regression /dep api00 /method = enter yr_rnd /method = test(mealcat1 mealcat2).
| Model Summary | ||||
|---|---|---|---|---|
| Model | R | R Square | Adjusted R Square | Std. Error of the Estimate |
| 1 | .475a | .226 | .224 | 125.300 |
| 2 | .876b | .767 | .765 | 68.893 |
| a. Predictors: (Constant), year round school b. Predictors: (Constant), year round school, mealcat2, mealcat1 |
||||
| ANOVA(d) | ||||||||
| Model | Sum of Squares | df | Mean Square | F | Sig. | R Square Change | ||
| 1 | Regression | 1825000.563 | 1 | 1825000.563 | 116.241 | .000a | ||
| Residual | 6248671.435 | 398 | 15700.179 | |||||
| Total | 8073671.997 | 399 | ||||||
| 2 | Subset Tests | mealcat1, mealcat2 | 4369143.740 | 2 | 2184571.870 | 460.270 | .000b | .541 |
| Regression | 6194144.303 | 3 | 2064714.768 | 435.017 | .000c | |||
| Residual | 1879527.694 | 396 | 4746.282 | |||||
| Total | 8073671.997 | 399 | ||||||
| a. Predictors: (Constant), year round school b. Tested against the full model. c. Predictors in the Full Model: (Constant), year round school, mealcat2, mealcat1. d. Dependent Variable: api 2000 |
||||||||
| Coefficients(a) | ||||||
|---|---|---|---|---|---|---|
| Unstandardized Coefficients | Standardized Coefficients | |||||
| Model | B | Std. Error | Beta | t | Sig. | |
| 1 | (Constant) | 684.539 | 7.140 | 95.878 | .000 | |
| year round school | -160.506 | 14.887 | -.475 | -10.782 | .000 | |
| 2 | (Constant) | 526.330 | 7.585 | 69.395 | .000 | |
| year round school | -42.960 | 9.362 | -.127 | -4.589 | .000 | |
| mealcat1 | 281.683 | 9.446 | .930 | 29.821 | .000 | |
| mealcat2 | 117.946 | 9.189 | .390 | 12.836 | .000 | |
| a. Dependent Variable: api 2000 | ||||||
| Excluded Variables(b) | ||||||
|---|---|---|---|---|---|---|
| Collinearity Statistics | ||||||
| Model | Beta In | t | Sig. | Partial Correlation | Tolerance | |
| 1 | mealcat1 | .697a | 23.132 | .000 | .758 | .914 |
| mealcat2 | -.138a | -3.106 | .002 | -.154 | .962 | |
| a. Predictors in the Model: (Constant), year round school b. Dependent Variable: api 2000 |
||||||
해석에 대해서 . . . .
| interpretation | |||
|---|---|---|---|
| mealcat=1 | mealcat=2 | mealcat=0 | |
| yr_rnd=0 | cell1 | cell2 | cell3 |
| yr_rnd=1 | cell4 | cell5 | cell6 |
| interpretation | |||
|---|---|---|---|
| mealcat=1 | mealcat=2 | mealcat=0 | |
| mealcat=1→1 | mealcat=2→1 | mealcat=3→mealcat1,2=0 | |
| yr_rnd=0 | cell1 | cell2 | cell3 |
| intercept + BMealCat1 | intercept + BMealCat2 | intercept | |
| yr_rnd=1 | cell4 | cell5 | cell6 |
| intercept + BMealCat1 + Byr_rnd | intercept + BMealCat2 + Byr_rnd | intercept + Byr_rnd |
|
glm api00 BY yr_rnd mealcat /DESIGN = yr_rnd mealcat /print=parameter TEST(LMATRIX).
continuous + categorical variables
regress /dep = api00 /method = enter yr_rnd some_col /save pre. * pre = predicted value (y hat). output: Model Summary(b) Model R R Square Adjusted R Square Std. Error of the Estimate 1 .507a .257 .253 122.951 a. Predictors: (Constant), parent some college, year round school b. Dependent Variable: api 2000 ANOVA(b) Model Sum of Squares df Mean Square F Sig. 1 Regression 2072201.839 2 1036100.919 68.539 .000a Residual 6001470.159 397 15117.053 Total 8073671.997 399 a. Predictors: (Constant), parent some college, year round school b. Dependent Variable: api 2000 Coefficients(a) Unstandardized Coefficients Standardized Coefficients Model B Std. Error Beta t Sig. 1 (Constant) 637.858 13.503 47.237 .000 year round school -149.159 14.875 -.442 -10.027 .000 parent some college 2.236 553 .178 4.044 .000 a. Dependent Variable: api 2000
COMPUTE filt=(yr_rnd=0). FILTER BY filt. regress /dep = api00 /method = enter some_col.
위의 명령어는 (spss) yr_rnd value → 0 인것을 선택하여, 이를 필터링하면 (고르면) → 1 이 되고
필터링되지 않은 케이스들은 버려지게 되어 필터링이 된 케이스들만 선택이 되어 분석에 사용됨을 뜻 한다. 즉, 위는 rn_rnd값이 0 인 케이스에 대해서만 simple regression을 하라는 것이다.
Model Summary Model R R Square Adjusted R Square Std. Error of the Estimate 1 .126a .016 .013 131.278 a. Predictors: (Constant), parent some college ANOVA(b) Model Sum of Squares df Mean Square F Sig. 1 Regression 84700.858 1 84700.858 4.915 .027a Residual 5273591.675 306 17233.960 Total 5358292.532 307 a. Predictors: (Constant), parent some college b. Dependent Variable: api 2000 Coefficients(a) Unstandardized Coefficients Standardized Coefficients Model B Std. Error Beta t Sig. 1 (Constant) 655.110 15.237 42.995 .000 parent some college 1.409 .636 .126 2.217 .027 a. Dependent Variable: api 2000
COMPUTE filt=(yr_rnd=1). FILTER BY filt. regress /dep = api00 /method = enter some_col. Model Summary Model R R Square Adjusted R Square Std. Error of the Estimate 1 .648a .420 .413 75.773 a. Predictors: (Constant), parent some college ANOVA(b) Model Sum of Squares df Mean Square F Sig. 1 Regression 373644.064 1 373644.064 65.078 .000a Residual 516734.838 90 5741.498 Total 890378.902 91 a. Predictors: (Constant), parent some college b. Dependent Variable: api 2000 Coefficients(a) Unstandardized Coefficients Standardized Coefficients Model B Std. Error Beta t Sig. 1 (Constant) 407.039 16.515 24.647 .000 parent some college 7.403 .918 .648 8.067 .000 a. Dependent Variable: api 2000
interaction effect
위의 두 regression은 yr_rnd 변인이 갖는 두 가지 특성에 대해서 따로 regression (api_00 ← some_col) 을 한 것이다. 이 결과, 두 집단의 regression 기울기 (coefficient)가 다르다는 것을 알았다. 즉, some_col의 api_00에 대한 영향력이 다르다는 것이다. 이는 각각의 상황(변인의 특성)에 따라서 동일한 독립변인이 역할을 달리하는 것으로 상호효과(interaction effect)라고 할 수 있다. 따라서, 두 기울기가 혹은 계수(coefficients)가 서로 다르다는 것을 검증한다면, 상호효과를 알아볼 수 있다.
아래는 새로운 변인을 만들어서 변인의 값으로 yr_rnd와 some_col값을 곱한 값을 대체한 것이다. 즉,
DV: api00
IV1: some_col
IV2: yr_rnd
IV3: yr_rnd * some_col = interaction effects
compute yrXsome = yr_rnd*some_col. execute.
그리고, 이 변인을 regression 공식에 이용한다.
regress /dep = api00 /method = enter some_col yr_rnd yrXsome /save pre.
output: Model Summary(b) Model R R Square Adjusted R Square Std. Error of the Estimate 1 .532a .283 .277 120.922 a. Predictors: (Constant), yrXsome, parent some college, year round school b. Dependent Variable: api 2000 ANOVA(b) Model Sum of Squares df Mean Square F Sig. 1 Regression 2283345.485 3 761115.162 52.053 .000a Residual 5790326.513 396 14622.037 Total 8073671.997 399 a. Predictors: (Constant), yrXsome, parent some college, year round school b. Dependent Variable: api 2000 Coefficients(a) Unstandardized Coefficients Standardized Coefficients Model B Std. Error Beta t Sig. 1 (Constant) 655.110 14.035 46.677 .000 parent some college 1.409 .586 .112 2.407 .017 year round school -248.071 29.859 -.735 -8.308 .000 __yrXsome__ 5.993 1.577 .330 3.800 .000 a. Dependent Variable: api 2000 Residuals Statistics(a) Minimum Maximum Mean Std. Deviation N Predicted Value 407.04 749.54 647.62 75.648 400 Residual -275.118 279.252 .000 120.466 400 Std. Predicted Value -3.180 1.347 .000 1.000 400 Std. Residual -2.275 2.309 .000 .996 400 a. Dependent Variable: api 2000
위에서 yr_rnd의 b 계수 값이 5.993으로 유의미하다고 판단된다 (t = 3.800, p < .001). 따라서 두 변인 간의 상호효과가 존재한다고 할 수 있다. 이를 다시 도표화해서 보면, 두 집단의 기울기가 서로 다르다는 것을 알 수 있다.
위의 테스트를 살펴보면, 두 개의 독립변인 중 하나는 종류변인이고 다른 하나는 숫자변인이다. 각 변인의 영향력에 대해서 regression을 통해서 알아보면서 두 변인의 상호작용까지 알아본 것이 된다. 이와 같은 절차는 FactorialAnova 에서 살펴본 것과 같다. 사실, 위의 연구문제(가설)를 ANOVA를 이용해서도 할 수 있다.
glm api00 BY yr_rnd WITH some_col /DESIGN = some_col yr_rnd yr_rnd*some_col. or UNIANOVA api00 BY yr_rnd WITH some_col /DESIGN=yr_rnd some_col some_col*yr_rnd /print=parameter .
using_dummy_variables.1781046721.txt.gz · Last modified: by hkimscil