From 5f77530f8fd2f3f9db2215c61baab65c479bce5b Mon Sep 17 00:00:00 2001 From: "Jason K. Moore" Date: Mon, 10 Aug 2015 14:34:02 -0700 Subject: [PATCH 1/2] Started a mixed effects example for one variable RCBD. --- examples/rcbd_one_var_year.R | 41 +++++++++++++++++++++++++ examples/wheat_yield_data_with_year.csv | 33 ++++++++++++++++++++ 2 files changed, 74 insertions(+) create mode 100644 examples/rcbd_one_var_year.R create mode 100644 examples/wheat_yield_data_with_year.csv diff --git a/examples/rcbd_one_var_year.R b/examples/rcbd_one_var_year.R new file mode 100644 index 0000000..8f6da6a --- /dev/null +++ b/examples/rcbd_one_var_year.R @@ -0,0 +1,41 @@ +# This is an example of a single variate RCBD that has an addition random +# effect. +# +# To run the script type: +# source('rcbd_one_var_year.R', print.eval = TRUE) +# +# The code that will be displayed to the user is bounded by this separator: +#-----------------------------------------------------------------------------# + +# Load the necessary libraries. +#-----------------------------------------------------------------------------# +library('nlme') # for lme() +#-----------------------------------------------------------------------------# + +sep <- function(n){ + # This function simply prints a line of "=" to the screen as a separator and + # will not show up in the app. It is only here for nice printing. + line <- paste0(paste(replicate(n, "="), collapse = ""), '\n') + cat(line) +} + +# Load the data. +#-----------------------------------------------------------------------------# +my.data <- read.csv("wheat_yield_data_with_year.csv") +#-----------------------------------------------------------------------------# + +# Set Block, Treatment, and Year as factors. +#-----------------------------------------------------------------------------# +my.data$Block <- as.factor(my.data$Block) +my.data$Treatment <- as.factor(my.data$Treatment) +my.data$Year <- as.factor(my.data$Year) +#-----------------------------------------------------------------------------# + +# Construct a mized effects model and print the summary results. +cat('Linear mixed effect model results\n') +sep(79) +#-----------------------------------------------------------------------------# +model <- lme(fixed = Yield ~ Treatment, random = ~1|Block/Year, data = my.data) +summary(model) +#-----------------------------------------------------------------------------# +sep(79) diff --git a/examples/wheat_yield_data_with_year.csv b/examples/wheat_yield_data_with_year.csv new file mode 100644 index 0000000..cc23036 --- /dev/null +++ b/examples/wheat_yield_data_with_year.csv @@ -0,0 +1,33 @@ +Year,Block,Treatment,Yield +1999,1,Variety A,4.776113643 +1999,2,Variety A,4.167309262 +1999,3,Variety A,7.449164952 +1999,4,Variety A,7.757967337 +1999,1,Variety B,8.6989031 +1999,2,Variety B,6.332884039 +1999,3,Variety B,9.673666482 +1999,4,Variety B,10.81933243 +1999,1,Variety C,9.663780964 +1999,2,Variety C,7.861877074 +1999,3,Variety C,12.72423637 +1999,4,Variety C,8.837111214 +1999,1,Variety D,5.968923961 +1999,2,Variety D,7.078335949 +1999,3,Variety D,9.021241207 +1999,4,Variety D,8.428939255 +2000,1,Variety A,3.91823493 +2000,2,Variety A,4.58761123 +2000,3,Variety A,7.64636576 +2000,4,Variety A,7.91800989 +2000,1,Variety B,9.06640113 +2000,2,Variety B,6.31342374 +2000,3,Variety B,7.57774909 +2000,4,Variety B,10.95260769 +2000,1,Variety C,9.41454649 +2000,2,Variety C,8.40584383 +2000,3,Variety C,13.03441029 +2000,4,Variety C,7.72061644 +2000,1,Variety D,6.65400835 +2000,2,Variety D,6.56413817 +2000,3,Variety D,8.33194168 +2000,4,Variety D,8.85568545 From d895b82dcda11b658d8d682aad2df50468bbf7ef Mon Sep 17 00:00:00 2001 From: msimmond Date: Tue, 11 Aug 2015 10:51:58 -0700 Subject: [PATCH 2/2] block is nested within year (Year/Block) --- examples/rcbd_one_var_year.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/rcbd_one_var_year.R b/examples/rcbd_one_var_year.R index 8f6da6a..3f2dfa0 100644 --- a/examples/rcbd_one_var_year.R +++ b/examples/rcbd_one_var_year.R @@ -35,7 +35,7 @@ my.data$Year <- as.factor(my.data$Year) cat('Linear mixed effect model results\n') sep(79) #-----------------------------------------------------------------------------# -model <- lme(fixed = Yield ~ Treatment, random = ~1|Block/Year, data = my.data) +model <- lme(fixed = Yield ~ Treatment, random = ~1|Year/Block, data = my.data) summary(model) #-----------------------------------------------------------------------------# sep(79)