Generalized Linear Mixed-Effects for 2-level Fractional Factorial Design

1 Generalized Linear Mixed-Effects Models: lme4
Reference

1 Generalized Linear Mixed-Effects Models: lme4

Reed Johnson et al. (2013)
Franzén, Dinnétz, and Hammer (2016)
Grilli and Rampichini (2015)
Bates et al. (2015)
ML is better for unbalanced data, but it produces biased results.
REML is unbiased, but it cannot be used when comparing two nested models with a likelihood ratio test.

1.1 Specifying the random effects in multilevel models: beyond standard assumptions

Examplesof mixed-effects model formulas: sourced from @Bates2015

Examplesof mixed-effects model formulas: sourced from Bates et al. (2015)

Profiles for Full Factorial and Fractional Factorial Designs: sourced from @Ueniversitesi2015

Profiles for Full Factorial and Fractional Factorial Designs: sourced from Üniversitesi et al. (2015)

library(readxl)
ChoiceExperiment = read_excel("ChoiceExperiment.xls")
ChoiceExperimentDF = as.data.frame(ChoiceExperiment)

1.2 transform five attributes as factor levels

ChoiceExperimentDF$CHOICE = as.factor(ChoiceExperimentDF$CHOICE)

ChoiceExperimentDF$SUBSIDIES = as.factor(ChoiceExperimentDF$SUBSIDIES)
ChoiceExperimentDF$TIMEFRAME = as.factor(ChoiceExperimentDF$TIMEFRAME)
ChoiceExperimentDF$SUPPORT = as.factor(ChoiceExperimentDF$SUPPORT)
ChoiceExperimentDF$COMPENSATION = as.factor(ChoiceExperimentDF$COMPENSATION)
ChoiceExperimentDF$CEILING = as.factor(ChoiceExperimentDF$CEILING)

ChoiceExperimentDF$ID = as.factor(ChoiceExperimentDF$ID)

str(ChoiceExperimentDF)

'data.frame':   116 obs. of  7 variables:
 $ ID          : Factor w/ 29 levels "4","14","23",..: 2 2 2 2 3 3 3 3 5 5 ...
 $ CHOICE      : Factor w/ 2 levels "0","1": 1 2 2 1 2 2 2 2 1 2 ...
 $ SUBSIDIES   : Factor w/ 2 levels "0","1": 2 2 1 1 2 2 1 1 2 2 ...
 $ TIMEFRAME   : Factor w/ 2 levels "0","1": 2 1 2 1 2 1 2 1 2 1 ...
 $ SUPPORT     : Factor w/ 2 levels "0","1": 1 2 2 2 1 2 2 2 1 2 ...
 $ COMPENSATION: Factor w/ 2 levels "0","1": 1 2 1 1 1 2 1 1 1 2 ...
 $ CEILING     : Factor w/ 2 levels "0","1": 1 2 2 1 1 2 2 1 1 2 ...

1.3 Fractional Fractorial Design for Mixed Effect Model

total number of sampling size = 116
number of replications per each ID = 4
number of unique ID case = 116/4 = 29

1.4 Seting up a model: Generalized linear mixed model

Model.01 <- glmer(CHOICE ~ SUBSIDIES + TIMEFRAME + SUPPORT + COMPENSATION + CEILING + 
    (1 | ID), family = binomial(link = "logit"), REML = T, ChoiceExperimentDF)
summary(Model.01)

Generalized linear mixed model fit by maximum likelihood (Laplace
  Approximation) [glmerMod]
 Family: binomial  ( logit )
Formula: CHOICE ~ SUBSIDIES + TIMEFRAME + SUPPORT + COMPENSATION + CEILING +  
    (1 | ID)
   Data: ChoiceExperimentDF

     AIC      BIC   logLik deviance df.resid 
   137.1    156.3    -61.5    123.1      109 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.9987 -0.6302  0.2026  0.5866  2.0604 

Random effects:
 Groups Name        Variance Std.Dev.
 ID     (Intercept) 0.5822   0.763   
Number of obs: 116, groups:  ID, 29

Fixed effects:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept)   -2.15464    0.66119  -3.259 0.001119 ** 
SUBSIDIES1     1.19939    0.49915   2.403 0.016268 *  
TIMEFRAME1     0.04066    0.47564   0.085 0.931882    
SUPPORT1       0.15320    0.49831   0.307 0.758516    
COMPENSATION1  2.10216    0.54718   3.842 0.000122 ***
CEILING1       1.33760    0.49703   2.691 0.007120 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) SUBSID TIMEFR SUPPOR COMPEN
SUBSIDIES1  -0.552                            
TIMEFRAME1  -0.384 -0.041                     
SUPPORT1    -0.385  0.224 -0.066              
COMPENSATIO -0.510  0.183  0.088  0.021       
CEILING1    -0.453  0.181 -0.037 -0.102  0.208

Anova(Model.01)

Analysis of Deviance Table (Type II Wald chisquare tests)

Response: CHOICE
               Chisq Df Pr(>Chisq)    
SUBSIDIES     5.7737  1  0.0162680 *  
TIMEFRAME     0.0073  1  0.9318825    
SUPPORT       0.0945  1  0.7585163    
COMPENSATION 14.7593  1  0.0001221 ***
CEILING       7.2425  1  0.0071196 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Model.02 <- glmer(CHOICE ~ SUBSIDIES + COMPENSATION + CEILING + (1 | ID), family = binomial(link = "logit"), 
    REML = T, ChoiceExperimentDF)
summary(Model.02)

Generalized linear mixed model fit by maximum likelihood (Laplace
  Approximation) [glmerMod]
 Family: binomial  ( logit )
Formula: CHOICE ~ SUBSIDIES + COMPENSATION + CEILING + (1 | ID)
   Data: ChoiceExperimentDF

     AIC      BIC   logLik deviance df.resid 
   133.2    146.9    -61.6    123.2      111 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.8797 -0.6125  0.2096  0.6088  2.1136 

Random effects:
 Groups Name        Variance Std.Dev.
 ID     (Intercept) 0.5677   0.7534  
Number of obs: 116, groups:  ID, 29

Fixed effects:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept)    -2.0504     0.5488  -3.736 0.000187 ***
SUBSIDIES1      1.1688     0.4860   2.405 0.016178 *  
COMPENSATION1   2.0967     0.5469   3.834 0.000126 ***
CEILING1        1.3569     0.4936   2.749 0.005981 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) SUBSID COMPEN
SUBSIDIES1  -0.598              
COMPENSATIO -0.571  0.184       
CEILING1    -0.620  0.214  0.223

Anova(Model.02)

Analysis of Deviance Table (Type II Wald chisquare tests)

Response: CHOICE
               Chisq Df Pr(>Chisq)    
SUBSIDIES     5.7834  1  0.0161784 *  
COMPENSATION 14.6970  1  0.0001262 ***
CEILING       7.5559  1  0.0059814 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

1.5 Diagnostic Plots

plot(Model.02, type = c("p", "smooth"), main = "Residuals vs Fitted values")

qqmath(Model.02, id = 0.05, main = "Q-Q plot")

1.6 Effects Plots

eff1 <- allEffects(Model.02)
plot(eff1)

plot(eff1, type = "link")

e1 <- Effect(c("SUBSIDIES"), Model.01)
e2 <- Effect(c("TIMEFRAME"), Model.01)
e3 <- Effect(c("SUPPORT"), Model.01)
e4 <- Effect(c("COMPENSATION"), Model.01)
e5 <- Effect(c("CEILING"), Model.02)

par(mar = c(5, 4, 2, 2))
a <- summary(e1)
plotCI(1:2, a$effect, ui = a$upper, li = a$lower, xaxt = "n", xlim = c(1, 10), ylim = c(0, 
    1), xlab = "", ylab = "Willing to create a wetland", las = 1, cex.lab = 0.7, 
    cex.axis = 0.5, bty = "n", gap = 0.02, col = "red")
axis(1, 1:2, c("Current", "Higher"), cex.axis = 0.5)
text(1.5, 0.5, "*")
mtext("Subsidy", 1, line = 2, at = 1.5, cex = 0.7)

b <- summary(e2)
par(new = TRUE)
plotCI(3:4, b$effect, ui = b$upper, li = b$lower, xaxt = "n", yaxt = "n", xlim = c(1, 
    10), ylim = c(0, 1), xlab = "", ylab = "", las = 1, cex.axis = 0.7, bty = "n", 
    gap = 0.02, col = "blue")
axis(1, 3:4, c("Current", "Longer"), cex.axis = 0.5)
text(3.5, 0.5, "ns", cex = 0.7)
mtext("Time frame", 1, line = 2, at = 3.5, cex = 0.7)

c <- summary(e3)
par(new = TRUE)
plotCI(5:6, c$effect, ui = c$upper, li = c$lower, xaxt = "n", yaxt = "n", xlim = c(1, 
    10), ylim = c(0, 1), xlab = "", ylab = "", las = 1, cex.axis = 0.7, bty = "n", 
    gap = 0.02, col = "blue")
axis(1, 5:6, c("No", "Yes"), cex.axis = 0.5)
text(5.5, 0.5, "ns", cex = 0.7)
mtext("Support", 1, line = 2, at = 5.5, cex = 0.7)

d <- summary(e4)
par(new = TRUE)
plotCI(7:8, d$effect, ui = d$upper, li = d$lower, xaxt = "n", yaxt = "n", xlim = c(1, 
    10), ylim = c(0, 1), xlab = "", ylab = "", las = 1, cex.axis = 0.7, bty = "n", 
    gap = 0.02, col = "red")
axis(1, 7:8, c("Current", "Higher"), cex.axis = 0.5)
text(7.5, 0.5, "***")
mtext("Compensation", 1, line = 2, at = 7.5, cex = 0.7)

e <- summary(e5)
par(new = TRUE)
plotCI(9:10, e$effect, ui = e$upper, li = e$lower, xaxt = "n", yaxt = "n", xlim = c(1, 
    10), ylim = c(0, 1), xlab = "", ylab = "", las = 1, cex.axis = 0.7, bty = "n", 
    gap = 0.02, col = "red")
axis(1, 9:10, c("Current", "Higher"), cex.axis = 0.5)
text(9.5, 0.5, "**")
mtext("Ceiling", 1, line = 2, at = 9.5, cex = 0.7)

grid()

anova(Model.02, Model.01)

Data: ChoiceExperimentDF
Models:
Model.02: CHOICE ~ SUBSIDIES + COMPENSATION + CEILING + (1 | ID)
Model.01: CHOICE ~ SUBSIDIES + TIMEFRAME + SUPPORT + COMPENSATION + CEILING + 
Model.01:     (1 | ID)
         Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
Model.02  5 133.17 146.94 -61.585   123.17                         
Model.01  7 137.06 156.34 -61.532   123.06 0.1068      2      0.948

1.7 Results:

Model_02 is better than Model_01

Reference

Bates, Douglas, Martin Mächler, Ben Bolker, and Steve Walker. 2015. “Fitting Linear Mixed-Effects Models Using Lme4.” Journal of Statistical Software; Vol 1, Issue 1 (2015). https://www.jstatsoft.org/v067/i01.

Franzén, Frida, Patrik Dinnétz, and Monica Hammer. 2016. “Factors Affecting Farmers’ Willingness to Participate in Eutrophication Mitigation – a Case Study of Preferences for Wetland Creation in Sweden.” Ecological Economics 130 (October): 8–15. http://www.sciencedirect.com/science/article/pii/S0921800916305961.

Grilli, Leonardo, and Carla Rampichini. 2015. “Specification of Random Effects in Multilevel Models: A Review.” Quality & Quantity 49 (3): 967–76. https://doi.org/10.1007/s11135-014-0060-5.

Reed Johnson, F., Emily Lancsar, Deborah Marshall, Vikram Kilambi, Axel Mühlbacher, Dean A. Regier, Brian W. Bresnahan, Barbara Kanninen, and John F. P. Bridges. 2013. “Constructing Experimental Designs for Discrete-Choice Experiments: Report of the Ispor Conjoint Analysis Experimental Design Good Research Practices Task Force.” Value in Health 16 (1): 3–13. http://www.sciencedirect.com/science/article/pii/S1098301512041629.

Üniversitesi, Trakya, İktisadi İdari, Bilimler Fakültesi, Nihat Taş, Nil Kodaz Engizek, Emrah Önder, and Güler Önder. 2015. “HEALTH Education Planning in Marketing Perspective Using Conjoint Analysis*.” Trakya Üniversitesi İktisadi Ve İdari Bilimler Fakültesi E-Dergi 4 (January): 40–66.

DK C.

2020-03-01