User Tools

Site Tools


c:ma:anova_note

Differences

This shows you the differences between two versions of the page.

Link to this comparison view

Both sides previous revisionPrevious revision
Next revision
Previous revision
c:ma:anova_note [2024/09/19 07:39] hkimscilc:ma:anova_note [2024/09/23 10:54] (current) – [ANOVA e.g.1] hkimscil
Line 4: Line 4:
 <code> <code>
 set.seed(10) set.seed(10)
-sa <- rnorm2(100, 11010+sa <- rnorm2(100, 7.62
-sb <- rnorm2(100, 11310)+sb <- rnorm2(100, 8.32)
  
-sa.mean <- mean(sa) +mean.sa <- mean(sa) 
-sb.mean <- mean(sb)+mean.sb <- mean(sb)
  
 na <- length(sa) na <- length(sa)
Line 31: Line 31:
  
 se <- sqrt((vp/na)+(vp/nb)) se <- sqrt((vp/na)+(vp/nb))
-diff <- sa.mean - sb.mean +diff <- mean.sa - mean.sb
 t.cal <- diff/se t.cal <- diff/se
 se se
Line 49: Line 49:
 sall <- data.frame(sab,group) sall <- data.frame(sab,group)
 sall sall
-attach(sall)+ 
 +attach(sall)
  
 df.tot <- (na+nb)-1 df.tot <- (na+nb)-1
-ss.tot <- var(sall$sab)*df.tot +ss.tot <- var(sab)*df.tot 
-var.tot <- var(sall$sab)+var.tot <- var(sab)
 df.tot df.tot
 ss.tot ss.tot
 var.tot var.tot
 var.tot*df.tot var.tot*df.tot
 +
 +# for uniqueN function data.table required
 +# install.packages("data.table"
 +# library(data.table) 
  
 uniqueN(sall, by = c("group")) uniqueN(sall, by = c("group"))
Line 76: Line 81:
  
 mean.grand <- mean(sab) mean.grand <- mean(sab)
-mean.sa <- mean(sa) +mean.sa <- mean(sa) 
-mean.sb <- mean(sb)+mean.sb <- mean(sb)
 mean.grand mean.grand
 mean.sa mean.sa
Line 115: Line 120:
  
 t.test(sa,sb) t.test(sa,sb)
-str(t.test(sa,sb))+str(t.test(sa,sb))
 t.test(sa,sb)$statistic t.test(sa,sb)$statistic
 t.test(sa,sb)$statistic^2 t.test(sa,sb)$statistic^2
Line 126: Line 131:
 ====== ANOVA e.g.1 ====== ====== ANOVA e.g.1 ======
 <code r> <code r>
-# ANOVA 
 set.seed(1024) set.seed(1024)
-<-30 +na <-50 
-s1 <- round(rnorm(n5.6, 1),0) +nb <- 50 
-s2 <- round(rnorm(n5.9, 1),0) +nc <- 50 
-s3 <- round(rnorm(n6.51),0)+ 
 +s1 <- round(rnorm(na11.6, 2),0) 
 +s2 <- round(rnorm(nb12.9, 2),0) 
 +s3 <- round(rnorm(nc13.42),0)
 s1 s1
 s2 s2
 s3 s3
  
-sam <- c(s1,s2,s3) +sabc <- c(s1,s2,s3) 
-sam +sabc 
-group <- c(rep("g1",n), rep("g2",n), rep("g3",n))+group <- c(rep("g1",na), rep("g2",nb), rep("g3",nc))
 group group
  
-sa <- data.frame(sam,group) +sall <- data.frame(sabc, group) 
-sa +sall 
-attach(sa+attach(sall
-aov.sa <- aov(sam ~ group, data=sa+aov.sall <- aov(sabc ~ group, data=sall
-summary(aov.sa)+summary(aov.sall)
  
-df.tot <- (3*n)-1 +df.total <- length(sabc) - 1 
-ss.tot <- var(sam)*df.tot +ss.total <- var(sabc)*df.total 
-var.tot <- var(sam+var.total <- var(sabc
-df.tot +df.total 
-ss.tot +ss.total 
-var.tot+var.total 
 + 
 +df.s1 <- na-1 
 +df.s2 <- nb-1 
 +df.s3 <- nc-1
  
-df.bet <- 3 - 1 
-df.s1 <- n-1 
-df.s2 <- n-1 
-df.s3 <- n-1 
 df.within <- df.s1+df.s2+df.s3 df.within <- df.s1+df.s2+df.s3
 df.within df.within
Line 170: Line 177:
 ss.within ss.within
  
-mean.grand <- mean(sam)+ 
 +# for uniqueN function data.table required 
 +# install.packages("data.table")  
 +library(data.table)  
 +uniqueN(sall, by = c("group")) 
 +nofg <- uniqueN(sall, by = c("group")) 
 +nofg 
 +df.between <- nofg - 1 
 +df.between  
 + 
 +mean.grand <- mean(sabc)
 mean.s1 <- mean(s1) mean.s1 <- mean(s1)
 mean.s2 <- mean(s2) mean.s2 <- mean(s2)
Line 179: Line 196:
 mean.s3 mean.s3
  
-diff.g1 <- * (mean.grand - mean.s1)^2 +error.g1 <- na * (mean.grand - mean.s1)^2 
-diff.g2 <- * (mean.grand - mean.s2)^2 +error.g2 <- nb * (mean.grand - mean.s2)^2 
-diff.g3 <- * (mean.grand - mean.s3)^2+error.g3 <- nc * (mean.grand - mean.s3)^2
  
-diff.g1  +error.g1  
-diff.g2 +error.g2 
-diff.g3+error.g3
  
-ss.bet <-  diff.g1 + diff.g2 + diff.g3 +ss.between <-  error.g1 + error.g2 + error.g3 
-ss.bet+ss.between
  
 # sumup # sumup
-ss.tot +ss.total 
-ss.bet+ss.between
 ss.within ss.within
-ss.bet+ss.within+ss.between + ss.within 
 +ss.total
  
-df.tot  +df.total 
-df.bet+df.between
 df.within df.within
-df.bet + df.within+df.between + df.within
  
-ms.bet <- ss.bet / df.bet +ms.between <- ss.between / df.between 
-ms.wit <- ss.within / df.within +ms.within <- ss.within / df.within 
-ms.bet +ms.between 
-ms.wit +ms.within
- +
-fvalue <- ms.bet/ms.wit+
  
 +fvalue <- ms.between/ms.within
 fvalue  fvalue 
  
-f.res <- (aov(sa$sam~sa$group))+# fvalue에서 판단했을 때 그 판단이 잘못일 확률 
 +1 - pf(fvalue, df.between, df.within) 
 + 
 + 
 +f.res <- aov(sabc ~ group, data=sall)
 summary(f.res) summary(f.res)
 +
 +# for regression 
 +r.res <- lm(sabc ~ group, data=sall)
 +summary(r.res)
 +anova(r.res)
 +summary(r.res)$r.square
 +
 +ss.total
 +ss.between
 +ss.within 
 +
 +# this is r.square value
 +ss.between/ss.total
 +</code>
 +
 +<code>
 +> set.seed(1024)
 +> na <-50
 +> nb <- 50
 +> nc <- 50
 +
 +> s1 <- round(rnorm(na, 11.6, 2),0)
 +> s2 <- round(rnorm(nb, 12.9, 2),0)
 +> s3 <- round(rnorm(nc, 13.4, 2),0)
 +> s1
 + [1] 10 11  8 10 12  7 11 16 14 13  8  9 14 12 11 15  9 11 10  9 11 13 15 14 11 11 12 13 10 13 10 12
 +[33] 14 11 14 11 14 11 10 11 11 14 10  9 14 13 10 10  7 13
 +> s2
 + [1] 12 12 16 14 12 12 12 12 11 12 13 12 10 13 15 12 12 18 14 13 13 13 17 13 13 14 11 13 13 11 14 15
 +[33] 11 12 10 14 13 12 14 15 13 10 10 17 12 14 14 16 13 12
 +> s3
 + [1] 13 12 12 16 14 14  9 12 11 14 15 13  7 17 16 12 12 12  9 11 12 16 17 13 18 14 14 12 15 11 13 12
 +[33] 13 10 14 10 15 10 11 14 10 14 14 13 13 13 11 11 12 11
 +
 +> sabc <- c(s1,s2,s3)
 +> sabc
 +  [1] 10 11  8 10 12  7 11 16 14 13  8  9 14 12 11 15  9 11 10  9 11 13 15 14 11 11 12 13 10 13 10 12
 + [33] 14 11 14 11 14 11 10 11 11 14 10  9 14 13 10 10  7 13 12 12 16 14 12 12 12 12 11 12 13 12 10 13
 + [65] 15 12 12 18 14 13 13 13 17 13 13 14 11 13 13 11 14 15 11 12 10 14 13 12 14 15 13 10 10 17 12 14
 + [97] 14 16 13 12 13 12 12 16 14 14  9 12 11 14 15 13  7 17 16 12 12 12  9 11 12 16 17 13 18 14 14 12
 +[129] 15 11 13 12 13 10 14 10 15 10 11 14 10 14 14 13 13 13 11 11 12 11
 +> group <- c(rep("g1",na), rep("g2",nb), rep("g3",nc))
 +> group
 +  [1] "g1" "g1" "g1" "g1" "g1" "g1" "g1" "g1" "g1" "g1" "g1" "g1" "g1" "g1" "g1" "g1" "g1" "g1" "g1"
 + [20] "g1" "g1" "g1" "g1" "g1" "g1" "g1" "g1" "g1" "g1" "g1" "g1" "g1" "g1" "g1" "g1" "g1" "g1" "g1"
 + [39] "g1" "g1" "g1" "g1" "g1" "g1" "g1" "g1" "g1" "g1" "g1" "g1" "g2" "g2" "g2" "g2" "g2" "g2" "g2"
 + [58] "g2" "g2" "g2" "g2" "g2" "g2" "g2" "g2" "g2" "g2" "g2" "g2" "g2" "g2" "g2" "g2" "g2" "g2" "g2"
 + [77] "g2" "g2" "g2" "g2" "g2" "g2" "g2" "g2" "g2" "g2" "g2" "g2" "g2" "g2" "g2" "g2" "g2" "g2" "g2"
 + [96] "g2" "g2" "g2" "g2" "g2" "g3" "g3" "g3" "g3" "g3" "g3" "g3" "g3" "g3" "g3" "g3" "g3" "g3" "g3"
 +[115] "g3" "g3" "g3" "g3" "g3" "g3" "g3" "g3" "g3" "g3" "g3" "g3" "g3" "g3" "g3" "g3" "g3" "g3" "g3"
 +[134] "g3" "g3" "g3" "g3" "g3" "g3" "g3" "g3" "g3" "g3" "g3" "g3" "g3" "g3" "g3" "g3" "g3"
 +
 +> sall <- data.frame(sabc, group)
 +> sall
 +    sabc group
 +1     10    g1
 +2     11    g1
 +3      8    g1
 +4     10    g1
 +5     12    g1
 +6      7    g1
 +7     11    g1
 +8     16    g1
 +9     14    g1
 +10    13    g1
 +11        g1
 +12        g1
 +13    14    g1
 +14    12    g1
 +15    11    g1
 +16    15    g1
 +17        g1
 +18    11    g1
 +19    10    g1
 +20        g1
 +21    11    g1
 +22    13    g1
 +23    15    g1
 +24    14    g1
 +25    11    g1
 +26    11    g1
 +27    12    g1
 +28    13    g1
 +29    10    g1
 +30    13    g1
 +31    10    g1
 +32    12    g1
 +33    14    g1
 +34    11    g1
 +35    14    g1
 +36    11    g1
 +37    14    g1
 +38    11    g1
 +39    10    g1
 +40    11    g1
 +41    11    g1
 +42    14    g1
 +43    10    g1
 +44        g1
 +45    14    g1
 +46    13    g1
 +47    10    g1
 +48    10    g1
 +49        g1
 +50    13    g1
 +51    12    g2
 +52    12    g2
 +53    16    g2
 +54    14    g2
 +55    12    g2
 +56    12    g2
 +57    12    g2
 +58    12    g2
 +59    11    g2
 +60    12    g2
 +61    13    g2
 +62    12    g2
 +63    10    g2
 +64    13    g2
 +65    15    g2
 +66    12    g2
 +67    12    g2
 +68    18    g2
 +69    14    g2
 +70    13    g2
 +71    13    g2
 +72    13    g2
 +73    17    g2
 +74    13    g2
 +75    13    g2
 +76    14    g2
 +77    11    g2
 +78    13    g2
 +79    13    g2
 +80    11    g2
 +81    14    g2
 +82    15    g2
 +83    11    g2
 +84    12    g2
 +85    10    g2
 +86    14    g2
 +87    13    g2
 +88    12    g2
 +89    14    g2
 +90    15    g2
 +91    13    g2
 +92    10    g2
 +93    10    g2
 +94    17    g2
 +95    12    g2
 +96    14    g2
 +97    14    g2
 +98    16    g2
 +99    13    g2
 +100   12    g2
 +101   13    g3
 +102   12    g3
 +103   12    g3
 +104   16    g3
 +105   14    g3
 +106   14    g3
 +107    9    g3
 +108   12    g3
 +109   11    g3
 +110   14    g3
 +111   15    g3
 +112   13    g3
 +113    7    g3
 +114   17    g3
 +115   16    g3
 +116   12    g3
 +117   12    g3
 +118   12    g3
 +119    9    g3
 +120   11    g3
 +121   12    g3
 +122   16    g3
 +123   17    g3
 +124   13    g3
 +125   18    g3
 +126   14    g3
 +127   14    g3
 +128   12    g3
 +129   15    g3
 +130   11    g3
 +131   13    g3
 +132   12    g3
 +133   13    g3
 +134   10    g3
 +135   14    g3
 +136   10    g3
 +137   15    g3
 +138   10    g3
 +139   11    g3
 +140   14    g3
 +141   10    g3
 +142   14    g3
 +143   14    g3
 +144   13    g3
 +145   13    g3
 +146   13    g3
 +147   11    g3
 +148   11    g3
 +149   12    g3
 +150   11    g3
 +> attach(sall)
 +The following objects are masked _by_ .GlobalEnv:
 +
 +    group, sabc
 +
 +> aov.sall <- aov(sabc ~ group, data=sall)
 +> summary(aov.sall)
 +             Df Sum Sq Mean Sq F value   Pr(>F)    
 +group           68.7   34.33   8.049 0.000482 ***
 +Residuals   147  626.9    4.26                     
 +---
 +Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
 +
 +> df.total <- length(sabc) - 1
 +> ss.total <- var(sabc)*df.total
 +> var.total <- var(sabc)
 +> df.total
 +[1] 149
 +> ss.total
 +[1] 695.5733
 +> var.total
 +[1] 4.668277
 +
 +> df.s1 <- na-1
 +> df.s2 <- nb-1
 +> df.s3 <- nc-1
 +
 +> df.within <- df.s1+df.s2+df.s3
 +> df.within
 +[1] 147
 +
 +> ss.s1 <- var(s1)*df.s1
 +> ss.s2 <- var(s2)*df.s2
 +> ss.s3 <- var(s3)*df.s3
 +> ss.s1
 +[1] 222.32
 +> ss.s2
 +[1] 160.98
 +> ss.s3
 +[1] 243.62
 +> ss.within <- ss.s1+ss.s2+ss.s3
 +> ss.within
 +[1] 626.92
 +
 +
 +> # for uniqueN function data.table required
 +> # install.packages("data.table"
 +> library(data.table) 
 +data.table 1.14.8 using 8 threads (see ?getDTthreads).  Latest news: r-datatable.com
 +> uniqueN(sall, by = c("group"))
 +[1] 3
 +> nofg <- uniqueN(sall, by = c("group"))
 +> nofg
 +[1] 3
 +> df.between <- nofg - 1
 +> df.between 
 +[1] 2
 +
 +> mean.grand <- mean(sabc)
 +> mean.s1 <- mean(s1)
 +> mean.s2 <- mean(s2)
 +> mean.s3 <- mean(s3)
 +> mean.grand
 +[1] 12.38667
 +> mean.s1
 +[1] 11.44
 +> mean.s2
 +[1] 12.98
 +> mean.s3
 +[1] 12.74
 +
 +> error.g1 <- na * (mean.grand - mean.s1)^2
 +> error.g2 <- nb * (mean.grand - mean.s2)^2
 +> error.g3 <- nc * (mean.grand - mean.s3)^2
 +
 +> error.g1 
 +[1] 44.80889
 +> error.g2
 +[1] 17.60222
 +> error.g3
 +[1] 6.242222
 +
 +> ss.between <-  error.g1 + error.g2 + error.g3
 +> ss.between
 +[1] 68.65333
 +
 +> # sumup
 +> ss.total
 +[1] 695.5733
 +> ss.between
 +[1] 68.65333
 +> ss.within
 +[1] 626.92
 +> ss.between + ss.within
 +[1] 695.5733
 +> ss.total
 +[1] 695.5733
 +
 +> df.total
 +[1] 149
 +> df.between
 +[1] 2
 +> df.within
 +[1] 147
 +> df.between + df.within
 +[1] 149
 +
 +> ms.between <- ss.between / df.between
 +> ms.within <- ss.within / df.within
 +> ms.between
 +[1] 34.32667
 +> ms.within
 +[1] 4.264762
 +
 +> fvalue <- ms.between/ms.within
 +> fvalue 
 +[1] 8.048906
 +
 +> # fvalue에서 판단했을 때 그 판단이 잘못일 확률
 +> 1 - pf(fvalue, df.between, df.within)
 +[1] 0.0004818216
 +
 +
 +> f.res <- aov(sabc ~ group, data=sall)
 +> summary(f.res)
 +             Df Sum Sq Mean Sq F value   Pr(>F)    
 +group           68.7   34.33   8.049 0.000482 ***
 +Residuals   147  626.9    4.26                     
 +---
 +Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
 +
 +> # for regression 
 +> r.res <- lm(sabc ~ group, data=sall)
 +> summary(r.res)
 +
 +Call:
 +lm(formula = sabc ~ group, data = sall)
 +
 +Residuals:
 +   Min     1Q Median     3Q    Max 
 + -5.74  -1.44  -0.21   1.26   5.26 
 +
 +Coefficients:
 +            Estimate Std. Error t value Pr(>|t|)    
 +(Intercept)  11.4400     0.2921  39.171  < 2e-16 ***
 +groupg2       1.5400     0.4130   3.729 0.000274 ***
 +groupg3       1.3000     0.4130   3.148 0.001994 ** 
 +---
 +Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
 +
 +Residual standard error: 2.065 on 147 degrees of freedom
 +Multiple R-squared:  0.0987, Adjusted R-squared:  0.08644 
 +F-statistic: 8.049 on 2 and 147 DF,  p-value: 0.0004818
 +
 +> anova(r.res)
 +Analysis of Variance Table
 +
 +Response: sabc
 +           Df Sum Sq Mean Sq F value    Pr(>F)    
 +group        68.65  34.327  8.0489 0.0004818 ***
 +Residuals 147 626.92   4.265                      
 +---
 +Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
 +> summary(r.res)$r.square
 +[1] 0.09870035
 +
 +> ss.total
 +[1] 695.5733
 +> ss.between
 +[1] 68.65333
 +> ss.within 
 +[1] 626.92
 +
 +> # this is r.square value
 +> ss.between/ss.total
 +[1] 0.09870035
 +
 </code> </code>
 ====== ANOVA e.g. 2 ====== ====== ANOVA e.g. 2 ======
c/ma/anova_note.1726699171.txt.gz · Last modified: 2024/09/19 07:39 by hkimscil

Donate Powered by PHP Valid HTML5 Valid CSS Driven by DokuWiki