update

Author: ysaito8015

ex_assign_and_eval/assign_and_eval.R

@@ -0,0 +1,85 @@
+# assign function
+a <- 10
+assign("a", 10) # equal
+# paste function
+x <- paste("Hello", "World", sep = "") # "HelloWorld"
+x <- paste("Hello", "World", sep = " ") # "Hello World"
+
+
+
+for(i in 1:6) { #-- Create objects  'r.1', 'r.2', ... 'r.6' --
+      nam <- paste("r", i, sep = ".")
+        assign(nam, 1:i) # r.1 -> 1, r.2 -> 1,2, ... r.6 -> 1,2,3,4,5,6
+}
+print(ls(pattern = "^r..$"))
+# [1] "r.1" "r.2" "r.3" "r.4" "r.5" "r.6"
+
+# create vecters
+# create objects by for loop
+
+for( i in seq(5)){
+  assign(paste("Test", i, sep = ""),
+         # sample func: if replace = T and sampling from over 200, Walket's alias method (Ripley, 1987) is used
+         sample(10:100, size = 10, replace = TRUE),
+         env = .GlobalEnv) # set as global variable
+}
+
+# create dataframes by for loop
+for( i in seq(5)) {
+  assign(paste("TestDF", i, sep = ""),
+         data.frame(TEST1 = sample(10:100, size = 10, replace = TRUE),
+                    TEST2 = sample(200:300, size = 10, replace = TRUE)),
+         env = .GlobalEnv)
+}
+
+print(ls())
+
+#  [1] "i"       "Test1"   "Test2"   "Test3"   "Test4"   "Test5"   "TestDF1"
+#  [8] "TestDF2" "TestDF3" "TestDF4" "TestDF5"
+
+# eval(parse(text = "commands"))
+# print Test1 to Test5 with labels
+for( i in seq(5)) {
+  cat(eval(parse(text = paste("Test", i, sep = ""))),
+      labels = paste("Test", i, ":", sep = ""), fill = TRUE)
+}
+
+# Test1: 86 10 47 73 17 59 42 12 77 99
+# Test2: 56 27 90 52 38 86 53 31 99 85
+# Test3: 68 37 88 50 11 56 41 30 58 60
+# Test4: 61 73 73 16 31 97 32 21 31 62
+# Test5: 31 92 56 16 94 62 92 13 16 50
+
+# Plot TestDF1 to TestDF5 in a canvas
+layout(matrix(seq(5), rep(1, 5)), 1, 5)
+ColSet <- c("#d9bb9c", "#4b61ba", "#deb7a0", "#a87963", "#28231e")
+
+for( i in seq(5)) {
+  plot(eval(parse(text = paste("TestDF", i, "[, 1]", sep = ""))),
+       eval(parse(text = paste("TestDF", i, "[, 2]", sep = ""))),
+       xlab = "", ylab = "", main = paste("TestDF", i, sep = ""),
+       col = ColSet[i], pch = 15, cex = 4)
+}
+
+# > paste("TestDF", 1, "[, 1]", sep = "")
+# [1] "TestDF1[, 1]"
+# > parse(text = paste("TestDF", 1, "[, 1]", sep = ""))
+# expression(TestDF1[, 1])
+# > eval(parse(text = paste("TestDF", 1, "[, 1]", sep = "")))
+#  [1] 78 62 95 68 95 36 87 55 18 29
+
+# eval(substitute(do.call()))
+
+# yarrr
+library(yarrr)
+
+par(mfrow = c(2,2))
+loop.vector <- 1:4
+
+for (i in loop.vector) {
+  x <- examscores[, i]
+  hist(x,
+       main = paste("Question", i),
+       xlab = "Scores",
+       xlim = c(0, 100))
+}

pseudoreplication_and_random_effect/barplotting_with_SE.R

@@ -0,0 +1,24 @@
+rm(list = ls())
+yieldData <- read.table("./yieldData.dat", header = TRUE)
+
+library(lme4)
+lmer.m <- lmer(Yield ~ Treatment + (1|Plot), yieldData)
+
+library(emmeans)
+lsm.lmer.m <- lsmeans(lmer.m, specs="Treatment")
+
+summ.lsm.lmer.m <- summary(lsm.lmer.m)
+lsm <- summ.lsm.lmer.m$lsmean
+ui <- summ.lsm.lmer.m$lsmean + summ.lsm.lmer.m$SE
+li <- summ.lsm.lmer.m$lsmean - summ.lsm.lmer.m$SE
+
+par(las = 1,
+    cex = 1.5,
+    family = "sans",
+    lwd = 2)
+x <- barplot(lsm,
+             ylim = c(0, 80),
+             names.arg = c("Control", "Pesticide"),
+             xlab = "Treatment",
+             ylab = "Yield/plant (mg)")
+arrows(x, ui, x, li, angle = 90, code = 3)

pseudoreplication_and_random_effect/least-square-means.R

@@ -4,7 +4,7 @@ yieldData <- read.table("./yieldData.dat", header = TRUE)
 library(lme4)
 lmer.m <- lmer(Yield ~ Treatment + (1|Plot), yieldData)
 
-library(lsmeans)
+library(emmeans)
 print(lsm.lmer.m <- lsmeans(lmer.m, specs="Treatment"))
 
 # Loading required package: emmeans

repeated_measures_ANOVA_with_measures_over_time/least-square-means.R

@@ -0,0 +1,107 @@
+rm(list = ls())
+growthData <- read.table("./growthData.dat", header = TRUE)
+growthData$Day <- as.factor(growthData$Day)
+attach(growthData)
+
+library(lme4)
+library(emmeans)
+
+growth.lmer.m <- lmer(log(Mass) ~ Food * Day + (1|ID), growthData)
+
+# 同じ日の中で, 処理間で体重が異なるか
+print(lsmeans(growth.lmer.m, specs = c("Food", "Day")))
+
+# Loading required package: Matrix
+#  Food Day lsmean     SE   df lower.CL upper.CL
+#  A    1     3.89 0.0367 12.2     3.81     3.97
+#  B    1     3.93 0.0367 12.2     3.85     4.01
+#  A    3     4.06 0.0420 13.9     3.97     4.15
+#  B    3     4.34 0.0367 12.2     4.26     4.42
+#  A    5     4.24 0.0420 13.9     4.15     4.33
+#  B    5     4.74 0.0420 13.9     4.65     4.83
+# 
+# Degrees-of-freedom method: kenward-roger 
+# Results are given on the log (not the response) scale. 
+# Confidence level used: 0.95 
+
+cat("\n")
+# すべての組み合わせで事後検定を行う
+print(pairs(lsmeans(growth.lmer.m, specs = c("Food", "Day"))))
+
+#  contrast  estimate     SE    df t.ratio p.value
+#  A 1 - B 1   -0.036 0.0520 12.16  -0.693  0.9794
+#  A 1 - A 3   -0.168 0.0456  9.61  -3.680  0.0384
+#  A 1 - B 3   -0.450 0.0520 12.16  -8.658  <.0001
+#  A 1 - A 5   -0.348 0.0456  9.61  -7.628  0.0002
+#  A 1 - B 5   -0.852 0.0558 13.19 -15.282  <.0001
+#  B 1 - A 3   -0.132 0.0558 13.17  -2.360  0.2380
+#  B 1 - B 3   -0.414 0.0408  9.05 -10.153  <.0001
+#  B 1 - A 5   -0.312 0.0558 13.17  -5.585  0.0009
+#  B 1 - B 5   -0.816 0.0455  9.57 -17.930  <.0001
+#  A 3 - B 3   -0.282 0.0558 13.17  -5.054  0.0023
+#  A 3 - A 5   -0.180 0.0511 10.47  -3.522  0.0442
+#  A 3 - B 5   -0.685 0.0594 13.88 -11.527  <.0001
+#  B 3 - A 5    0.102 0.0558 13.17   1.829  0.4814
+#  B 3 - B 5   -0.402 0.0455  9.57  -8.840  0.0001
+#  A 5 - B 5   -0.505 0.0594 13.88  -8.496  <.0001
+# 
+# Degrees-of-freedom method: kenward-roger 
+# Results are given on the log (not the response) scale. 
+# P value adjustment: tukey method for comparing a family of 6 estimates 
+
+cat("\n")
+# Day で場合分けし, Food に注目する
+print(lsmeans(growth.lmer.m, specs = "Food", by = "Day"))
+
+# Day = 1:
+#  Food lsmean     SE   df lower.CL upper.CL
+#  A      3.89 0.0367 12.2     3.81     3.97
+#  B      3.93 0.0367 12.2     3.85     4.01
+# 
+# Day = 3:
+#  Food lsmean     SE   df lower.CL upper.CL
+#  A      4.06 0.0420 13.9     3.97     4.15
+#  B      4.34 0.0367 12.2     4.26     4.42
+# 
+# Day = 5:
+#  Food lsmean     SE   df lower.CL upper.CL
+#  A      4.24 0.0420 13.9     4.15     4.33
+#  B      4.74 0.0420 13.9     4.65     4.83
+# 
+# Degrees-of-freedom method: kenward-roger 
+# Results are given on the log (not the response) scale. 
+# Confidence level used: 0.95 
+
+cat("\n")
+# pairs 関数で囲むと, Day で場合分けした水準ごとに(水準内), 繰り返し事後検定を実行する
+print(pairs(lsmeans(growth.lmer.m, specs = "Food", by = "Day"), adjust = "tukey"))
+
+# Day = 1:
+#  contrast estimate     SE   df t.ratio p.value
+#  A - B      -0.036 0.0520 12.2  -0.693  0.5014
+# 
+# Day = 3:
+#  contrast estimate     SE   df t.ratio p.value
+#  A - B      -0.282 0.0558 13.2  -5.054  0.0002
+# 
+# Day = 5:
+#  contrast estimate     SE   df t.ratio p.value
+#  A - B      -0.505 0.0594 13.9  -8.496  <.0001
+# 
+# Degrees-of-freedom method: kenward-roger 
+# Results are given on the log (not the response) scale. 
+
+cat("\n")
+# 3 つの比較全体, 水準間の p 値の調整を実行する
+print(rbind(pairs(lsmeans(growth.lmer.m, specs = "Food", by = "Day")), adjust = "tukey"))
+
+# Note: adjust = "tukey" was changed to "sidak"
+# because "tukey" is only appropriate for one set of pairwise comparisons
+#  Day contrast estimate     SE   df t.ratio p.value
+#  1   A - B      -0.036 0.0520 12.2  -0.693  0.8760
+#  3   A - B      -0.282 0.0558 13.2  -5.054  0.0006
+#  5   A - B      -0.505 0.0594 13.9  -8.496  <.0001
+# 
+# Degrees-of-freedom method: kenward-roger 
+# Results are given on the log (not the response) scale. 
+# P value adjustment: sidak method for 3 tests 

repeated_measures_ANOVA_with_measures_over_time/repeated-measures-ANOVA.R

@@ -0,0 +1,25 @@
+rm(list = ls())
+growthData <- read.table("./growthData.dat", header = TRUE)
+growthData$Day <- as.factor(growthData$Day)
+attach(growthData)
+
+library(lme4)
+library(car)
+
+growth.lmer.m <- lmer(log(Mass) ~ Food * Day + (1|ID), growthData)
+print(Anova(growth.lmer.m, test = "F"))
+
+# The following objects are masked from growthData (pos = 3):
+# 
+#     Day, Food, ID, Mass
+# 
+# Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
+# 
+# Response: log(Mass)
+#                F Df Df.res    Pr(>F)    
+# Food      35.249  1 5.7562 0.0011830 **   # 個体差を無視した, two-way ANOVA の分析値より p 値が下がっている
+# Day      165.941  2 9.5603 3.777e-08 ***
+# Food:Day  27.124  2 9.6315 0.0001104 ***
+# Df.res, 分母自由度が整数ではなくなっている
+# ---
+# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1