,
Adeela Munawar [aut, ctb],
Walter W. Stroup [aut, ctb],
Kent M. Eskridge [aut, ctb]| Title: | R Codes and Datasets for Generalized Linear Mixed Models: Modern Concepts, Methods and Applications by Walter W. Stroup |
|---|---|
| Description: | R Codes and Datasets for Stroup, W. W. (2012). Generalized Linear Mixed Models Modern Concepts, Methods and Applications, CRC Press. |
| Authors: | Muhammad Yaseen [aut, cre, cph]
,
Adeela Munawar [aut, ctb],
Walter W. Stroup [aut, ctb],
Kent M. Eskridge [aut, ctb] |
| Maintainer: | Muhammad Yaseen <[email protected]> |
| License: | GPL-3 |
| Version: | 0.3.0 |
| Built: | 2025-01-30 04:07:57 UTC |
| Source: | https://github.com/myaseen208/stroupglmm |
Exam2.B.2 is used to visualize the effect of glm model statement with binomial data with logit and probit links.
data(DataExam2.B.2)data(DataExam2.B.2)
A data.frame with 11 rows and 3 variables.
x independent variable
n bernouli trials (bernouli outcomes on each individual)
y number of successes on each individual
Muhammad Yaseen ([email protected])
Adeela Munawar ([email protected])
Stroup, W. W. (2012).Generalized linear mixed models: modern concepts, methods and applications. CRC press.
data(DataExam2.B.2)data(DataExam2.B.2)
Exam2.B.3 is used to illustrate one way treatment design with Gaussian observations.
data(DataExam2.B.3)data(DataExam2.B.3)
A data.frame with 6 rows and 2 variables.
trt treatments as factor with number 1 to 3
y response variable
Muhammad Yaseen ([email protected])
Adeela Munawar ([email protected])
Stroup, W. W. (2012).Generalized linear mixed models: modern concepts, methods and applications. CRC press.
data(DataExam2.B.3)data(DataExam2.B.3)
Exam2.B.4 is used to illustrate one way treatment design with Binomial observations.
data(DataExam2.B.4)data(DataExam2.B.4)
A data.frame with 6 rows and 4 variables.
obs number of observations
trt three treatments with class factor
Nij number of bernouli trials on each individual
y number of successes on each individual
Muhammad Yaseen ([email protected])
Adeela Munawar ([email protected])
Stroup, W. W. (2012).Generalized linear mixed models: modern concepts, methods and applications. CRC press.
data(DataExam2.B.4)data(DataExam2.B.4)
Exam2.B.7 is related to multi batch regression data assuming different forms of linear models with factorial experiment.
data(DataExam2.B.7)data(DataExam2.B.7)
A data.frame with 16 rows and 4 variables.
Rep number of replications
a factor with two levels 1 and 2
b factor with two levels 1 and 2
y response variable
Muhammad Yaseen ([email protected])
Adeela Munawar ([email protected])
Stroup, W. W. (2012).Generalized linear mixed models: modern concepts, methods and applications. CRC press.
data(DataExam2.B.7)data(DataExam2.B.7)
DataSet3.1 is used for linear and generalized linear models
data(DataSet3.1)data(DataSet3.1)
A data.frame with 20 rows and 5 variables.
trt two treatment 0 and 1
rep unit of observation or observation ID
Y is continuous & may be assumed Gaussian
N is the number of obs
F is the number of "successes"(N and F specify a binomial response)
Muhammad Yaseen ([email protected])
Adeela Munawar ([email protected])
Stroup, W. W. (2012).Generalized linear mixed models: modern concepts, methods and applications. CRC press.
data(DataSet3.1)data(DataSet3.1)
DataSet3.2 Multi-Location, 4 Treatment Randomized Block
data(DataSet3.2)data(DataSet3.2)
A data.frame with 32 rows and 10 variables.
trt two treatment 0 and 1
loc four locations used as blocks
Y is Gaussian response variable
Nbin subjects at each Loc x Trt for binomial response
S1 and S2 are two binomial response variables
count1 and count 2 used later
A and B are factors with level 0 and 1
Muhammad Yaseen ([email protected])
Adeela Munawar ([email protected])
Stroup, W. W. (2012).Generalized linear mixed models: modern concepts, methods and applications. CRC press.
data(DataSet3.2)data(DataSet3.2)
Exam1.2 is used to see types of model effects by plotting regression data
data(DataSet3.3)data(DataSet3.3)
A data.frame with 36 rows and 6 variables.
X Each batch observed at several times:0,3,6,12,24,36,48 months
Y continuous variable observed at each level of X
Fav number of successes
N isndependent bernoulli trials
Batch Batches as 1, 2, 3, 4
Count binomial response variable
Muhammad Yaseen ([email protected])
Adeela Munawar ([email protected])
Stroup, W. W. (2012).Generalized linear mixed models: modern concepts, methods and applications. CRC press.
data(DataSet3.3)data(DataSet3.3)
DataSet4.1 comes from Cochran and Cox (1957) Experimental Design
data(DataSet4.1)data(DataSet4.1)
A data.frame with 60 rows and 3 variables.
blocks 15 blocks in an incomplete block desgin
trt treatments representing incomplete block desgin
y is continuous & may be assumed Gaussian
Muhammad Yaseen ([email protected])
Adeela Munawar ([email protected])
Stroup, W. W. (2012).Generalized linear mixed models: modern concepts, methods and applications. CRC press.
Cochran, W. G., & Cox, G. M. (1957). Experimental designs.
data(DataSet4.1)data(DataSet4.1)
DataSet5.1 is used for polynomial multiple regression
data(DataSet5.1)data(DataSet5.1)
A data.frame with 14 rows and 3 variables.
X is predictor variable with level 0, 1, 2, 4, 8, 12, 16
N is the number of independent bernoulli trials for a given observation
F is the number of "successes"(N and F specify a binomial response)
Muhammad Yaseen ([email protected])
Adeela Munawar ([email protected])
Stroup, W. W. (2012).Generalized linear mixed models: modern concepts, methods and applications. CRC press.
data(DataSet5.1)data(DataSet5.1)
DataSet5.2 is used for three factor orthogonal main effects only design with sequential fitting of predictors
data(DataSet5.2)data(DataSet5.2)
A data.frame with 9 rows and 4 variables.
a is predictor variable with level 0, 1
b is predictor variable with level 0, 1
c is predictor variable with level 0, 1
y response variable
Muhammad Yaseen ([email protected])
Adeela Munawar ([email protected])
Stroup, W. W. (2012).Generalized linear mixed models: modern concepts, methods and applications. CRC press.
data(DataSet5.2)data(DataSet5.2)
Data for Example 7.1 from Generalized Linear Mixed Models: Modern Concepts, Methods and Applications by Walter W. Stroup
data(DataSet7.1)data(DataSet7.1)
Muhammad Yaseen ([email protected])
Adeela Munawar ([email protected])
Stroup, W. W. (2012).Generalized linear mixed models: modern concepts, methods and applications. CRC press.
data(DataSet7.1)data(DataSet7.1)
Data for Example 7.2 from Generalized Linear Mixed Models: Modern Concepts, Methods and Applications by Walter W. Stroup
data(DataSet7.2)data(DataSet7.2)
Muhammad Yaseen ([email protected])
Adeela Munawar ([email protected])
Stroup, W. W. (2012).Generalized linear mixed models: modern concepts, methods and applications. CRC press.
data(DataSet7.2)data(DataSet7.2)
Data for Example 7.3 from Generalized Linear Mixed Models: Modern Concepts, Methods and Applications by Walter W. Stroup
data(DataSet7.3)data(DataSet7.3)
Muhammad Yaseen ([email protected])
Adeela Munawar ([email protected])
Stroup, W. W. (2012).Generalized linear mixed models: modern concepts, methods and applications. CRC press.
data(DataSet7.3)data(DataSet7.3)
Data for Example 7.4 from Generalized Linear Mixed Models: Modern Concepts, Methods and Applications by Walter W. Stroup
data(DataSet7.4)data(DataSet7.4)
Muhammad Yaseen ([email protected])
Adeela Munawar ([email protected])
Stroup, W. W. (2012).Generalized linear mixed models: modern concepts, methods and applications. CRC press.
data(DataSet7.4)data(DataSet7.4)
Data for Example 7.4 from Generalized Linear Mixed Models: Modern Concepts, Methods and Applications by Walter W. Stroup
data(DataSet7.4rsm)data(DataSet7.4rsm)
Muhammad Yaseen ([email protected])
Adeela Munawar ([email protected])
Stroup, W. W. (2012).Generalized linear mixed models: modern concepts, methods and applications. CRC press.
data(DataSet7.4rsm)data(DataSet7.4rsm)
Data for Example 7.6 from Generalized Linear Mixed Models: Modern Concepts, Methods and Applications by Walter W. Stroup
data(DataSet7.6)data(DataSet7.6)
Muhammad Yaseen ([email protected])
Adeela Munawar ([email protected])
Stroup, W. W. (2012).Generalized linear mixed models: modern concepts, methods and applications. CRC press.
data(DataSet7.6)data(DataSet7.6)
Data for Example 7.7 from Generalized Linear Mixed Models: Modern Concepts, Methods and Applications by Walter W. Stroup
data(DataSet7.7)data(DataSet7.7)
Muhammad Yaseen ([email protected])
Adeela Munawar ([email protected])
Stroup, W. W. (2012).Generalized linear mixed models: modern concepts, methods and applications. CRC press.
data(DataSet7.7)data(DataSet7.7)
DataSet8.1 is used for Nested factorial structure
data(DataSet8.1)data(DataSet8.1)
A data.frame with 30 rows and 4 variables.
block 10 blocks
trt 6 treatments nested within sets
set 2 sets
y is a Gaussian response variable
Muhammad Yaseen ([email protected]) Adeela Munawar ([email protected])
Stroup, W. W. (2012).Generalized linear Mixed Models: Modern Concepts, Methods and Applications. CRC press.
data(DataSet8.1)data(DataSet8.1)
DataSet8.2 is used for Incomplete strip-plot ( 3 cross 3 factorial).
data(DataSet8.2)data(DataSet8.2)
A data.frame with 36 rows and 6 variables.
block 9 blocks each consisting of 2 rows and 2 coloumns
a is a factor with 3 levels assigned at random to rows
b is a factor with 3 levels assigned at random to columns
y is a Gaussian response variable
Muhammad Yaseen ([email protected]) Adeela Munawar ([email protected])
Stroup, W. W. (2012).Generalized linear Mixed Models: Modern Concepts, Methods and Applications. CRC press.
data(DataSet8.2)data(DataSet8.2)
DataSet8.3 is used for Response surface design with incomplete blocking
data(DataSet8.3)data(DataSet8.3)
A data.frame with 28 rows and 4 variables.
block with 7 blocks
a is a factor with 3 levels 0,-1 and 1
b is a factor with 3 levels 0,-1 and 1
c is a factor with 3 levels 0,-1 and 1
y is a Gaussian response variable
Muhammad Yaseen ([email protected]) Adeela Munawar ([email protected])
Stroup, W. W. (2012).Generalized linear Mixed Models: Modern Concepts, Methods and Applications. CRC press.
data(DataSet8.3)data(DataSet8.3)
DataSet8.4 is used for Multifactor treatment and Multilevel design structures
data(DataSet8.4)data(DataSet8.4)
A data.frame with 36 rows and 6 variables.
block 9 blocks each consisting of 2 rows and 2 coloumns
a is a factor with 3 levels assigned at random to rows
b is a factor with 3 levels assigned at random to columns
y is a Gaussian response variable
Muhammad Yaseen ([email protected]) Adeela Munawar ([email protected])
Stroup, W. W. (2012).Generalized linear Mixed Models: Modern Concepts, Methods and Applications. CRC press.
data(DataSet8.4)data(DataSet8.4)
DataSet9.1 is used for One-way random effects only model
data(DataSet9.1)data(DataSet9.1)
A data.frame with 24 rows and 2 variables.
a is a factor with 12 levels
y is a Gaussian response variable
Muhammad Yaseen ([email protected]) Adeela Munawar ([email protected])
Stroup, W. W. (2012).Generalized linear Mixed Models: Modern Concepts, Methods and Applications. CRC press.
data(DataSet9.1)data(DataSet9.1)
DataSet9.2 is used for Two way random effects nested model
data(DataSet9.2)data(DataSet9.2)
A data.frame with 28 rows and 3 variables with levels of b nested within levels of.
a is a factor with 7 levels
b is a factor with 2 levels
y is a Gaussian response variable
Muhammad Yaseen ([email protected]) Adeela Munawar ([email protected])
Stroup, W. W. (2012).Generalized linear Mixed Models: Modern Concepts, Methods and Applications. CRC press.
data(DataSet9.2)data(DataSet9.2)
DataSet9.4 is used for Relationship between BLUP and Fixed Effect Estimators
data(DataSet9.4)data(DataSet9.4)
A data.frame with 32 rows and 3 variables
a is a factor with 2 levels
b is a factor with 8 levels
y is a Gaussian response variable
Muhammad Yaseen ([email protected]) Adeela Munawar ([email protected])
Stroup, W. W. (2012).Generalized linear Mixed Models: Modern Concepts, Methods and Applications. CRC press.
data(DataSet9.4)data(DataSet9.4)
Exam1.1 is used for inspecting probability distribution and to define a plausible process through linear models and generalized linear models.
Muhammad Yaseen ([email protected])
Adeela Munawar ([email protected])
Stroup, W. W. (2012). Generalized Linear Mixed Models: Modern Concepts, Methods and Applications. CRC Press.
#------------------------------------------------------------- ## Linear Model and results discussed in Article 1.2.1 after Table1.1 #------------------------------------------------------------- data(Table1.1) Exam1.1.lm1 <- lm(formula = y/Nx ~ x, data = Table1.1) summary(Exam1.1.lm1 ) library(parameters) model_parameters(Exam1.1.lm1) #------------------------------------------------------------- ## GLM fitting with logit link (family=binomial) #------------------------------------------------------------- Exam1.1.glm1 <- glm( formula = y/Nx ~ x , family = binomial(link = "logit") , data = Table1.1 ) ## this glm() function gives warning message of non-integer success summary(Exam1.1.glm1) model_parameters(Exam1.1.glm1) #------------------------------------------------------------- ## GLM fitting with logit link (family = Quasibinomial) #------------------------------------------------------------- Exam1.1.glm2 <- glm( formula = y/Nx~x , family = quasibinomial(link = "logit") , data = Table1.1 ) ## problem of "warning message of non-integer success" is overome by using quasibinomial family summary(Exam1.1.glm2) model_parameters(Exam1.1.glm2) #------------------------------------------------------------- ## GLM fitting with survey package(produces same result as using quasi binomial family in glm) #------------------------------------------------------------- library(survey) design <- svydesign(ids = ~1, data = Table1.1) Exam1.1.svyglm <- svyglm( formula = y/Nx~x , design = design , family = quasibinomial(link = "logit") ) summary(Exam1.1.svyglm) model_parameters(Exam1.1.svyglm) #------------------------------------------------------------- ## Figure 1.1 #------------------------------------------------------------- Newdata <- data.frame( Table1.1 , LM = Exam1.1.lm1$fitted.values , GLM = Exam1.1.glm1$fitted.values , QB = Exam1.1.glm2$fitted.values , SM = Exam1.1.svyglm$fitted.values ) #------------------------------------------------------------- ## One Method to plot Figure1.1 #------------------------------------------------------------- library(ggplot2) Figure1.1 <- ggplot( data = Newdata , mapping = aes(x = x, y = y/Nx) ) + geom_point ( mapping = aes(colour = "black") ) + geom_point ( data = Newdata , mapping = aes(x = x, y = LM, colour = "blue"), shape = 2 ) + geom_line( data = Newdata , mapping = aes(x = x, y = LM, colour = "blue") ) + geom_point ( data = Newdata , mapping = aes(x = x, y = GLM, colour ="red"), shape = 3 ) + geom_smooth ( data = Newdata , mapping = aes(x = x, y = GLM, colour = "red") , stat = "smooth" ) + theme_bw() + scale_colour_manual ( values = c("black", "blue", "red"), labels = c("observed", "LM", "GLM") ) + guides ( colour = guide_legend(title = "Plot") ) + labs ( title = "Linear Model vs Logistic Model" ) + labs ( y = "p" ) print(Figure1.1) #------------------------------------------------------------- ## Another way to plot Figure 1.1 #------------------------------------------------------------- newdata <- data.frame( P = c( Table1.1$y/Table1.1$Nx , Exam1.1.lm1$fitted.values , Exam1.1.glm1$fitted.values ) , X = rep(Table1.1$x, 3) , group = rep(c('Obs','LM','GLM'), each = length(Table1.1$x)) ) Figure1.1 <- ggplot( data = newdata , mapping = aes(x = X , y = P) ) + geom_point( mapping = aes(x = X , y = P, colour = group , shape=group) ) + geom_smooth( data = subset(x = newdata, group == "LM") , mapping = aes(x=X,y=P) , col = "green" ) + geom_smooth( data = subset(x = newdata, group=="GLM") , mapping = aes(x = X , y = P) , col = "red" ) + theme_bw() + labs( title = "Linear Model vs Logistic Model" ) print(Figure1.1) #------------------------------------------------------------- ## Correlation among p and fitted values using Gaussian link #------------------------------------------------------------- (lmCor <- cor(Table1.1$y/Table1.1$Nx, Exam1.1.lm1$fitted.values)) #------------------------------------------------------------- ## Correlation among p and fitted values using quasi binomial link #------------------------------------------------------------- (glmCor <- cor(Table1.1$y/Table1.1$Nx, Exam1.1.glm1$fitted.values))#------------------------------------------------------------- ## Linear Model and results discussed in Article 1.2.1 after Table1.1 #------------------------------------------------------------- data(Table1.1) Exam1.1.lm1 <- lm(formula = y/Nx ~ x, data = Table1.1) summary(Exam1.1.lm1 ) library(parameters) model_parameters(Exam1.1.lm1) #------------------------------------------------------------- ## GLM fitting with logit link (family=binomial) #------------------------------------------------------------- Exam1.1.glm1 <- glm( formula = y/Nx ~ x , family = binomial(link = "logit") , data = Table1.1 ) ## this glm() function gives warning message of non-integer success summary(Exam1.1.glm1) model_parameters(Exam1.1.glm1) #------------------------------------------------------------- ## GLM fitting with logit link (family = Quasibinomial) #------------------------------------------------------------- Exam1.1.glm2 <- glm( formula = y/Nx~x , family = quasibinomial(link = "logit") , data = Table1.1 ) ## problem of "warning message of non-integer success" is overome by using quasibinomial family summary(Exam1.1.glm2) model_parameters(Exam1.1.glm2) #------------------------------------------------------------- ## GLM fitting with survey package(produces same result as using quasi binomial family in glm) #------------------------------------------------------------- library(survey) design <- svydesign(ids = ~1, data = Table1.1) Exam1.1.svyglm <- svyglm( formula = y/Nx~x , design = design , family = quasibinomial(link = "logit") ) summary(Exam1.1.svyglm) model_parameters(Exam1.1.svyglm) #------------------------------------------------------------- ## Figure 1.1 #------------------------------------------------------------- Newdata <- data.frame( Table1.1 , LM = Exam1.1.lm1$fitted.values , GLM = Exam1.1.glm1$fitted.values , QB = Exam1.1.glm2$fitted.values , SM = Exam1.1.svyglm$fitted.values ) #------------------------------------------------------------- ## One Method to plot Figure1.1 #------------------------------------------------------------- library(ggplot2) Figure1.1 <- ggplot( data = Newdata , mapping = aes(x = x, y = y/Nx) ) + geom_point ( mapping = aes(colour = "black") ) + geom_point ( data = Newdata , mapping = aes(x = x, y = LM, colour = "blue"), shape = 2 ) + geom_line( data = Newdata , mapping = aes(x = x, y = LM, colour = "blue") ) + geom_point ( data = Newdata , mapping = aes(x = x, y = GLM, colour ="red"), shape = 3 ) + geom_smooth ( data = Newdata , mapping = aes(x = x, y = GLM, colour = "red") , stat = "smooth" ) + theme_bw() + scale_colour_manual ( values = c("black", "blue", "red"), labels = c("observed", "LM", "GLM") ) + guides ( colour = guide_legend(title = "Plot") ) + labs ( title = "Linear Model vs Logistic Model" ) + labs ( y = "p" ) print(Figure1.1) #------------------------------------------------------------- ## Another way to plot Figure 1.1 #------------------------------------------------------------- newdata <- data.frame( P = c( Table1.1$y/Table1.1$Nx , Exam1.1.lm1$fitted.values , Exam1.1.glm1$fitted.values ) , X = rep(Table1.1$x, 3) , group = rep(c('Obs','LM','GLM'), each = length(Table1.1$x)) ) Figure1.1 <- ggplot( data = newdata , mapping = aes(x = X , y = P) ) + geom_point( mapping = aes(x = X , y = P, colour = group , shape=group) ) + geom_smooth( data = subset(x = newdata, group == "LM") , mapping = aes(x=X,y=P) , col = "green" ) + geom_smooth( data = subset(x = newdata, group=="GLM") , mapping = aes(x = X , y = P) , col = "red" ) + theme_bw() + labs( title = "Linear Model vs Logistic Model" ) print(Figure1.1) #------------------------------------------------------------- ## Correlation among p and fitted values using Gaussian link #------------------------------------------------------------- (lmCor <- cor(Table1.1$y/Table1.1$Nx, Exam1.1.lm1$fitted.values)) #------------------------------------------------------------- ## Correlation among p and fitted values using quasi binomial link #------------------------------------------------------------- (glmCor <- cor(Table1.1$y/Table1.1$Nx, Exam1.1.glm1$fitted.values))
Exam1.2 is used to see types of model effects by plotting regression data
Muhammad Yaseen ([email protected])
Adeela Munawar ([email protected])
Stroup, W. W. (2012). Generalized Linear Mixed Models: Modern Concepts, Methods and Applications. CRC Press.
#------------------------------------------------------------- ## Plot of multi-batch regression data discussed in Article 1.3 #------------------------------------------------------------- data(Table1.1) Table1.2$Batch <- factor(x = Table1.2$Batch) library(ggplot2) Plot <- ggplot(data = Table1.2, mapping = aes(y = Y, x = X, colour = Batch, shape = Batch)) + geom_point() + geom_smooth(method = "lm", fill = NA) + labs(title = "Plot of Multi Batch Regression data") + theme_bw() Plot#------------------------------------------------------------- ## Plot of multi-batch regression data discussed in Article 1.3 #------------------------------------------------------------- data(Table1.1) Table1.2$Batch <- factor(x = Table1.2$Batch) library(ggplot2) Plot <- ggplot(data = Table1.2, mapping = aes(y = Y, x = X, colour = Batch, shape = Batch)) + geom_point() + geom_smooth(method = "lm", fill = NA) + labs(title = "Plot of Multi Batch Regression data") + theme_bw() Plot
Exam2.B.1 is used to visualize the effect of lm model statement with Gaussian data and their design matrix
Muhammad Yaseen ([email protected])
Adeela Munawar ([email protected])
Stroup, W. W. (2012). Generalized Linear Mixed Models: Modern Concepts, Methods and Applications. CRC Press.
#----------------------------------------------------------------------------------- ## Linear Model discussed in Example 2.B.1 using simple regression data of Table1.1 #----------------------------------------------------------------------------------- data(Table1.1) Exam2.B.1.lm1 <- lm(formula = y~x, data = Table1.1) summary(Exam2.B.1.lm1) library(parameters) model_parameters(Exam2.B.1.lm1) DesignMatrix.lm1 <- model.matrix (object = Exam2.B.1.lm1) DesignMatrix.lm1#----------------------------------------------------------------------------------- ## Linear Model discussed in Example 2.B.1 using simple regression data of Table1.1 #----------------------------------------------------------------------------------- data(Table1.1) Exam2.B.1.lm1 <- lm(formula = y~x, data = Table1.1) summary(Exam2.B.1.lm1) library(parameters) model_parameters(Exam2.B.1.lm1) DesignMatrix.lm1 <- model.matrix (object = Exam2.B.1.lm1) DesignMatrix.lm1
Exam2.B.2 is used to visualize the effect of glm model statement with binomial data with logit and probit links.
Muhammad Yaseen ([email protected])
Adeela Munawar ([email protected])
Stroup, W. W. (2012). Generalized Linear Mixed Models: Modern Concepts, Methods and Applications. CRC Press.
#----------------------------------------------------------------------------------- ## probitit Model discussed in Example 2.B.2 using DataExam2.B.2 ## Default link is logit ## using fmaily = binomial gives warning message of no-integer successes #----------------------------------------------------------------------------------- data(DataExam2.B.2) Exam2.B.2glm <- glm(formula = y/n~x, family = quasibinomial(link = "probit"), data = DataExam2.B.2) summary(Exam2.B.2glm) library(parameters) model_parameters(Exam2.B.2glm)#----------------------------------------------------------------------------------- ## probitit Model discussed in Example 2.B.2 using DataExam2.B.2 ## Default link is logit ## using fmaily = binomial gives warning message of no-integer successes #----------------------------------------------------------------------------------- data(DataExam2.B.2) Exam2.B.2glm <- glm(formula = y/n~x, family = quasibinomial(link = "probit"), data = DataExam2.B.2) summary(Exam2.B.2glm) library(parameters) model_parameters(Exam2.B.2glm)
Exam2.B.3 is used to illustrate one way treatment design with Gaussian observations.
Muhammad Yaseen ([email protected])
Adeela Munawar ([email protected])
Stroup, W. W. (2012). Generalized Linear Mixed Models: Modern Concepts, Methods and Applications. CRC Press.
#----------------------------------------------------------------------------------- ## Means Model discussed in Example 2.B.3 using DataExam2.B.3 #----------------------------------------------------------------------------------- Exam2.B.3.lm1 <- lm(formula = y ~ trt, data = DataExam2.B.3) summary(Exam2.B.3.lm1) #----------------------------------------------------------------------------------- ## Effectss Model discussed in Example 2.B.3 using DataExam2.B.3 #----------------------------------------------------------------------------------- Exam2.B.3.lm2 <- lm(formula = y ~ 0 + trt, data = DataExam2.B.3) summary(Exam2.B.3.lm2) library(parameters) model_parameters(Exam2.B.3.lm2)#----------------------------------------------------------------------------------- ## Means Model discussed in Example 2.B.3 using DataExam2.B.3 #----------------------------------------------------------------------------------- Exam2.B.3.lm1 <- lm(formula = y ~ trt, data = DataExam2.B.3) summary(Exam2.B.3.lm1) #----------------------------------------------------------------------------------- ## Effectss Model discussed in Example 2.B.3 using DataExam2.B.3 #----------------------------------------------------------------------------------- Exam2.B.3.lm2 <- lm(formula = y ~ 0 + trt, data = DataExam2.B.3) summary(Exam2.B.3.lm2) library(parameters) model_parameters(Exam2.B.3.lm2)
Exam2.B.4 is used to illustrate one way treatment design with Binomial observations.
Muhammad Yaseen ([email protected])
Adeela Munawar ([email protected])
Stroup, W. W. (2012). Generalized Linear Mixed Models: Modern Concepts, Methods and Applications. CRC Press.
#----------------------------------------------------------------------------------- ## logit Model discussed in Example 2.B.2 using DataExam2.B.4 ## Default link is logit ## using fmaily=binomial gives warning message of no-integer successes #----------------------------------------------------------------------------------- data(DataExam2.B.4) DataExam2.B.4$trt <- factor(x = DataExam2.B.4$trt) Exam2.B.4glm <- glm( formula = Yij/Nij ~ trt , family = quasibinomial(link = "probit") , data = DataExam2.B.4 ) summary(Exam2.B.4glm) library(parameters) model_parameters(Exam2.B.4glm)#----------------------------------------------------------------------------------- ## logit Model discussed in Example 2.B.2 using DataExam2.B.4 ## Default link is logit ## using fmaily=binomial gives warning message of no-integer successes #----------------------------------------------------------------------------------- data(DataExam2.B.4) DataExam2.B.4$trt <- factor(x = DataExam2.B.4$trt) Exam2.B.4glm <- glm( formula = Yij/Nij ~ trt , family = quasibinomial(link = "probit") , data = DataExam2.B.4 ) summary(Exam2.B.4glm) library(parameters) model_parameters(Exam2.B.4glm)
Exam2.B.5 is related to multi batch regression data assuming different forms of linear models.
Muhammad Yaseen ([email protected])
Adeela Munawar ([email protected])
Stroup, W. W. (2012). Generalized Linear Mixed Models: Modern Concepts, Methods and Applications. CRC Press.
#----------------------------------------------------------------------------------- ## Nested Model with no intercept #----------------------------------------------------------------------------------- data(Table1.2) Table1.2$Batch <- factor(x = Table1.2$Batch) Exam2.B.5.lm1 <- lm(formula = Y ~ 0 + Batch + Batch/X, data = Table1.2) DesignMatrix.lm1 <- model.matrix (object = Exam2.B.5.lm1) DesignMatrix.lm1 #----------------------------------------------------------------------------------- ## Interaction Model with intercept #----------------------------------------------------------------------------------- Exam2.B.5.lm2 <-lm(formula = Y ~ Batch + X + Batch*X, data = Table1.2) DesignMatrix.lm2 <- model.matrix (object = Exam2.B.5.lm2) DesignMatrix.lm2 #----------------------------------------------------------------------------------- ## Interaction Model with no intercept #----------------------------------------------------------------------------------- Exam2.B.5.lm3 <- lm(formula = Y ~ 0 + Batch + Batch*X, data = Table1.2) DesignMatrix.lm3 <- model.matrix(object = Exam2.B.5.lm3) DesignMatrix.lm3 #----------------------------------------------------------------------------------- ## Interaction Model with intercept but omitting X term as main effect #----------------------------------------------------------------------------------- Exam2.B.5.lm4 <- lm(formula = Y ~ Batch + Batch*X, data = Table1.2) DesignMatrix.lm4 <- model.matrix(object = Exam2.B.5.lm4) DesignMatrix.lm4#----------------------------------------------------------------------------------- ## Nested Model with no intercept #----------------------------------------------------------------------------------- data(Table1.2) Table1.2$Batch <- factor(x = Table1.2$Batch) Exam2.B.5.lm1 <- lm(formula = Y ~ 0 + Batch + Batch/X, data = Table1.2) DesignMatrix.lm1 <- model.matrix (object = Exam2.B.5.lm1) DesignMatrix.lm1 #----------------------------------------------------------------------------------- ## Interaction Model with intercept #----------------------------------------------------------------------------------- Exam2.B.5.lm2 <-lm(formula = Y ~ Batch + X + Batch*X, data = Table1.2) DesignMatrix.lm2 <- model.matrix (object = Exam2.B.5.lm2) DesignMatrix.lm2 #----------------------------------------------------------------------------------- ## Interaction Model with no intercept #----------------------------------------------------------------------------------- Exam2.B.5.lm3 <- lm(formula = Y ~ 0 + Batch + Batch*X, data = Table1.2) DesignMatrix.lm3 <- model.matrix(object = Exam2.B.5.lm3) DesignMatrix.lm3 #----------------------------------------------------------------------------------- ## Interaction Model with intercept but omitting X term as main effect #----------------------------------------------------------------------------------- Exam2.B.5.lm4 <- lm(formula = Y ~ Batch + Batch*X, data = Table1.2) DesignMatrix.lm4 <- model.matrix(object = Exam2.B.5.lm4) DesignMatrix.lm4
Exam2.B.6 is related to multi batch regression data assuming different forms of linear models keeping batch effect random.
Muhammad Yaseen ([email protected])
Adeela Munawar ([email protected])
Stroup, W. W. (2012). Generalized Linear Mixed Models: Modern Concepts, Methods and Applications. CRC Press.
#----------------------------------------------------------------------------------- ## Nested Model with no intercept #----------------------------------------------------------------------------------- data(Table1.2) Table1.2$Batch <- factor(x = Table1.2$Batch) library(nlme) Exam2.B.6fm1 <- lme( fixed = Y ~ X , data = Table1.2 , random = list(Batch = pdDiag(~1), X = pdDiag(~1)) , method = c("REML", "ML")[1] ) Exam2.B.6fm1 library(broom.mixed) tidy(Exam2.B.6fm1)#----------------------------------------------------------------------------------- ## Nested Model with no intercept #----------------------------------------------------------------------------------- data(Table1.2) Table1.2$Batch <- factor(x = Table1.2$Batch) library(nlme) Exam2.B.6fm1 <- lme( fixed = Y ~ X , data = Table1.2 , random = list(Batch = pdDiag(~1), X = pdDiag(~1)) , method = c("REML", "ML")[1] ) Exam2.B.6fm1 library(broom.mixed) tidy(Exam2.B.6fm1)
Exam2.B.7 is related to multi batch regression data assuming different forms of linear models with factorial experiment.
Muhammad Yaseen ([email protected])
Adeela Munawar ([email protected])
Stroup, W. W. (2012). Generalized Linear Mixed Models: Modern Concepts, Methods and Applications. CRC Press.
#----------------------------------------------------------------------------------- ## Classical main effects and Interaction Model #----------------------------------------------------------------------------------- data(DataExam2.B.7) DataExam2.B.7$a <- factor(x = DataExam2.B.7$a) DataExam2.B.7$b <- factor(x = DataExam2.B.7$b) Exam2.B.7.lm1 <- lm(formula = y~ a + b + a*b, data = DataExam2.B.7) #----------------------------------------------------------------------------------- ## One way treatment effects model #----------------------------------------------------------------------------------- DesignMatrix.lm1 <- model.matrix (object = Exam2.B.7.lm1) DesignMatrix2.B.7.2 <- DesignMatrix.lm1[,!colnames(DesignMatrix.lm1) %in% c("a2","b")] lmfit2 <- lm.fit(x = DesignMatrix2.B.7.2, y = DataExam2.B.7$y) Coefficientslmfit2 <- coef( object = lmfit2) Coefficientslmfit2 #----------------------------------------------------------------------------------- ## One way treatment effects model without intercept #----------------------------------------------------------------------------------- DesignMatrix2.B.7.3 <- as.matrix(DesignMatrix.lm1[,!colnames(DesignMatrix.lm1) %in% c("(Intercept)","a2","b")]) lmfit3 <- lm.fit(x = DesignMatrix2.B.7.3, y = DataExam2.B.7$y) Coefficientslmfit3 <- coef( object = lmfit3) Coefficientslmfit3 #----------------------------------------------------------------------------------- ## Nested Model (both models give the same result) #----------------------------------------------------------------------------------- Exam2.B.7.lm4 <- lm(formula = y~ a + a/b, data = DataExam2.B.7) summary(Exam2.B.7.lm4) Exam2.B.7.lm4 <- lm(formula = y~ a + a*b, data = DataExam2.B.7) summary(Exam2.B.7.lm4)#----------------------------------------------------------------------------------- ## Classical main effects and Interaction Model #----------------------------------------------------------------------------------- data(DataExam2.B.7) DataExam2.B.7$a <- factor(x = DataExam2.B.7$a) DataExam2.B.7$b <- factor(x = DataExam2.B.7$b) Exam2.B.7.lm1 <- lm(formula = y~ a + b + a*b, data = DataExam2.B.7) #----------------------------------------------------------------------------------- ## One way treatment effects model #----------------------------------------------------------------------------------- DesignMatrix.lm1 <- model.matrix (object = Exam2.B.7.lm1) DesignMatrix2.B.7.2 <- DesignMatrix.lm1[,!colnames(DesignMatrix.lm1) %in% c("a2","b")] lmfit2 <- lm.fit(x = DesignMatrix2.B.7.2, y = DataExam2.B.7$y) Coefficientslmfit2 <- coef( object = lmfit2) Coefficientslmfit2 #----------------------------------------------------------------------------------- ## One way treatment effects model without intercept #----------------------------------------------------------------------------------- DesignMatrix2.B.7.3 <- as.matrix(DesignMatrix.lm1[,!colnames(DesignMatrix.lm1) %in% c("(Intercept)","a2","b")]) lmfit3 <- lm.fit(x = DesignMatrix2.B.7.3, y = DataExam2.B.7$y) Coefficientslmfit3 <- coef( object = lmfit3) Coefficientslmfit3 #----------------------------------------------------------------------------------- ## Nested Model (both models give the same result) #----------------------------------------------------------------------------------- Exam2.B.7.lm4 <- lm(formula = y~ a + a/b, data = DataExam2.B.7) summary(Exam2.B.7.lm4) Exam2.B.7.lm4 <- lm(formula = y~ a + a*b, data = DataExam2.B.7) summary(Exam2.B.7.lm4)
Exam3.2 used binomial data, two treatment samples
Muhammad Yaseen ([email protected])
Adeela Munawar ([email protected])
Stroup, W. W. (2012). Generalized Linear Mixed Models: Modern Concepts, Methods and Applications. CRC Press.
#------------------------------------------------------------- ## Linear Model and results discussed in Article 1.2.1 after Table1.1 #------------------------------------------------------------- data(DataSet3.1) DataSet3.1$trt <- factor(x = DataSet3.1$trt) Exam3.2.glm <- glm(formula = F/N~trt, family = quasibinomial(link = "logit"), data = DataSet3.1) summary(Exam3.2.glm) library(parameters) model_parameters(Exam3.2.glm) #------------------------------------------------------------- ## Individula least squares treatment means #------------------------------------------------------------- library(emmeans) emmeans(object = Exam3.2.glm, specs = "trt") emmeans(object = Exam3.2.glm, specs = "trt", type = "response") #--------------------------------------------------- ## Over all mean #--------------------------------------------------- library(phia) list3.2 <- list(trt = c("0" = 0.5, "1" = 0.5 )) testFactors(model = Exam3.2.glm, levels = list3.2 ) #--------------------------------------------------- ## Repairwise treatment means estimate #--------------------------------------------------- contrast(emmeans(object = Exam3.2.glm, specs = "trt")) contrast(emmeans(object = Exam3.2.glm, specs = "trt", type = "response"))#------------------------------------------------------------- ## Linear Model and results discussed in Article 1.2.1 after Table1.1 #------------------------------------------------------------- data(DataSet3.1) DataSet3.1$trt <- factor(x = DataSet3.1$trt) Exam3.2.glm <- glm(formula = F/N~trt, family = quasibinomial(link = "logit"), data = DataSet3.1) summary(Exam3.2.glm) library(parameters) model_parameters(Exam3.2.glm) #------------------------------------------------------------- ## Individula least squares treatment means #------------------------------------------------------------- library(emmeans) emmeans(object = Exam3.2.glm, specs = "trt") emmeans(object = Exam3.2.glm, specs = "trt", type = "response") #--------------------------------------------------- ## Over all mean #--------------------------------------------------- library(phia) list3.2 <- list(trt = c("0" = 0.5, "1" = 0.5 )) testFactors(model = Exam3.2.glm, levels = list3.2 ) #--------------------------------------------------- ## Repairwise treatment means estimate #--------------------------------------------------- contrast(emmeans(object = Exam3.2.glm, specs = "trt")) contrast(emmeans(object = Exam3.2.glm, specs = "trt", type = "response"))
Exam3.3 use RCBD data with fixed location effect and different forms of estimable functions are shown in this example.
Muhammad Yaseen ([email protected])
Adeela Munawar ([email protected])
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)) levels(DataSet3.2$trt) levels(DataSet3.2$loc) Exam3.3.lm1 <- lm(formula = Y~ trt + loc, data = DataSet3.2) summary( Exam3.3.lm1 ) #------------------------------------------------------------- ## Individual least squares treatment means #------------------------------------------------------------- library(emmeans) (Lsm3.3 <- emmeans(object = Exam3.3.lm1, specs = ~trt)) #--------------------------------------------------- ## Pairwise treatment means estimate #--------------------------------------------------- contrast(object = Lsm3.3 , method = "pairwise") #--------------------------------------------------- ## Revpairwise treatment means estimate #--------------------------------------------------- contrast(object = Lsm3.3, method = "revpairwise") #------------------------------------------------------- ## LSM Trt0 (This term is used in Walter Stroups' book) #------------------------------------------------------- contrast( object = emmeans(object = Exam3.3.lm1, specs = ~ trt) , list(trt = c(0, 1, 0, 0)) ) library(phia) testFactors(model = Exam3.3.lm1, levels = list(trt = c("0" = 1))) #------------------------------------------------------- ## LSM Trt0 alt(This term is used in Walter Stroups' book) #------------------------------------------------------- # contrast( # object = emmeans(object = Exam3.3.lm1, specs = ~ trt + loc) # , list( # trt = c(0, 1, 0, 0) # , loc = c(1, 0, 0, 0, 0, 0, 0, 0) # ) # ) # # # list3.3.2 <- # list( # trt = c("0" = 1 ) # , loc = c("1" = 0, "2" = 0,"3" = 0,"4" = 0,"5" = 0,"6" = 0,"7" = 0) # ) # testFactors(model = Exam3.3.lm1, levels = list3.3.2) #------------------------------------------------------- ## Trt0 Vs Trt1 #------------------------------------------------------- contrast( emmeans(object = Exam3.3.lm1, specs = ~trt) , list(trt = c(0, 1, -1, 0)) ) testFactors(model = Exam3.3.lm1, levels = list(trt = c("0" = 1, "1" = -1))) #------------------------------------------------------- ## average Trt0 + Trt1 #------------------------------------------------------- contrast( emmeans(object = Exam3.3.lm1, specs = ~trt) , list(trt = c(0, 1/2, 1/2, 0)) ) testFactors(model = Exam3.3.lm1, levels = list(trt = c("0" = 0.5 , "1" = 0.5))) #------------------------------------------------------- ## average Trt0+2+3 #------------------------------------------------------- contrast( emmeans(object = Exam3.3.lm1, specs = ~trt) , list(trt = c(1/3, 1/3, 0, 1/3)) ) testFactors(model = Exam3.3.lm1, levels = list(trt = c("0" = 1/3,"2" = 1/3,"3" = 1/3))) #------------------------------------------------------- ## Trt 2 Vs 3 difference #------------------------------------------------------- contrast( emmeans(object = Exam3.3.lm1, specs = ~trt) , list(trt = c(-1, 0, 0, 1)) ) testFactors(model = Exam3.3.lm1, levels = list(trt = c("2" = 1,"3" = -1))) #------------------------------------------------------- ## Trt 1 Vs 2 difference #------------------------------------------------------- contrast( emmeans(object = Exam3.3.lm1, specs = ~trt) , list(trt = c(0, 0, 1, -1)) ) testFactors(model = Exam3.3.lm1, levels = list(trt = c("1" = 1,"2" = -1))) #------------------------------------------------------- ## Trt 1 Vs 3 difference #------------------------------------------------------- contrast( emmeans(object = Exam3.3.lm1, specs = ~trt) , list(trt = c(-1, 0, 1, 0)) ) testFactors(model = Exam3.3.lm1, levels = list(trt = c("1" = 1,"3" = -1))) #------------------------------------------------------- ## Average trt0+1 vs Average Trt2+3 #------------------------------------------------------- contrast( emmeans(object = Exam3.3.lm1, specs = ~trt) , list(trt = c(-1/2, 1/2, 1/2, -1/2)) ) testFactors(model = Exam3.3.lm1, levels = list(trt = c("0" = 0.5,"1" = 0.5,"2" = -0.5,"3" = -0.5))) #------------------------------------------------------- ## Trt1 vs Average Trt0+1+2 #------------------------------------------------------- contrast( emmeans(object = Exam3.3.lm1, specs = ~trt) , list(trt = c(1/3, 1/3, -1, 1/3)) ) testFactors(model = Exam3.3.lm1, levels = list(trt = c("0" = 1/3,"1" = -1,"2" = 1/3,"3" = 1/3)))#----------------------------------------------------------------------------------- ## 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)) levels(DataSet3.2$trt) levels(DataSet3.2$loc) Exam3.3.lm1 <- lm(formula = Y~ trt + loc, data = DataSet3.2) summary( Exam3.3.lm1 ) #------------------------------------------------------------- ## Individual least squares treatment means #------------------------------------------------------------- library(emmeans) (Lsm3.3 <- emmeans(object = Exam3.3.lm1, specs = ~trt)) #--------------------------------------------------- ## Pairwise treatment means estimate #--------------------------------------------------- contrast(object = Lsm3.3 , method = "pairwise") #--------------------------------------------------- ## Revpairwise treatment means estimate #--------------------------------------------------- contrast(object = Lsm3.3, method = "revpairwise") #------------------------------------------------------- ## LSM Trt0 (This term is used in Walter Stroups' book) #------------------------------------------------------- contrast( object = emmeans(object = Exam3.3.lm1, specs = ~ trt) , list(trt = c(0, 1, 0, 0)) ) library(phia) testFactors(model = Exam3.3.lm1, levels = list(trt = c("0" = 1))) #------------------------------------------------------- ## LSM Trt0 alt(This term is used in Walter Stroups' book) #------------------------------------------------------- # contrast( # object = emmeans(object = Exam3.3.lm1, specs = ~ trt + loc) # , list( # trt = c(0, 1, 0, 0) # , loc = c(1, 0, 0, 0, 0, 0, 0, 0) # ) # ) # # # list3.3.2 <- # list( # trt = c("0" = 1 ) # , loc = c("1" = 0, "2" = 0,"3" = 0,"4" = 0,"5" = 0,"6" = 0,"7" = 0) # ) # testFactors(model = Exam3.3.lm1, levels = list3.3.2) #------------------------------------------------------- ## Trt0 Vs Trt1 #------------------------------------------------------- contrast( emmeans(object = Exam3.3.lm1, specs = ~trt) , list(trt = c(0, 1, -1, 0)) ) testFactors(model = Exam3.3.lm1, levels = list(trt = c("0" = 1, "1" = -1))) #------------------------------------------------------- ## average Trt0 + Trt1 #------------------------------------------------------- contrast( emmeans(object = Exam3.3.lm1, specs = ~trt) , list(trt = c(0, 1/2, 1/2, 0)) ) testFactors(model = Exam3.3.lm1, levels = list(trt = c("0" = 0.5 , "1" = 0.5))) #------------------------------------------------------- ## average Trt0+2+3 #------------------------------------------------------- contrast( emmeans(object = Exam3.3.lm1, specs = ~trt) , list(trt = c(1/3, 1/3, 0, 1/3)) ) testFactors(model = Exam3.3.lm1, levels = list(trt = c("0" = 1/3,"2" = 1/3,"3" = 1/3))) #------------------------------------------------------- ## Trt 2 Vs 3 difference #------------------------------------------------------- contrast( emmeans(object = Exam3.3.lm1, specs = ~trt) , list(trt = c(-1, 0, 0, 1)) ) testFactors(model = Exam3.3.lm1, levels = list(trt = c("2" = 1,"3" = -1))) #------------------------------------------------------- ## Trt 1 Vs 2 difference #------------------------------------------------------- contrast( emmeans(object = Exam3.3.lm1, specs = ~trt) , list(trt = c(0, 0, 1, -1)) ) testFactors(model = Exam3.3.lm1, levels = list(trt = c("1" = 1,"2" = -1))) #------------------------------------------------------- ## Trt 1 Vs 3 difference #------------------------------------------------------- contrast( emmeans(object = Exam3.3.lm1, specs = ~trt) , list(trt = c(-1, 0, 1, 0)) ) testFactors(model = Exam3.3.lm1, levels = list(trt = c("1" = 1,"3" = -1))) #------------------------------------------------------- ## Average trt0+1 vs Average Trt2+3 #------------------------------------------------------- contrast( emmeans(object = Exam3.3.lm1, specs = ~trt) , list(trt = c(-1/2, 1/2, 1/2, -1/2)) ) testFactors(model = Exam3.3.lm1, levels = list(trt = c("0" = 0.5,"1" = 0.5,"2" = -0.5,"3" = -0.5))) #------------------------------------------------------- ## Trt1 vs Average Trt0+1+2 #------------------------------------------------------- contrast( emmeans(object = Exam3.3.lm1, specs = ~trt) , list(trt = c(1/3, 1/3, -1, 1/3)) ) testFactors(model = Exam3.3.lm1, levels = list(trt = c("0" = 1/3,"1" = -1,"2" = 1/3,"3" = 1/3)))
Exam3.5 fixed location, factorial treatment structure, Gaussian response
Muhammad Yaseen ([email protected])
Adeela Munawar ([email protected])
Stroup, W. W. (2012). Generalized Linear Mixed Models: Modern Concepts, Methods and Applications. CRC Press.
data(DataSet3.2) DataSet3.2$A <- factor(x = DataSet3.2$A) DataSet3.2$B <- factor(x = DataSet3.2$B) DataSet3.2$loc <- factor(x = DataSet3.2$loc, level = c(8, 1, 2, 3, 4, 5, 6, 7)) Exam3.5.lm <- lm(formula = Y~ A + B +loc, data = DataSet3.2) Exam3.5.lm ##---a0 marginal mean library(emmeans) contrast( object = emmeans(object = Exam3.5.lm, specs = ~ B) , list(trt = c(1, 0)) ) library(phia) testFactors(model = Exam3.5.lm, levels = list(B = c("0" = 1,"1" = 0) )) ##---b0 marginal mean testFactors(model = Exam3.5.lm, levels=list(B = c("0" = 1, "1" = 0))) ##---Simple effect of A on B0 testInteractions(model = Exam3.5.lm, custom = list(B = c("0" = 1,"1" = 0)), across = "B") ##---Simple effect of B on A0 testInteractions(model = Exam3.5.lm, custom = list(A = c("0" = 1, "1" = 0)), across = "A") ##---Simple Effect of A over B testInteractions(model = Exam3.5.lm, fixed = "A", across = "B") ##---Simple Effect of B over A testInteractions(model = Exam3.5.lm, fixed = "B", across = "A") #------------------------------------------------------------- ## Individula least squares treatment means #------------------------------------------------------------- emmeans(object = Exam3.5.lm, specs = ~A*B)data(DataSet3.2) DataSet3.2$A <- factor(x = DataSet3.2$A) DataSet3.2$B <- factor(x = DataSet3.2$B) DataSet3.2$loc <- factor(x = DataSet3.2$loc, level = c(8, 1, 2, 3, 4, 5, 6, 7)) Exam3.5.lm <- lm(formula = Y~ A + B +loc, data = DataSet3.2) Exam3.5.lm ##---a0 marginal mean library(emmeans) contrast( object = emmeans(object = Exam3.5.lm, specs = ~ B) , list(trt = c(1, 0)) ) library(phia) testFactors(model = Exam3.5.lm, levels = list(B = c("0" = 1,"1" = 0) )) ##---b0 marginal mean testFactors(model = Exam3.5.lm, levels=list(B = c("0" = 1, "1" = 0))) ##---Simple effect of A on B0 testInteractions(model = Exam3.5.lm, custom = list(B = c("0" = 1,"1" = 0)), across = "B") ##---Simple effect of B on A0 testInteractions(model = Exam3.5.lm, custom = list(A = c("0" = 1, "1" = 0)), across = "A") ##---Simple Effect of A over B testInteractions(model = Exam3.5.lm, fixed = "A", across = "B") ##---Simple Effect of B over A testInteractions(model = Exam3.5.lm, fixed = "B", across = "A") #------------------------------------------------------------- ## Individula least squares treatment means #------------------------------------------------------------- emmeans(object = Exam3.5.lm, specs = ~A*B)
Exam3.9 used to differentiate conditional and marginal binomial models with and without interaction for S2 variable.
Muhammad Yaseen ([email protected])
Adeela Munawar ([email protected])
Stroup, W. W. (2012). Generalized Linear Mixed Models: Modern Concepts, Methods and Applications. CRC Press.
#----------------------------------------------------------------------------------- ## Binomial conditional GLMM without interaction, logit link #----------------------------------------------------------------------------------- library(MASS) DataSet3.2$trt <- factor( x = DataSet3.2$trt ) DataSet3.2$loc <- factor( x = DataSet3.2$loc ) Exam3.9.fm1 <- glmmPQL( fixed = S2/Nbin~trt , random = ~1|loc , family = quasibinomial(link = "logit") , data = DataSet3.2 , niter = 10 , verbose = TRUE ) summary(Exam3.9.fm1) library(parameters) model_parameters(Exam3.9.fm1) #------------------------------------------------------------- ## treatment means #------------------------------------------------------------- library(emmeans) emmeans(object = Exam3.9.fm1, specs = ~trt, type = "response") emmeans(object = Exam3.9.fm1, specs = ~trt, type = "link") emmeans(object = Exam3.9.fm1, specs = ~trt, type = "logit") ##--- Normal Approximation library(nlme) Exam3.9fm2 <- lme( fixed = S2/Nbin~trt , data = DataSet3.2 , random = ~1|loc , method = c("REML", "ML")[1] ) Exam3.9fm2 model_parameters(Exam3.9fm2) emmeans(object = Exam3.9fm2, specs = ~trt) ##---Binomial GLMM with interaction Exam3.9fm3 <- glmmPQL( fixed = S2/Nbin~trt , random = ~1|trt/loc , family = quasibinomial(link = "logit") , data = DataSet3.2 , niter = 10 , verbose = TRUE ) summary(Exam3.9fm3) model_parameters(Exam3.9fm3) emmeans(object = Exam3.9fm3, specs = ~trt) ##---Binomial Marginal GLMM(assuming compound symmetry) Exam3.9fm4 <- glmmPQL( fixed = S2/Nbin~trt , random = ~1|loc , family = quasibinomial(link = "logit") , data = DataSet3.2 , correlation = corCompSymm(form = ~1|loc) , niter = 10 , verbose = TRUE ) summary(Exam3.9fm4) model_parameters(Exam3.9fm4) emmeans(object = Exam3.9fm4, specs = ~trt)#----------------------------------------------------------------------------------- ## Binomial conditional GLMM without interaction, logit link #----------------------------------------------------------------------------------- library(MASS) DataSet3.2$trt <- factor( x = DataSet3.2$trt ) DataSet3.2$loc <- factor( x = DataSet3.2$loc ) Exam3.9.fm1 <- glmmPQL( fixed = S2/Nbin~trt , random = ~1|loc , family = quasibinomial(link = "logit") , data = DataSet3.2 , niter = 10 , verbose = TRUE ) summary(Exam3.9.fm1) library(parameters) model_parameters(Exam3.9.fm1) #------------------------------------------------------------- ## treatment means #------------------------------------------------------------- library(emmeans) emmeans(object = Exam3.9.fm1, specs = ~trt, type = "response") emmeans(object = Exam3.9.fm1, specs = ~trt, type = "link") emmeans(object = Exam3.9.fm1, specs = ~trt, type = "logit") ##--- Normal Approximation library(nlme) Exam3.9fm2 <- lme( fixed = S2/Nbin~trt , data = DataSet3.2 , random = ~1|loc , method = c("REML", "ML")[1] ) Exam3.9fm2 model_parameters(Exam3.9fm2) emmeans(object = Exam3.9fm2, specs = ~trt) ##---Binomial GLMM with interaction Exam3.9fm3 <- glmmPQL( fixed = S2/Nbin~trt , random = ~1|trt/loc , family = quasibinomial(link = "logit") , data = DataSet3.2 , niter = 10 , verbose = TRUE ) summary(Exam3.9fm3) model_parameters(Exam3.9fm3) emmeans(object = Exam3.9fm3, specs = ~trt) ##---Binomial Marginal GLMM(assuming compound symmetry) Exam3.9fm4 <- glmmPQL( fixed = S2/Nbin~trt , random = ~1|loc , family = quasibinomial(link = "logit") , data = DataSet3.2 , correlation = corCompSymm(form = ~1|loc) , niter = 10 , verbose = TRUE ) summary(Exam3.9fm4) model_parameters(Exam3.9fm4) emmeans(object = Exam3.9fm4, specs = ~trt)
Exam4.1 REML vs ML criterion is used keeping block effects random
Muhammad Yaseen ([email protected])
Adeela Munawar ([email protected])
Stroup, W. W. (2012). Generalized Linear Mixed Models: Modern Concepts, Methods and Applications. CRC Press.
DataSet4.1$trt <- factor(x = DataSet4.1$trt) DataSet4.1$block <- factor(x = DataSet4.1$block) #---REML estimates on page 138(article 4.4.3.3) library(lmerTest) Exam4.1REML <- lmer(formula = y~ trt +( 1|block ), data = DataSet4.1) library(parameters) model_parameters(Exam4.1REML) print(VarCorr(x = Exam4.1REML), comp = c("Variance")) ##---ML estimates on page 138(article 4.4.3.3) Exam4.1ML <- lmer(formula = y ~ trt + (1|block), data = DataSet4.1, REML = FALSE) model_parameters(Exam4.1ML) print(VarCorr(x = Exam4.1ML), comp = c("Variance")) Exam4.1.lm <- lm(formula = y~ trt + block, data = DataSet4.1) anova(object = Exam4.1.lm)DataSet4.1$trt <- factor(x = DataSet4.1$trt) DataSet4.1$block <- factor(x = DataSet4.1$block) #---REML estimates on page 138(article 4.4.3.3) library(lmerTest) Exam4.1REML <- lmer(formula = y~ trt +( 1|block ), data = DataSet4.1) library(parameters) model_parameters(Exam4.1REML) print(VarCorr(x = Exam4.1REML), comp = c("Variance")) ##---ML estimates on page 138(article 4.4.3.3) Exam4.1ML <- lmer(formula = y ~ trt + (1|block), data = DataSet4.1, REML = FALSE) model_parameters(Exam4.1ML) print(VarCorr(x = Exam4.1ML), comp = c("Variance")) Exam4.1.lm <- lm(formula = y~ trt + block, data = DataSet4.1) anova(object = Exam4.1.lm)
Exam5.1 is used to show polynomial multiple regression with binomial response
Muhammad Yaseen ([email protected])
Adeela Munawar ([email protected])
Stroup, W. W. (2012). Generalized Linear Mixed Models: Modern Concepts, Methods and Applications. CRC Press.
##---Sequential Fit of the logit Model Exam5.1.glm.1 <- glm( formula = F/N~ X , family = quasibinomial(link = "logit") , data = DataSet5.1 ) summary(Exam5.1.glm.1) library(parameters) model_parameters(Exam5.1.glm.1) ## confint.default() produce Wald Confidence interval as SAS produces ##---Likelihood Ratio test for Model 1 anova(object = Exam5.1.glm.1, test = "Chisq") library(aod) WaldExam5.1.glm.1 <- wald.test( Sigma = vcov(object = Exam5.1.glm.1) , b = coef(object = Exam5.1.glm.1) , Terms = 2 , L = NULL , H0 = NULL , df = NULL , verbose = FALSE ) ##---Sequential Fit of the logit Model quadratic terms involved Exam5.1.glm.2 <- glm( formula = F/N~ X + I(X^2) , family = quasibinomial(link = "logit") , data = DataSet5.1 ) summary( Exam5.1.glm.2 ) model_parameters( Exam5.1.glm.2 ) ##---Likelihood Ratio test for Model Exam5.1.glm.2 anova(object = Exam5.1.glm.2, test = "Chisq") WaldExam5.1.glm.2 <- wald.test( Sigma = vcov(object = Exam5.1.glm.2) , b = coef(object = Exam5.1.glm.2) , Terms = 3 , L = NULL , H0 = NULL , df = NULL , verbose = FALSE ) ##---Sequential Fit of the logit Model 5th power terms involved Exam5.1.glm.3 <- glm( formula = F/N~ X + I(X^2) + I(X^3) + I(X^4) + I(X^5) , family = quasibinomial(link = "logit") , data = DataSet5.1 ) summary(Exam5.1.glm.3) model_parameters(Exam5.1.glm.3) ## confint.default() produce Wald Confidence interval as SAS produces ##---Likelihood Ratio test for Model 1 anova(object = Exam5.1.glm.3, test = "Chisq") WaldExam5.1.glm.3 <- wald.test( Sigma = vcov(object = Exam5.1.glm.3) , b = coef(object = Exam5.1.glm.3) , Terms = 6 , L = NULL , H0 = NULL , df = NULL , verbose = FALSE )##---Sequential Fit of the logit Model Exam5.1.glm.1 <- glm( formula = F/N~ X , family = quasibinomial(link = "logit") , data = DataSet5.1 ) summary(Exam5.1.glm.1) library(parameters) model_parameters(Exam5.1.glm.1) ## confint.default() produce Wald Confidence interval as SAS produces ##---Likelihood Ratio test for Model 1 anova(object = Exam5.1.glm.1, test = "Chisq") library(aod) WaldExam5.1.glm.1 <- wald.test( Sigma = vcov(object = Exam5.1.glm.1) , b = coef(object = Exam5.1.glm.1) , Terms = 2 , L = NULL , H0 = NULL , df = NULL , verbose = FALSE ) ##---Sequential Fit of the logit Model quadratic terms involved Exam5.1.glm.2 <- glm( formula = F/N~ X + I(X^2) , family = quasibinomial(link = "logit") , data = DataSet5.1 ) summary( Exam5.1.glm.2 ) model_parameters( Exam5.1.glm.2 ) ##---Likelihood Ratio test for Model Exam5.1.glm.2 anova(object = Exam5.1.glm.2, test = "Chisq") WaldExam5.1.glm.2 <- wald.test( Sigma = vcov(object = Exam5.1.glm.2) , b = coef(object = Exam5.1.glm.2) , Terms = 3 , L = NULL , H0 = NULL , df = NULL , verbose = FALSE ) ##---Sequential Fit of the logit Model 5th power terms involved Exam5.1.glm.3 <- glm( formula = F/N~ X + I(X^2) + I(X^3) + I(X^4) + I(X^5) , family = quasibinomial(link = "logit") , data = DataSet5.1 ) summary(Exam5.1.glm.3) model_parameters(Exam5.1.glm.3) ## confint.default() produce Wald Confidence interval as SAS produces ##---Likelihood Ratio test for Model 1 anova(object = Exam5.1.glm.3, test = "Chisq") WaldExam5.1.glm.3 <- wald.test( Sigma = vcov(object = Exam5.1.glm.3) , b = coef(object = Exam5.1.glm.3) , Terms = 6 , L = NULL , H0 = NULL , df = NULL , verbose = FALSE )
Exam5.2 three factor main effects only design
Muhammad Yaseen ([email protected])
Adeela Munawar ([email protected])
Stroup, W. W. (2012). Generalized Linear Mixed Models: Modern Concepts, Methods and Applications. CRC Press.
DataSet5.2$a <- factor( x = DataSet5.2$a) DataSet5.2$b <- factor( x = DataSet5.2$b) DataSet5.2$c <- factor(x = DataSet5.2$c) ##---first adding factor a in model Exam5.2.lm1 <- lm(formula = y~ a, data = DataSet5.2) summary(Exam5.2.lm1) library(parameters) model_parameters(Exam5.2.lm1) library(emmeans) ##---A first emmeans(object = Exam5.2.lm1, specs = ~a) contrast(emmeans(object = Exam5.2.lm1, specs = ~a), method = "pairwise") anova(object = Exam5.2.lm1) ##---then adding factor b in model Exam5.2.lm2 <- lm(formula = y~ a + b, data = DataSet5.2) summary(Exam5.2.lm2) model_parameters(Exam5.2.lm2) emmeans(object = Exam5.2.lm2, specs = ~b) contrast(emmeans(object = Exam5.2.lm2, specs = ~b), method = "pairwise") anova(object = Exam5.2.lm2) ##---then adding factor c in model Exam5.2.lm3 <- lm(formula = y~ a + b + c, data = DataSet5.2) summary(Exam5.2.lm3) model_parameters(Exam5.2.lm3) emmeans(object = Exam5.2.lm3, specs = ~c) contrast(emmeans(object = Exam5.2.lm3, specs = ~c), method = "pairwise") anova(object = Exam5.2.lm3) ##---Now Change the order and add b first in model Exam5.2.lm4 <- lm(formula = y ~ b, data = DataSet5.2) summary(Exam5.2.lm4) model_parameters(Exam5.2.lm4) emmeans(object = Exam5.2.lm4, specs = ~b) contrast(emmeans(object = Exam5.2.lm4, specs = ~b), method = "pairwise") anova(object = Exam5.2.lm4) ##---then adding factor a in model Exam5.2.lm5 <- lm(formula = y ~ b + a, data = DataSet5.2) summary(Exam5.2.lm5) model_parameters(Exam5.2.lm5) emmeans(object = Exam5.2.lm5, specs = ~a) contrast(emmeans(object = Exam5.2.lm5, specs = ~a), method = "pairwise") anova(object = Exam5.2.lm5)DataSet5.2$a <- factor( x = DataSet5.2$a) DataSet5.2$b <- factor( x = DataSet5.2$b) DataSet5.2$c <- factor(x = DataSet5.2$c) ##---first adding factor a in model Exam5.2.lm1 <- lm(formula = y~ a, data = DataSet5.2) summary(Exam5.2.lm1) library(parameters) model_parameters(Exam5.2.lm1) library(emmeans) ##---A first emmeans(object = Exam5.2.lm1, specs = ~a) contrast(emmeans(object = Exam5.2.lm1, specs = ~a), method = "pairwise") anova(object = Exam5.2.lm1) ##---then adding factor b in model Exam5.2.lm2 <- lm(formula = y~ a + b, data = DataSet5.2) summary(Exam5.2.lm2) model_parameters(Exam5.2.lm2) emmeans(object = Exam5.2.lm2, specs = ~b) contrast(emmeans(object = Exam5.2.lm2, specs = ~b), method = "pairwise") anova(object = Exam5.2.lm2) ##---then adding factor c in model Exam5.2.lm3 <- lm(formula = y~ a + b + c, data = DataSet5.2) summary(Exam5.2.lm3) model_parameters(Exam5.2.lm3) emmeans(object = Exam5.2.lm3, specs = ~c) contrast(emmeans(object = Exam5.2.lm3, specs = ~c), method = "pairwise") anova(object = Exam5.2.lm3) ##---Now Change the order and add b first in model Exam5.2.lm4 <- lm(formula = y ~ b, data = DataSet5.2) summary(Exam5.2.lm4) model_parameters(Exam5.2.lm4) emmeans(object = Exam5.2.lm4, specs = ~b) contrast(emmeans(object = Exam5.2.lm4, specs = ~b), method = "pairwise") anova(object = Exam5.2.lm4) ##---then adding factor a in model Exam5.2.lm5 <- lm(formula = y ~ b + a, data = DataSet5.2) summary(Exam5.2.lm5) model_parameters(Exam5.2.lm5) emmeans(object = Exam5.2.lm5, specs = ~a) contrast(emmeans(object = Exam5.2.lm5, specs = ~a), method = "pairwise") anova(object = Exam5.2.lm5)
Exam5.3 Inference using empirical standard error with different Bias connection
Muhammad Yaseen ([email protected])
Adeela Munawar ([email protected])
Stroup, W. W. (2012). Generalized Linear Mixed Models: Modern Concepts, Methods and Applications. CRC Press.
data(DataSet4.1) DataSet4.1$trt <- factor(x = DataSet4.1$trt) DataSet4.1$block <- factor( x = DataSet4.1$block) ##---REML estimates on page 172 library(lmerTest) Exam5.3REML <- lmerTest::lmer(formula = y ~ trt + (1|block), data = DataSet4.1, REML = TRUE) Exam5.3REML library(parameters) model_parameters(Exam5.3REML) ##---Standard Error Type "Model Based" with no Bias Connection anova(object = Exam5.3REML) anova(object = Exam5.3REML, ddf = "Satterthwaite") ##---Standard Error Type "Model Based" with "Kenward-Roger approximation" Bias Connection anova(object = Exam5.3REML, ddf = "Kenward-Roger") ##---ML estimates on page 172 Exam5.3ML <- lmerTest::lmer(formula = y ~ trt + ( 1|block ), data = DataSet4.1, REML = FALSE) Exam5.3ML library(parameters) model_parameters(Exam5.3ML) ##---Standard Error Type "Model Based" with no Bias Connection anova(object = Exam5.3ML ) anova(object = Exam5.3ML, ddf = "Satterthwaite")data(DataSet4.1) DataSet4.1$trt <- factor(x = DataSet4.1$trt) DataSet4.1$block <- factor( x = DataSet4.1$block) ##---REML estimates on page 172 library(lmerTest) Exam5.3REML <- lmerTest::lmer(formula = y ~ trt + (1|block), data = DataSet4.1, REML = TRUE) Exam5.3REML library(parameters) model_parameters(Exam5.3REML) ##---Standard Error Type "Model Based" with no Bias Connection anova(object = Exam5.3REML) anova(object = Exam5.3REML, ddf = "Satterthwaite") ##---Standard Error Type "Model Based" with "Kenward-Roger approximation" Bias Connection anova(object = Exam5.3REML, ddf = "Kenward-Roger") ##---ML estimates on page 172 Exam5.3ML <- lmerTest::lmer(formula = y ~ trt + ( 1|block ), data = DataSet4.1, REML = FALSE) Exam5.3ML library(parameters) model_parameters(Exam5.3ML) ##---Standard Error Type "Model Based" with no Bias Connection anova(object = Exam5.3ML ) anova(object = Exam5.3ML, ddf = "Satterthwaite")
Exam7.1 explains multifactor models with all factors qualitative
Muhammad Yaseen ([email protected])
Adeela Munawar ([email protected])
Stroup, W. W. (2012). Generalized Linear Mixed Models: Modern Concepts, Methods and Applications. CRC Press.
@seealso
DataSet7.1
library(emmeans) library(car) data(DataSet7.1) DataSet7.1$a <- factor(x = DataSet7.1$a) DataSet7.1$b <- factor(x = DataSet7.1$b) Exam7.1.lm1 <- lm(formula = y ~ a + b + a*b, data = DataSet7.1) summary(Exam7.1.lm1) library(parameters) model_parameters(Exam7.1.lm1) anova(Exam7.1.lm1) ##---Result obtained as in SLICE statement in SAS for a0 & a1 library(phia) testInteractions( model = Exam7.1.lm1 , custom = list(a = c("0" = 1)) , across = "b" ) testInteractions( model = Exam7.1.lm1 , custom = list(a = c("1" = 1)) , across = "b" ) ##---Interaction plot emmip( object = Exam7.1.lm1 , formula = a~b , ylab = "y Lsmeans" , main = "Lsmeans for a*b" ) #------------------------------------------------------------- ## Individula least squares treatment means #------------------------------------------------------------- emmeans(object = Exam7.1.lm1, specs = ~a*b) ##---Simpe effects comparison of interaction by a ## (but it doesn't give the same p-value as in article 7.4.2 page#215) emmeans(object = Exam7.1.lm1, specs = pairwise~b|a)$contrasts pairs(emmeans(object = Exam7.1.lm1, specs = ~b|a), simple = "each", combine = TRUE) pairs(emmeans(object = Exam7.1.lm1, specs = ~b|a), simple = "a") pairs(emmeans(object = Exam7.1.lm1, specs = ~b|a), simple = "b") pairs(emmeans(object = Exam7.1.lm1, specs = ~b|a)) contrast(emmeans(object = Exam7.1.lm1, specs = ~b|a)) emmeans(object = Exam7.1.lm1, specs = pairwise~b|a) emmeans(object = Exam7.1.lm1, specs = pairwise~b|a)$contrasts ##---Alternative method of pairwise comparisons by ## applying contrast ## coefficient (gives the same p-value as in 7.4.2) contrast( emmeans(object = Exam7.1.lm1, specs = ~a*b) , list ( c1 = c(1, 0, -1, 0, 0, 0) , c2 = c(1, 0, 0, 0, -1, 0) , c3 = c(0, 0, 1, 0, -1, 0) , c4 = c(0, 1, 0, -1, 0, 0) , c5 = c(0, 1, 0, 0, 0, -1) , c6 = c(0, 1, 0, 0, -1, 0) ) ) ##---Nested Model (page 216)---- Exam7.1.lm2 <- lm(formula = y ~ a + a %in% b, data = DataSet7.1) summary(Exam7.1.lm2) model_parameters(Exam7.1.lm2) anova(Exam7.1.lm2) car::linearHypothesis(Exam7.1.lm2, c("a0:b1 = a0:b2")) car::linearHypothesis(Exam7.1.lm2, c("a1:b1 = a1:b2")) ##---Bonferroni's adjusted p-values emmeans(object = Exam7.1.lm2, specs = pairwise~b|a, adjust = "bonferroni")$contrasts ##--- Alternative method of pairwise comparisons by ## applying contrast coefficient with Bonferroni adjustment contrast( emmeans(object = Exam7.1.lm1, specs = ~a*b) , list ( c1 = c(1, 0, -1, 0, 0, 0) , c2 = c(1, 0, 0, 0, -1, 0) , c3 = c(0, 0, 1, 0, -1, 0) , c4 = c(0, 1, 0, -1, 0, 0) , c5 = c(0, 1, 0, 0, 0, -1) , c6 = c(0, 1, 0, 0, -1, 0) ) , adjust = "bonferroni" )library(emmeans) library(car) data(DataSet7.1) DataSet7.1$a <- factor(x = DataSet7.1$a) DataSet7.1$b <- factor(x = DataSet7.1$b) Exam7.1.lm1 <- lm(formula = y ~ a + b + a*b, data = DataSet7.1) summary(Exam7.1.lm1) library(parameters) model_parameters(Exam7.1.lm1) anova(Exam7.1.lm1) ##---Result obtained as in SLICE statement in SAS for a0 & a1 library(phia) testInteractions( model = Exam7.1.lm1 , custom = list(a = c("0" = 1)) , across = "b" ) testInteractions( model = Exam7.1.lm1 , custom = list(a = c("1" = 1)) , across = "b" ) ##---Interaction plot emmip( object = Exam7.1.lm1 , formula = a~b , ylab = "y Lsmeans" , main = "Lsmeans for a*b" ) #------------------------------------------------------------- ## Individula least squares treatment means #------------------------------------------------------------- emmeans(object = Exam7.1.lm1, specs = ~a*b) ##---Simpe effects comparison of interaction by a ## (but it doesn't give the same p-value as in article 7.4.2 page#215) emmeans(object = Exam7.1.lm1, specs = pairwise~b|a)$contrasts pairs(emmeans(object = Exam7.1.lm1, specs = ~b|a), simple = "each", combine = TRUE) pairs(emmeans(object = Exam7.1.lm1, specs = ~b|a), simple = "a") pairs(emmeans(object = Exam7.1.lm1, specs = ~b|a), simple = "b") pairs(emmeans(object = Exam7.1.lm1, specs = ~b|a)) contrast(emmeans(object = Exam7.1.lm1, specs = ~b|a)) emmeans(object = Exam7.1.lm1, specs = pairwise~b|a) emmeans(object = Exam7.1.lm1, specs = pairwise~b|a)$contrasts ##---Alternative method of pairwise comparisons by ## applying contrast ## coefficient (gives the same p-value as in 7.4.2) contrast( emmeans(object = Exam7.1.lm1, specs = ~a*b) , list ( c1 = c(1, 0, -1, 0, 0, 0) , c2 = c(1, 0, 0, 0, -1, 0) , c3 = c(0, 0, 1, 0, -1, 0) , c4 = c(0, 1, 0, -1, 0, 0) , c5 = c(0, 1, 0, 0, 0, -1) , c6 = c(0, 1, 0, 0, -1, 0) ) ) ##---Nested Model (page 216)---- Exam7.1.lm2 <- lm(formula = y ~ a + a %in% b, data = DataSet7.1) summary(Exam7.1.lm2) model_parameters(Exam7.1.lm2) anova(Exam7.1.lm2) car::linearHypothesis(Exam7.1.lm2, c("a0:b1 = a0:b2")) car::linearHypothesis(Exam7.1.lm2, c("a1:b1 = a1:b2")) ##---Bonferroni's adjusted p-values emmeans(object = Exam7.1.lm2, specs = pairwise~b|a, adjust = "bonferroni")$contrasts ##--- Alternative method of pairwise comparisons by ## applying contrast coefficient with Bonferroni adjustment contrast( emmeans(object = Exam7.1.lm1, specs = ~a*b) , list ( c1 = c(1, 0, -1, 0, 0, 0) , c2 = c(1, 0, 0, 0, -1, 0) , c3 = c(0, 0, 1, 0, -1, 0) , c4 = c(0, 1, 0, -1, 0, 0) , c5 = c(0, 1, 0, 0, 0, -1) , c6 = c(0, 1, 0, 0, -1, 0) ) , adjust = "bonferroni" )
Exam7.2 explains multifactor models with some factors qualitative and some quantitative(Equal slopes ANCOVA)
Muhammad Yaseen ([email protected])
Adeela Munawar ([email protected])
Stroup, W. W. (2012). Generalized Linear Mixed Models: Modern Concepts, Methods and Applications. CRC Press.
@seealso
DataSet7.2
library(emmeans) library(car) library(ggplot2) data(DataSet7.2) DataSet7.2$trt <- factor( x = DataSet7.2$trt ) ##----ANCOVA(Equal slope Model) Exam7.2fm1 <- aov(formula = y ~ trt*x, data = DataSet7.2) car::Anova(mod = Exam7.2fm1 , type = "III") ##---ANCOVA(without interaction because of non significant slope effect) Exam7.2fm2 <- aov(formula = y ~ trt + x, data = DataSet7.2) car::Anova(mod = Exam7.2fm2 , type = "III") ##---Ls means for 2nd model emmeans(object = Exam7.2fm2, specs = ~trt) ##---Anova without covariate Exam7.2fm3 <- aov(formula = y ~ trt, data = DataSet7.2) car::Anova(mod = Exam7.2fm3, type = "III") ##---Ls means for 3rd model emmeans(object = Exam7.2fm3, specs = ~trt) ##---Box Plot of Covariate by treatment Plot <- ggplot( data = DataSet7.2 , mapping = aes(x = factor(trt), y = x) ) + geom_boxplot(width = 0.5) + coord_flip() + geom_point() + stat_summary( fun = "mean" , geom = "point" , shape = 23 , size = 2 , fill = "red" ) + theme_bw() + ggtitle("Covariate by treatment Box Plot") + xlab("Treatment") print(Plot)library(emmeans) library(car) library(ggplot2) data(DataSet7.2) DataSet7.2$trt <- factor( x = DataSet7.2$trt ) ##----ANCOVA(Equal slope Model) Exam7.2fm1 <- aov(formula = y ~ trt*x, data = DataSet7.2) car::Anova(mod = Exam7.2fm1 , type = "III") ##---ANCOVA(without interaction because of non significant slope effect) Exam7.2fm2 <- aov(formula = y ~ trt + x, data = DataSet7.2) car::Anova(mod = Exam7.2fm2 , type = "III") ##---Ls means for 2nd model emmeans(object = Exam7.2fm2, specs = ~trt) ##---Anova without covariate Exam7.2fm3 <- aov(formula = y ~ trt, data = DataSet7.2) car::Anova(mod = Exam7.2fm3, type = "III") ##---Ls means for 3rd model emmeans(object = Exam7.2fm3, specs = ~trt) ##---Box Plot of Covariate by treatment Plot <- ggplot( data = DataSet7.2 , mapping = aes(x = factor(trt), y = x) ) + geom_boxplot(width = 0.5) + coord_flip() + geom_point() + stat_summary( fun = "mean" , geom = "point" , shape = 23 , size = 2 , fill = "red" ) + theme_bw() + ggtitle("Covariate by treatment Box Plot") + xlab("Treatment") print(Plot)
Exam7.3 explains multifactor models with some factors qualitative and some quantitative(Unequal slopes ANCOVA)
Muhammad Yaseen ([email protected])
Adeela Munawar ([email protected])
Stroup, W. W. (2012). Generalized Linear Mixed Models: Modern Concepts, Methods and Applications. CRC Press.
@seealso
DataSet7.3
library(car) library(ggplot2) library(emmeans) data(DataSet7.3) DataSet7.3$trt <- factor(x = DataSet7.3$trt ) ##----ANCOVA(Unequal slope Model) Exam7.3fm1 <- aov(formula = y ~ trt*x, data = DataSet7.3) car::Anova( mod = Exam7.3fm1 , type = "III") Plot <- ggplot( data = DataSet7.3 , mapping = aes(x = factor(trt), y = x) ) + geom_boxplot(width = 0.5) + coord_flip() + geom_point() + stat_summary( fun = "mean" , geom = "point" , shape = 23 , size = 2 , fill = "red" ) + theme_bw() + ggtitle("Covariate by treatment Box Plot") + xlab("Treatment") print(Plot) ##----ANCOVA(Unequal slope Model without intercept at page 224) Exam7.3fm2 <- lm(formula = y ~ 0 + trt/x, data = DataSet7.3) summary(Exam7.3fm2) library(parameters) model_parameters(Exam7.3fm2) ##--Lsmeans treatment at x=7 & 12 at page 225 emmeans(object = Exam7.3fm2, specs = ~trt|x, at = list(x = c(7, 12)))library(car) library(ggplot2) library(emmeans) data(DataSet7.3) DataSet7.3$trt <- factor(x = DataSet7.3$trt ) ##----ANCOVA(Unequal slope Model) Exam7.3fm1 <- aov(formula = y ~ trt*x, data = DataSet7.3) car::Anova( mod = Exam7.3fm1 , type = "III") Plot <- ggplot( data = DataSet7.3 , mapping = aes(x = factor(trt), y = x) ) + geom_boxplot(width = 0.5) + coord_flip() + geom_point() + stat_summary( fun = "mean" , geom = "point" , shape = 23 , size = 2 , fill = "red" ) + theme_bw() + ggtitle("Covariate by treatment Box Plot") + xlab("Treatment") print(Plot) ##----ANCOVA(Unequal slope Model without intercept at page 224) Exam7.3fm2 <- lm(formula = y ~ 0 + trt/x, data = DataSet7.3) summary(Exam7.3fm2) library(parameters) model_parameters(Exam7.3fm2) ##--Lsmeans treatment at x=7 & 12 at page 225 emmeans(object = Exam7.3fm2, specs = ~trt|x, at = list(x = c(7, 12)))
Exam7.6.2.1 Nonlinear Mean Models ( Quantitative by quantitative models)
Muhammad Yaseen ([email protected])
Adeela Munawar ([email protected])
Stroup, W. W. (2012). Generalized Linear Mixed Models: Modern Concepts, Methods and Applications. CRC Press.
@seealso
DataSet7.6
library(scatterplot3d) data(DataSet7.6) library(dplyr) library(magrittr) DataSet7.6 <- DataSet7.6 %>% mutate( logx1 = ifelse(test = x1 == 0, yes = log(x1 + 0.1), no = log(x1)) , logx2 = ifelse(test = x2 == 0, yes = log(x2 + 0.1), no = log(x2)) ) DataSet7.6 Exam7.6.2.1.lm <- lm(formula = response ~ x1*x2 + logx1*logx2 , data = DataSet7.6) summary(Exam7.6.2.1.lm) library(parameters) model_parameters(Exam7.6.2.1.lm) ##---3D Scatter plot ( page#232) attach(DataSet7.6) ( ScatterPlot1 <- scatterplot3d( x = x1 , y = x2 , z = response , color = response , main = " 3D Scatter plot of response") ) ##--- scatter plot with regression plane by using Hoerl function ( page#233) grid.lines <- 5 x1.pred <- seq(min(x1), max(x1), length.out = grid.lines) x2.pred <- seq(min(x2), max(x2), length.out = grid.lines) x1x2 <- expand.grid( x = x1.pred, y = x2.pred) z.pred <- matrix(data = predict(Exam7.6.2.1.lm, newdata = x1x2), nrow = grid.lines , ncol = grid.lines) (ScatterPlot2 <- scatterplot3d( x = x1 , y = x2 , z = response , pch = 20 , phi = 25 , theta = 30 , ticktype = "detailed" , xlab = "x1" , ylab = "x2" , zlab = "response" , add = FALSE , surf = list(x = x1.pred , y = x2.pred , z = z.pred , facets = NA ) , plot = TRUE , main = "Fitted Response Surface by Hoerl Function" ) )library(scatterplot3d) data(DataSet7.6) library(dplyr) library(magrittr) DataSet7.6 <- DataSet7.6 %>% mutate( logx1 = ifelse(test = x1 == 0, yes = log(x1 + 0.1), no = log(x1)) , logx2 = ifelse(test = x2 == 0, yes = log(x2 + 0.1), no = log(x2)) ) DataSet7.6 Exam7.6.2.1.lm <- lm(formula = response ~ x1*x2 + logx1*logx2 , data = DataSet7.6) summary(Exam7.6.2.1.lm) library(parameters) model_parameters(Exam7.6.2.1.lm) ##---3D Scatter plot ( page#232) attach(DataSet7.6) ( ScatterPlot1 <- scatterplot3d( x = x1 , y = x2 , z = response , color = response , main = " 3D Scatter plot of response") ) ##--- scatter plot with regression plane by using Hoerl function ( page#233) grid.lines <- 5 x1.pred <- seq(min(x1), max(x1), length.out = grid.lines) x2.pred <- seq(min(x2), max(x2), length.out = grid.lines) x1x2 <- expand.grid( x = x1.pred, y = x2.pred) z.pred <- matrix(data = predict(Exam7.6.2.1.lm, newdata = x1x2), nrow = grid.lines , ncol = grid.lines) (ScatterPlot2 <- scatterplot3d( x = x1 , y = x2 , z = response , pch = 20 , phi = 25 , theta = 30 , ticktype = "detailed" , xlab = "x1" , ylab = "x2" , zlab = "response" , add = FALSE , surf = list(x = x1.pred , y = x2.pred , z = z.pred , facets = NA ) , plot = TRUE , main = "Fitted Response Surface by Hoerl Function" ) )
Exam7.7 is an explaination of segmented regression
Muhammad Yaseen ([email protected])
Adeela Munawar ([email protected])
Stroup, W. W. (2012). Generalized Linear Mixed Models: Modern Concepts, Methods and Applications. CRC Press.
library(splines) library(ggplot2) DataSet7.7$a <- factor(x = DataSet7.7$a) knots <- c(0, 0, 0, 0, 10, 10, 20, 30, 40, 40, 40, 45, 45, 45, 50, 50, 50) bx <- splineDesign(knots = knots, x = DataSet7.7$x, outer.ok = TRUE) Exam7.7fm <- lm(formula = y ~ a*bx, data = DataSet7.7) anova(Exam7.7fm) Data <- data.frame(DataSet7.7, fit = Exam7.7fm$fit) ##---Estimated response surface by using segmented regression Plot <- ggplot(data = Data , mapping = aes(x = x, y = y, colour = a)) + geom_point() + geom_line(linewidth = 1) + ggtitle("Response surface by using segmented regression") print(Plot)library(splines) library(ggplot2) DataSet7.7$a <- factor(x = DataSet7.7$a) knots <- c(0, 0, 0, 0, 10, 10, 20, 30, 40, 40, 40, 45, 45, 45, 50, 50, 50) bx <- splineDesign(knots = knots, x = DataSet7.7$x, outer.ok = TRUE) Exam7.7fm <- lm(formula = y ~ a*bx, data = DataSet7.7) anova(Exam7.7fm) Data <- data.frame(DataSet7.7, fit = Exam7.7fm$fit) ##---Estimated response surface by using segmented regression Plot <- ggplot(data = Data , mapping = aes(x = x, y = y, colour = a)) + geom_point() + geom_line(linewidth = 1) + ggtitle("Response surface by using segmented regression") print(Plot)
Exam8.1 Nested factorial structure
Muhammad Yaseen ([email protected])
Adeela Munawar ([email protected])
Stroup, W. W. (2012). Generalized Linear Mixed Models: Modern Concepts, Methods and Applications. CRC Press.
data(DataSet8.1) DataSet8.1$block <- factor(x = DataSet8.1$block) DataSet8.1$set <- factor(x = DataSet8.1$set) DataSet8.1$trt <- factor(x = DataSet8.1$trt) library(lmerTest) Exam8.1Lmer <- lmer(y ~ set + trt %in% set + (1|set/block), DataSet8.1) summary(Exam8.1Lmer) anova(Exam8.1Lmer) library(emmeans) emmeans(object = Exam8.1Lmer, specs = ~trt|set) contrast(emmeans(object = Exam8.1Lmer, specs = ~trt|set), method = "pairwise", by = "set")data(DataSet8.1) DataSet8.1$block <- factor(x = DataSet8.1$block) DataSet8.1$set <- factor(x = DataSet8.1$set) DataSet8.1$trt <- factor(x = DataSet8.1$trt) library(lmerTest) Exam8.1Lmer <- lmer(y ~ set + trt %in% set + (1|set/block), DataSet8.1) summary(Exam8.1Lmer) anova(Exam8.1Lmer) library(emmeans) emmeans(object = Exam8.1Lmer, specs = ~trt|set) contrast(emmeans(object = Exam8.1Lmer, specs = ~trt|set), method = "pairwise", by = "set")
Exam8.2 Incomplete strip-plot
Muhammad Yaseen ([email protected])
Adeela Munawar ([email protected])
Stroup, W. W. (2012). Generalized Linear Mixed Models: Modern Concepts, Methods and Applications. CRC Press.
data(DataSet8.2) DataSet8.2$block <- factor(x = DataSet8.2$block) DataSet8.2$a <- factor(x = DataSet8.2$a) DataSet8.2$b <- factor(x = DataSet8.2$b) library(lmerTest) Exam8.2lmer <- lmer( formula = y ~ a*b + (1|block) + (1|block:a) + (1|block:b) , data = DataSet8.2 ) anova(Exam8.2lmer,ddf="Kenward-Roger") library(emmeans) emmeans(object = Exam8.2lmer, specs = ~a|b) emmip( object = emmeans(object = Exam8.2lmer, specs = ~a|b) , formula = a~b , ylab = "y Lsmeans" , main = "Lsmeans for a*b" ) ##---Simple effect comparisons of a*b Least Squares Means by a ( page # 254) emmeans(Exam8.2lmer, pairwise ~ b|a)data(DataSet8.2) DataSet8.2$block <- factor(x = DataSet8.2$block) DataSet8.2$a <- factor(x = DataSet8.2$a) DataSet8.2$b <- factor(x = DataSet8.2$b) library(lmerTest) Exam8.2lmer <- lmer( formula = y ~ a*b + (1|block) + (1|block:a) + (1|block:b) , data = DataSet8.2 ) anova(Exam8.2lmer,ddf="Kenward-Roger") library(emmeans) emmeans(object = Exam8.2lmer, specs = ~a|b) emmip( object = emmeans(object = Exam8.2lmer, specs = ~a|b) , formula = a~b , ylab = "y Lsmeans" , main = "Lsmeans for a*b" ) ##---Simple effect comparisons of a*b Least Squares Means by a ( page # 254) emmeans(Exam8.2lmer, pairwise ~ b|a)
Exam8.3 explains Response surface design with incomplete blocking
Muhammad Yaseen ([email protected])
Adeela Munawar ([email protected])
Stroup, W. W. (2012). Generalized Linear Mixed Models: Modern Concepts, Methods and Applications. CRC Press.
## Response Surface Design with incomplete blocking (page 255) data(DataSet8.3) DataSet8.3$block <- factor(x = DataSet8.3$block) DataSet8.3$aa <- factor(x = DataSet8.3$a) DataSet8.3$bb <- factor(x = DataSet8.3$b) DataSet8.3$cc <- factor(x = DataSet8.3$c) library(lmerTest) library(lattice) Exam8.3.fm1 <- lmer( y ~ aa:bb:cc + a + b + c + I(a^2) + I(b^2) + I(c^2) + a*b + a*c + b*c + (1|block) , data = DataSet8.3 ) ##--- page 256 anova(Exam8.3.fm1, ddf = "Kenward-Roger", type = 1) Exam8.3.fm2 <- lmer( y ~ a + b + c + I(a^2) + I(b^2) + I(c^2) + a*b + a*c + b*c + (1|block) , data = DataSet8.3 ) ##--- page 257 anova(Exam8.3.fm2, ddf = "Kenward-Roger", type = 1) ##--- page 257 Exam8.3.fm3 <- lmer( y ~ a + b + c + I(a^2) + I(b^2) + a*c + b*c + (1|block) , DataSet8.3 ) anova(Exam8.3.fm3, ddf = "Kenward-Roger", type = 1) ##--- scatter plot with regression plane by using Hoerl function ( page#233) a <- seq(from = -1, to = 1, by = 1) b <- seq(from = -1, to = 1, by = 1) c <- seq(from = -1, to = 1, by = 1) abc <- expand.grid(a = a, b = b, c = c) Yhat <- NULL for(i in 1:nrow(abc)) { Yhat[i] <- 50.08500 + 1.6*abc$a[i] + 1.69375*abc$b[i] + 0.51875*abc$c[i]- 3.30250*I((abc$a[i])^2)-3.51500*I((abc$b)^2)[i] - 0.52500*(abc$a)[i]*(abc$c)[i]-1.16250*(abc$b)[i]*(abc$c)[i] } Newdata <- data.frame(abc, Yhat) Plot1 <- wireframe(Yhat ~ b*a, data = subset(Newdata,c==-1), xlab = "b", ylab = "a", main = "Predicte response surface at C=-1", colorkey = FALSE, drape = TRUE, scales = list(arrows = FALSE),xlim=c(max(b),(min(b))), screen = list(z = -50, x =-70) ) Plot2 <- wireframe(Yhat ~ b*a, data = subset(Newdata,c==0), xlab = "b", ylab = "a", main = "Predicte response surface at C=0", colorkey = FALSE , drape = TRUE, scales = list(arrows = FALSE),xlim=c(max(b),(min(b))), screen = list(z = -50, x =-70) ) Plot3 <- wireframe(Yhat ~ b*a, data = subset(Newdata,c==1), xlab = "b", ylab = "a", main = "Predicte response surface at C=1", colorkey = FALSE, drape = TRUE, scales = list(arrows = FALSE),xlim=c(max(b),(min(b))), screen = list(z = -50, x =-70) ) print(Plot1) print(Plot2) print(Plot3)## Response Surface Design with incomplete blocking (page 255) data(DataSet8.3) DataSet8.3$block <- factor(x = DataSet8.3$block) DataSet8.3$aa <- factor(x = DataSet8.3$a) DataSet8.3$bb <- factor(x = DataSet8.3$b) DataSet8.3$cc <- factor(x = DataSet8.3$c) library(lmerTest) library(lattice) Exam8.3.fm1 <- lmer( y ~ aa:bb:cc + a + b + c + I(a^2) + I(b^2) + I(c^2) + a*b + a*c + b*c + (1|block) , data = DataSet8.3 ) ##--- page 256 anova(Exam8.3.fm1, ddf = "Kenward-Roger", type = 1) Exam8.3.fm2 <- lmer( y ~ a + b + c + I(a^2) + I(b^2) + I(c^2) + a*b + a*c + b*c + (1|block) , data = DataSet8.3 ) ##--- page 257 anova(Exam8.3.fm2, ddf = "Kenward-Roger", type = 1) ##--- page 257 Exam8.3.fm3 <- lmer( y ~ a + b + c + I(a^2) + I(b^2) + a*c + b*c + (1|block) , DataSet8.3 ) anova(Exam8.3.fm3, ddf = "Kenward-Roger", type = 1) ##--- scatter plot with regression plane by using Hoerl function ( page#233) a <- seq(from = -1, to = 1, by = 1) b <- seq(from = -1, to = 1, by = 1) c <- seq(from = -1, to = 1, by = 1) abc <- expand.grid(a = a, b = b, c = c) Yhat <- NULL for(i in 1:nrow(abc)) { Yhat[i] <- 50.08500 + 1.6*abc$a[i] + 1.69375*abc$b[i] + 0.51875*abc$c[i]- 3.30250*I((abc$a[i])^2)-3.51500*I((abc$b)^2)[i] - 0.52500*(abc$a)[i]*(abc$c)[i]-1.16250*(abc$b)[i]*(abc$c)[i] } Newdata <- data.frame(abc, Yhat) Plot1 <- wireframe(Yhat ~ b*a, data = subset(Newdata,c==-1), xlab = "b", ylab = "a", main = "Predicte response surface at C=-1", colorkey = FALSE, drape = TRUE, scales = list(arrows = FALSE),xlim=c(max(b),(min(b))), screen = list(z = -50, x =-70) ) Plot2 <- wireframe(Yhat ~ b*a, data = subset(Newdata,c==0), xlab = "b", ylab = "a", main = "Predicte response surface at C=0", colorkey = FALSE , drape = TRUE, scales = list(arrows = FALSE),xlim=c(max(b),(min(b))), screen = list(z = -50, x =-70) ) Plot3 <- wireframe(Yhat ~ b*a, data = subset(Newdata,c==1), xlab = "b", ylab = "a", main = "Predicte response surface at C=1", colorkey = FALSE, drape = TRUE, scales = list(arrows = FALSE),xlim=c(max(b),(min(b))), screen = list(z = -50, x =-70) ) print(Plot1) print(Plot2) print(Plot3)
Exam8.4 Multifactor treatment and Multilevel design structures
Muhammad Yaseen ([email protected])
Adeela Munawar ([email protected])
Stroup, W. W. (2012). Generalized Linear Mixed Models: Modern Concepts, Methods and Applications. CRC Press.
data(DataSet8.4) DataSet8.4$block <- factor(x = DataSet8.4$block) DataSet8.4$a <- factor(x = DataSet8.4$a) DataSet8.4$b <- factor(x = DataSet8.4$b) library(lmerTest) Exam8.4lmer <- lmer( y ~ a + b %in% a + (1|block) + (1|block:a) + (1|block:b) , data = DataSet8.4 ) anova(Exam8.4lmer, ddf = "Kenward-Roger") library(emmeans) emmeans(object = Exam8.4lmer, specs = ~a|b)data(DataSet8.4) DataSet8.4$block <- factor(x = DataSet8.4$block) DataSet8.4$a <- factor(x = DataSet8.4$a) DataSet8.4$b <- factor(x = DataSet8.4$b) library(lmerTest) Exam8.4lmer <- lmer( y ~ a + b %in% a + (1|block) + (1|block:a) + (1|block:b) , data = DataSet8.4 ) anova(Exam8.4lmer, ddf = "Kenward-Roger") library(emmeans) emmeans(object = Exam8.4lmer, specs = ~a|b)
Exam9.1 One-way random effects only model
Muhammad Yaseen ([email protected])
Adeela Munawar ([email protected])
Stroup, W. W. (2012). Generalized Linear Mixed Models: Modern Concepts, Methods and Applications. CRC Press.
data(DataSet9.1) DataSet9.1$a <- factor(x = DataSet9.1$a) ##---Random effects model library(lmerTest) Exam9.1lmer <- lmer( y ~ 1 + (1|a), data = DataSet9.1) summary(Exam9.1lmer) ##---fixed effects model Exam9.1lmer2 <- lm(y ~ a, data = DataSet9.1) summary(Exam9.1lmer2) #--------------------------------------------------- ## Over all mean narrow( page # 274) #--------------------------------------------------- library(emmeans) library(phia) list9.1 <- list(a = c( "1" = 1/12,"2" = 1/12 , "3" = 1/12,"4" = 1/12 , "5" = 1/12,"6" = 1/12 , "7" = 1/12,"8" = 1/12 , "9" = 1/12,"10" = 1/12 , "11" = 1/12,"12" = 1/12 )) phia::testFactors(model = Exam9.1lmer2, levels = list9.1) #---BLUP Estimates (Table 9.1) coef <- unlist(ranef(Exam9.1lmer)) BLUPa <- NULL for( i in 1:length(coef)) { BLUPa[i] <- (mean(DataSet9.1$y)+coef[i]) } print(BLUPa)data(DataSet9.1) DataSet9.1$a <- factor(x = DataSet9.1$a) ##---Random effects model library(lmerTest) Exam9.1lmer <- lmer( y ~ 1 + (1|a), data = DataSet9.1) summary(Exam9.1lmer) ##---fixed effects model Exam9.1lmer2 <- lm(y ~ a, data = DataSet9.1) summary(Exam9.1lmer2) #--------------------------------------------------- ## Over all mean narrow( page # 274) #--------------------------------------------------- library(emmeans) library(phia) list9.1 <- list(a = c( "1" = 1/12,"2" = 1/12 , "3" = 1/12,"4" = 1/12 , "5" = 1/12,"6" = 1/12 , "7" = 1/12,"8" = 1/12 , "9" = 1/12,"10" = 1/12 , "11" = 1/12,"12" = 1/12 )) phia::testFactors(model = Exam9.1lmer2, levels = list9.1) #---BLUP Estimates (Table 9.1) coef <- unlist(ranef(Exam9.1lmer)) BLUPa <- NULL for( i in 1:length(coef)) { BLUPa[i] <- (mean(DataSet9.1$y)+coef[i]) } print(BLUPa)
Exam9.2 Two way random effects nested model
Muhammad Yaseen ([email protected])
Adeela Munawar ([email protected])
Stroup, W. W. (2012). Generalized Linear Mixed Models: Modern Concepts, Methods and Applications. CRC Press.
data(DataSet9.2) DataSet9.2$a <- factor(x = DataSet9.2$a) DataSet9.2$b <- factor(x = DataSet9.2$b) library(lmerTest) Exam9.2lmer <- lmer(y ~ (1|b/a), data = DataSet9.2) summary(Exam9.2lmer) Exam9.2lmer2 <- lm(y ~ a + b %in% a, data = DataSet9.2) summary(Exam9.2lmer2) ##--- Over all mean library(phia) list9.2 <- list(a = c("1" = 1/7,"2" = 1/7 , "3" = 1/7,"4" = 1/7 , "5" = 1/7,"6" = 1/7 , "7" = 1/7 )) phia::testFactors(model = Exam9.2lmer2, levels = list9.2) #---BLUP Estimates coef <- unlist(ranef(Exam9.2lmer)$a) BLUPa <- NULL for(i in 1:length(coef)){ BLUPa[i] <- (mean(DataSet9.2$y) + coef[i]) } print(BLUPa) #---BLUP Estimates Narrow BLUPaNar <- NULL for( i in 1:length(coef)) { BLUPaNar[i] <- (mean(DataSet9.2$y) + coef[i]) } BLUPaNardata(DataSet9.2) DataSet9.2$a <- factor(x = DataSet9.2$a) DataSet9.2$b <- factor(x = DataSet9.2$b) library(lmerTest) Exam9.2lmer <- lmer(y ~ (1|b/a), data = DataSet9.2) summary(Exam9.2lmer) Exam9.2lmer2 <- lm(y ~ a + b %in% a, data = DataSet9.2) summary(Exam9.2lmer2) ##--- Over all mean library(phia) list9.2 <- list(a = c("1" = 1/7,"2" = 1/7 , "3" = 1/7,"4" = 1/7 , "5" = 1/7,"6" = 1/7 , "7" = 1/7 )) phia::testFactors(model = Exam9.2lmer2, levels = list9.2) #---BLUP Estimates coef <- unlist(ranef(Exam9.2lmer)$a) BLUPa <- NULL for(i in 1:length(coef)){ BLUPa[i] <- (mean(DataSet9.2$y) + coef[i]) } print(BLUPa) #---BLUP Estimates Narrow BLUPaNar <- NULL for( i in 1:length(coef)) { BLUPaNar[i] <- (mean(DataSet9.2$y) + coef[i]) } BLUPaNar
Exam9.4 Relationship between BLUP and Fixed Effect Estimators
Muhammad Yaseen ([email protected])
Adeela Munawar ([email protected])
Stroup, W. W. (2012). Generalized Linear Mixed Models: Modern Concepts, Methods and Applications. CRC Press.
data(DataSet9.4) DataSet9.4$a <- factor(x = DataSet9.4$a) DataSet9.4$b <- factor(x = DataSet9.4$b) library(lmerTest) Exam9.4lmer <- lmer(y ~ a + (1|b) + (1|a/b), data = DataSet9.4) summary(Exam9.4lmer) library(emmeans) emmeans(Exam9.4lmer, spec = ~a)data(DataSet9.4) DataSet9.4$a <- factor(x = DataSet9.4$a) DataSet9.4$b <- factor(x = DataSet9.4$b) library(lmerTest) Exam9.4lmer <- lmer(y ~ a + (1|b) + (1|a/b), data = DataSet9.4) summary(Exam9.4lmer) library(emmeans) emmeans(Exam9.4lmer, spec = ~a)
Table1.1 is used for inspecting probability distribution and to define a plausible process.
data(Table1.1)data(Table1.1)
A data.frame with 11 rows and 3 variables.
x independent variable
Nx bernouli trials (bernouli outcomes on each individual)
y number of successes on each individual
Muhammad Yaseen ([email protected])
Adeela Munawar ([email protected])
Stroup, W. W. (2012).Generalized linear mixed models: modern concepts, methods and applications. CRC press.
library(StroupGLMM) data(Table1.1)library(StroupGLMM) data(Table1.1)
Exam1.2 is used to see types of model effects by plotting regression data
data(Table1.2)data(Table1.2)
A data.frame with 36 rows and 5 variables.
X have 11 levels in varying intervals from 0 to 48 observed for multiple batches
Y continuous variable observed at each level of X
Fav number of successes
N number of bernoulli trials
Batch Batches as 1, 2, 3, 4
Muhammad Yaseen ([email protected])
Adeela Munawar ([email protected])
Stroup, W. W. (2012).Generalized linear mixed models: modern concepts, methods and applications. CRC press.
data(Table1.2)data(Table1.2)