Exam3.3.RdExam3.3 use RCBD data with fixed location effect and different forms of estimable functions are shown in this example.
Stroup, W. W. (2012). Generalized Linear Mixed Models: Modern Concepts, Methods and Applications. CRC Press.
#----------------------------------------------------------------------------------- ## linear model for Gaussian data #----------------------------------------------------------------------------------- data(DataSet3.2) DataSet3.2$trt <- factor(x = DataSet3.2$trt, level = c(3,0,1,2)) DataSet3.2$loc <- factor(x = DataSet3.2$loc, level = c(8, 1, 2, 3, 4, 5, 6, 7)) Exam3.3.lm1 <- lm( formula = Y~ trt+loc , data = DataSet3.2 # , subset # , weights # , na.action , method = "qr" , model = TRUE # , x = FALSE # , y = FALSE , qr = TRUE , singular.ok = TRUE , contrasts = NULL # , offset # , ... ) summary( Exam3.3.lm1 )#> #> Call: #> lm(formula = Y ~ trt + loc, data = DataSet3.2, method = "qr", #> model = TRUE, qr = TRUE, singular.ok = TRUE, contrasts = NULL) #> #> Residuals: #> Min 1Q Median 3Q Max #> -2.7750 -0.7875 0.3625 1.0813 2.3000 #> #> Coefficients: #> Estimate Std. Error t value Pr(>|t|) #> (Intercept) 28.2500 0.9945 28.407 < 2e-16 *** #> trt0 -2.1375 0.8481 -2.520 0.01988 * #> trt1 -2.9500 0.8481 -3.478 0.00224 ** #> trt2 0.2875 0.8481 0.339 0.73798 #> loc1 -0.9750 1.1994 -0.813 0.42538 #> loc2 -3.2500 1.1994 -2.710 0.01312 * #> loc3 -1.6750 1.1994 -1.397 0.17713 #> loc4 -0.3500 1.1994 -0.292 0.77329 #> loc5 0.8250 1.1994 0.688 0.49907 #> loc6 -0.3000 1.1994 -0.250 0.80492 #> loc7 -3.5750 1.1994 -2.981 0.00713 ** #> --- #> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 #> #> Residual standard error: 1.696 on 21 degrees of freedom #> Multiple R-squared: 0.6818, Adjusted R-squared: 0.5303 #> F-statistic: 4.5 on 10 and 21 DF, p-value: 0.001795 #>#------------------------------------------------------------- ## Individula least squares treatment means #------------------------------------------------------------- library(lsmeans) (Lsm3.3 <- lsmeans::lsmeans( object = Exam3.3.lm1 , specs = "trt" # , ... ) )#> trt lsmean SE df lower.CL upper.CL #> 3 27.0875 0.5996899 21 25.84038 28.33462 #> 0 24.9500 0.5996899 21 23.70288 26.19712 #> 1 24.1375 0.5996899 21 22.89038 25.38462 #> 2 27.3750 0.5996899 21 26.12788 28.62212 #> #> Results are averaged over the levels of: loc #> Confidence level used: 0.95#--------------------------------------------------- ## Pairwise treatment means estimate #--------------------------------------------------- contrast( object = Lsm3.3 , method = "pairwise")#> contrast estimate SE df t.ratio p.value #> 3 - 0 2.1375 0.8480896 21 2.520 0.0856 #> 3 - 1 2.9500 0.8480896 21 3.478 0.0111 #> 3 - 2 -0.2875 0.8480896 21 -0.339 0.9862 #> 0 - 1 0.8125 0.8480896 21 0.958 0.7742 #> 0 - 2 -2.4250 0.8480896 21 -2.859 0.0430 #> 1 - 2 -3.2375 0.8480896 21 -3.817 0.0051 #> #> Results are averaged over the levels of: loc #> P value adjustment: tukey method for comparing a family of 4 estimates#--------------------------------------------------- ## Repairwise treatment means estimate #--------------------------------------------------- ## contrast( object = Lsm3.3 , method = "repairwise") #------------------------------------------------------- ## LSM Trt0 (This term is used in Walter Stroups' book) #------------------------------------------------------- library(phia) list3.3.1 <- list(trt=c("0" = 1 ) ) Test3.3.1 <- summary(testFactors( model = Exam3.3.lm1 , levels = list3.3.1) ) #------------------------------------------------------- ## LSM Trt0 alt(This term is used in Walter Stroups' book) #------------------------------------------------------- list3.3.2 <- list(trt=c("0" = 1 ) , loc=c("1" = 0,"2" = 0,"3" = 0,"4" = 0,"5" = 0,"6" = 0,"7" = 0) ) Test3.3.2 <- summary(testFactors( model = Exam3.3.lm1 , levels = list3.3.2) ) #------------------------------------------------------- ## Trt0 Vs Trt1 #------------------------------------------------------- list3.3.3 <- list(trt=c("0" = 1,"1" = -1)) Test3.3.3 <- summary(testFactors( model = Exam3.3.lm1 , levels = list3.3.3) ) #------------------------------------------------------- ## average Trt0+1 #------------------------------------------------------- list3.3.4 <- list(trt=c("0" = 0.5 , "1" = 0.5)) Test3.3.4 <- summary(testFactors( model = Exam3.3.lm1 , levels = list3.3.4) ) #------------------------------------------------------- ## average Trt0+2+3 #------------------------------------------------------- list3.3.5 <- list(trt=c("0" = 0.33333,"2" = 0.33333,"3" = 0.33333)) Test3.3.5 <- summary(testFactors( model = Exam3.3.lm1 , levels = list3.3.5) ) #------------------------------------------------------- ## Trt 2 Vs 3 difference #------------------------------------------------------- list3.3.6 <- list(trt=c("2" = 1,"3" = -1)) Test3.3.6 <- summary(testFactors( model = Exam3.3.lm1 , levels = list3.3.6) ) #------------------------------------------------------- ## Trt 1 Vs 2 difference #------------------------------------------------------- list3.3.7 <- list(trt=c("1" = 1,"2" = -1)) Test3.3.7 <- summary(testFactors( model = Exam3.3.lm1 , levels = list3.3.7) ) #------------------------------------------------------- ## Trt 1 Vs 3 difference #------------------------------------------------------- list3.3.8 <- list(trt=c("1" = 1,"3" = -1)) Test3.3.8 <- summary(testFactors( model = Exam3.3.lm1 , levels = list3.3.8) ) #------------------------------------------------------- ## Average trt0+1 vs Average Trt2+3 #------------------------------------------------------- list3.3.9 <- list(trt=c("0" = 0.5,"1" = 0.5,"2" = -0.5,"3" = -0.5)) Test3.3.9 <- summary(testFactors( model = Exam3.3.lm1 , levels = list3.3.9) ) #------------------------------------------------------- ## Trt1 vs Average Trt0+1+2 #------------------------------------------------------- list3.3.10 <- list(trt=c("0" = 0.33333,"1" = -1,"2" = 0.33333,"3" = 0.33333)) Test3.3.10 <- summary(testFactors( model = Exam3.3.lm1 , levels = list3.3.10) ) #------------------------------------------------------- ## Sidak Multiplicity adjustment for p-values #------------------------------------------------------- library(mutoss) PValues3.3 <- c( Test3.3.3[[7]][1, 4] , Test3.3.6[[7]][1, 4] , Test3.3.7[[7]][1, 4] , Test3.3.8[[7]][1, 4] , Test3.3.9[[7]][1, 4] , Test3.3.10[[7]][1, 4] ) AdjPValues3.3 <- sidak(PValues3.3)