User Tools

Site Tools


using_dummy_variables

This is an old revision of the document!


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는 이를 확인한 것.

regressioncategory.jpg

위의 그래프는 두 그룹의 (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")
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'
> 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

regression-cat0.jpg

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

regression-cat1.jpg

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-cat2-interctionx.jpg
regression-cat2-interaction-each.jpg

위의 테스트를 살펴보면, 두 개의 독립변인 중 하나는 종류변인이고 다른 하나는 숫자변인이다. 각 변인의 영향력에 대해서 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.1781046370.txt.gz · Last modified: by hkimscil

Donate Powered by PHP Valid HTML5 Valid CSS Driven by DokuWiki