Skip to main content
Cross Validated

Return to Answer

added 14 characters in body
Source Link
DrJerryTAO
DrJerryTAO
  • 2.3k
  • 6
  • 22
library(dplyr)
Data <- data.frame(
  Rating = 1:30, 
  Condition = factor(c(rep("L", 10), rep("M", 10), rep("LM", 10)), levels = c(
    "L", "M", "LM")), 
  L01 = c(rep(1, 10), rep(0, 10), rep(1, 10)),
  M01 = c(rep(0, 10), rep(1, 10), rep(1, 10)), 
  LTrue = c(rep(TRUE, 10), rep(FALSE, 10), rep(TRUE, 10)),
  MTrue = c(rep(FALSE, 10), rep(TRUE, 10), rep(TRUE, 10)), 
  Lny = factor(c(rep("Yes", 10), rep("No", 10), rep("Yes", 10)), levels = c(
    "No", "Yes")),
  Mny = factor(c(rep("No", 10), rep("Yes", 10), rep("Yes", 10)), levels = c(
    "No", "Yes")))
Data |>
  group_by(Condition) |>
  summarise(mean = mean(Rating))
"# A tibble: 3 × 2
  Condition  mean
  <fct>     <dbl>
1 L           5.5
2 M          15.5
3 LM         25.5"

# Factor treatment contrast
summary(lm(Rating ~ Condition, data = Data))
"            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   5.5000     0.9574   5.745 4.16e-06 ***
ConditionM   10.0000     1.3540   7.385 6.05e-08 ***
ConditionLM  20.0000     1.3540  14.771 1.87e-14 ***
Residual standard error: 3.028 on 27 degrees of freedom
Multiple R-squared:  0.8899,    Adjusted R-squared:  0.8817 
F-statistic: 109.1 on 2 and 27 DF,  p-value: 1.162e-13"

# 0/1 encoding
summary(lm(Rating ~ 1 + L01 * M01, data = Data)) # weird coef
"            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   -4.500      1.658  -2.714   0.0114 *  
L01           10.000      1.354   7.385 6.05e-08 ***
M01           20.000      1.354  14.771 1.87e-14 ***
L01:M01           NA         NA      NA       NA
Residual standard error: 3.028 on 27 degrees of freedom
Multiple R-squared:  0.8899,    Adjusted R-squared:  0.8817 
F-statistic: 109.1 on 2 and 27 DF,  p-value: 1.162e-13"
summary(lm(Rating ~ 0 + L01 * M01, data = Data)) # wrong R2 
"        Estimate Std. Error t value Pr(>|t|)    
L01       5.5000     0.9574   5.745 4.16e-06 ***
M01      15.5000     0.9574  16.189 2.00e-15 ***
L01:M01   4.5000     1.6583   2.714   0.0114 *
Residual standard error: 3.028 on 27 degrees of freedom
Multiple R-squared:  0.9738,    Adjusted R-squared:  0.9709 
F-statistic: 334.8 on 3 and 27 DF,  p-value: < 2.2e-16"
summary(lm(Rating ~ M01 + L01 : M01, data = Data)) # correct formula
"            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   5.5000     0.9574   5.745 4.16e-06 ***
M01          10.0000     1.3540   7.385 6.05e-08 ***
M01:L01      10.0000     1.3540   7.385 6.05e-08 ***
Residual standard error: 3.028 on 27 degrees of freedom
Multiple R-squared:  0.8899,    Adjusted R-squared:  0.8817 
F-statistic: 109.1 on 2 and 27 DF,  p-value: 1.162e-13"

# Logical encoding
summary(lm(Rating ~ 1 + LTrue * MTrue, data = Data)) # weird coef
"                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)           -4.500      1.658  -2.714   0.0114 *  
LTrueTRUE             10.000      1.354   7.385 6.05e-08 ***
MTrueTRUE             20.000      1.354  14.771 1.87e-14 ***
LTrueTRUE:MTrueTRUE       NA         NA      NA       NA    
Residual standard error: 3.028 on 27 degrees of freedom
Multiple R-squared:  0.8899,    Adjusted R-squared:  0.8817 
F-statistic: 109.1 on 2 and 27 DF,  p-value: 1.162e-13"
summary(lm(Rating ~ 0 + LTrue * MTrue, data = Data)) # weird coef + wrong R2
"                    Estimate Std. Error t value Pr(>|t|)    
LTrueFALSE           -4.5000     1.6583  -2.714   0.0114 *  
LTrueTRUE             5.5000     0.9574   5.745 4.16e-06 ***
MTrueTRUE            20.0000     1.3540  14.771 1.87e-14 ***
LTrueTRUE:MTrueTRUE       NA         NA      NA       NA    
Residual standard error: 3.028 on 27 degrees of freedom
Multiple R-squared:  0.9738,    Adjusted R-squared:  0.9709 
F-statistic: 334.8 on 3 and 27 DF,  p-value: < 2.2e-16"
summary(lm(Rating ~ MTrue + LTrue : MTrue, data = Data)) # correct formula
"                     Estimate Std. Error t value Pr(>|t|)    
(Intercept)            5.5000     0.9574   5.745 4.16e-06 ***
MTrueTRUE             10.0000     1.3540   7.385 6.05e-08 ***
MTrueFALSE:LTrueTRUE       NA         NA      NA       NA    
MTrueTRUE:LTrueTRUE   10.0000     1.3540   7.385 6.05e-08 ***
Residual standard error: 3.028 on 27 degrees of freedom
Multiple R-squared:  0.8899,    Adjusted R-squared:  0.8817 
F-statistic: 109.1 on 2 and 27 DF,  p-value: 1.162e-13"
summary(lm(Rating ~ MTrue + I(LTrue * MTrue), data = Data)) # Removes NA above
"                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)        5.5000     0.9574   5.745 4.16e-06 ***
MTrueTRUE         10.0000     1.3540   7.385 6.05e-08 ***
I(LTrue * MTrue)  10.0000     1.3540   7.385 6.05e-08 ***
Residual standard error: 3.028 on 27 degrees of freedom
Multiple R-squared:  0.8899,    Adjusted R-squared:  0.8817 
F-statistic: 109.1 on 2 and 27 DF,  p-value: 1.162e-13"

# Binary factor encoding
summary(lm(Rating ~ 1 + Lny * Mny, data = Data)) # weird coef
"              Estimate Std. Error t value Pr(>|t|)    
(Intercept)     -4.500      1.658  -2.714   0.0114 *  
LnyYes          10.000      1.354   7.385 6.05e-08 ***
MnyYes          20.000      1.354  14.771 1.87e-14 ***
LnyYes:MnyYes       NA         NA      NA       NA    
Residual standard error: 3.028 on 27 degrees of freedom
Multiple R-squared:  0.8899,    Adjusted R-squared:  0.8817 
F-statistic: 109.1 on 2 and 27 DF,  p-value: 1.162e-13"
summary(lm(Rating ~ 0 + Lny * Mny, data = Data)) # weird coef + wrong R2
"              Estimate Std. Error t value Pr(>|t|)    
LnyNo          -4.5000     1.6583  -2.714   0.0114 *  
LnyYes          5.5000     0.9574   5.745 4.16e-06 ***
MnyYes         20.0000     1.3540  14.771 1.87e-14 ***
LnyYes:MnyYes       NA         NA      NA       NA    
Residual standard error: 3.028 on 27 degrees of freedom
Multiple R-squared:  0.9738,    Adjusted R-squared:  0.9709 
F-statistic: 334.8 on 3 and 27 DF,  p-value: < 2.2e-16"
summary(lm(Rating ~ Mny + Lny : Mny, data = Data)) # correct formula
"              Estimate Std. Error t value Pr(>|t|)    
(Intercept)     5.5000     0.9574   5.745 4.16e-06 ***
MnyYes         10.0000     1.3540   7.385 6.05e-08 ***
MnyNo:LnyYes        NA         NA      NA       NA    
MnyYes:LnyYes  10.0000     1.3540   7.385 6.05e-08 ***
Residual standard error: 3.028 on 27 degrees of freedom
Multiple R-squared:  0.8899,    Adjusted R-squared:  0.8817 
F-statistic: 109.1 on 2 and 27 DF,  p-value: 1.162e-13"
summary(lm(Rating ~ Mny + I(Lny : Mny), data = Data)) # NA still there
"                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)         5.5000     0.9574   5.745 4.16e-06 ***
MnyYes             10.0000     1.3540   7.385 6.05e-08 ***
I(Lny:Mny)Yes:No        NA         NA      NA       NA    
I(Lny:Mny)Yes:Yes  10.0000     1.3540   7.385 6.05e-08 ***
Residual standard error: 3.028 on 27 degrees of freedom
Multiple R-squared:  0.8899,    Adjusted R-squared:  0.8817 
F-statistic: 109.1 on 2 and 27 DF,  p-value: 1.162e-13"
summary(lm(Rating ~ Mny + I(Lny == "Yes" & Mny == "Yes"), data = Data)) # no NA
"                                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)                          5.5000     0.9574   5.745 4.16e-06 ***
MnyYes                              10.0000     1.3540   7.385 6.05e-08 ***
I(Lny == Yes & Mny == Yes)TRUE      10.0000     1.3540   7.385 6.05e-08 ***
Residual standard error: 3.028 on 27 degrees of freedom
Multiple R-squared:  0.8899,    Adjusted R-squared:  0.8817 
F-statistic: 109.1 on 2 and 27 DF,  p-value: 1.162e-13"

AIC is used for model comparison and selection. It does not test null hypothesis. To test specific hypotheses of whether one or more coefficients are significant, likelihood ratio test is a good choice. Its implementation in models with random effects have some caveats. Linear mixed models are usually fitted with the restricted maximum likelihood estimator, but this restricted maximum likelihood as the model results cannot be used in likelihood-ratio tests of fixed effects. We need to switch to maximum likelihood, usually by specifying the argument REML = FALSE or method = ML in the estimating function before testing with anova(). To test random effects, in contrast, we are usually interested in seeing if the standard deviation of the random term is zero or larger than zero. This test is on the boundary of the possible range of the parameter and the regular likelihood-ratio test such as anova() is inappropriate. This is because the test statistic under the null hypothesis $H_0: \sigma = 0$ is not $\chi^2(1)$ distributed but a mixture of two or more $\chi^2$ of different degrees of freedom, $\chi^2(0)/2 + \chi^2(1)/2$ in this one-variance case. It requires specific likelihood-ratio test methods for variance components, usually included in the package. See Molenberghs, G., & Verbeke, G. (2007). Likelihood ratio, score, and Wald tests in a constrained parameter space. The American StatisticianThe American Statistician, 61(1), 22–27. jstor.org/stable/27643833 https://www.jstor.org/stable/27643833.

library(dplyr)
Data <- data.frame(
  Rating = 1:30, 
  Condition = factor(c(rep("L", 10), rep("M", 10), rep("LM", 10)), levels = c(
    "L", "M", "LM")), 
  L01 = c(rep(1, 10), rep(0, 10), rep(1, 10)),
  M01 = c(rep(0, 10), rep(1, 10), rep(1, 10)), 
  LTrue = c(rep(TRUE, 10), rep(FALSE, 10), rep(TRUE, 10)),
  MTrue = c(rep(FALSE, 10), rep(TRUE, 10), rep(TRUE, 10)), 
  Lny = factor(c(rep("Yes", 10), rep("No", 10), rep("Yes", 10)), levels = c(
    "No", "Yes")),
  Mny = factor(c(rep("No", 10), rep("Yes", 10), rep("Yes", 10)), levels = c(
    "No", "Yes")))
Data |>
  group_by(Condition) |>
  summarise(mean = mean(Rating))
"# A tibble: 3 × 2
  Condition  mean
  <fct>     <dbl>
1 L           5.5
2 M          15.5
3 LM         25.5"

# Factor treatment contrast
summary(lm(Rating ~ Condition, data = Data))
"            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   5.5000     0.9574   5.745 4.16e-06 ***
ConditionM   10.0000     1.3540   7.385 6.05e-08 ***
ConditionLM  20.0000     1.3540  14.771 1.87e-14 ***
Residual standard error: 3.028 on 27 degrees of freedom
Multiple R-squared:  0.8899,    Adjusted R-squared:  0.8817 
F-statistic: 109.1 on 2 and 27 DF,  p-value: 1.162e-13"

# 0/1 encoding
summary(lm(Rating ~ 1 + L01 * M01, data = Data)) # weird coef
"            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   -4.500      1.658  -2.714   0.0114 *  
L01           10.000      1.354   7.385 6.05e-08 ***
M01           20.000      1.354  14.771 1.87e-14 ***
L01:M01           NA         NA      NA       NA
Residual standard error: 3.028 on 27 degrees of freedom
Multiple R-squared:  0.8899,    Adjusted R-squared:  0.8817 
F-statistic: 109.1 on 2 and 27 DF,  p-value: 1.162e-13"
summary(lm(Rating ~ 0 + L01 * M01, data = Data)) # wrong R2 
"        Estimate Std. Error t value Pr(>|t|)    
L01       5.5000     0.9574   5.745 4.16e-06 ***
M01      15.5000     0.9574  16.189 2.00e-15 ***
L01:M01   4.5000     1.6583   2.714   0.0114 *
Residual standard error: 3.028 on 27 degrees of freedom
Multiple R-squared:  0.9738,    Adjusted R-squared:  0.9709 
F-statistic: 334.8 on 3 and 27 DF,  p-value: < 2.2e-16"
summary(lm(Rating ~ M01 + L01 : M01, data = Data)) # correct formula
"            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   5.5000     0.9574   5.745 4.16e-06 ***
M01          10.0000     1.3540   7.385 6.05e-08 ***
M01:L01      10.0000     1.3540   7.385 6.05e-08 ***
Residual standard error: 3.028 on 27 degrees of freedom
Multiple R-squared:  0.8899,    Adjusted R-squared:  0.8817 
F-statistic: 109.1 on 2 and 27 DF,  p-value: 1.162e-13"

summary(lm(Rating ~ 1 + LTrue * MTrue, data = Data)) # weird coef
"                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)           -4.500      1.658  -2.714   0.0114 *  
LTrueTRUE             10.000      1.354   7.385 6.05e-08 ***
MTrueTRUE             20.000      1.354  14.771 1.87e-14 ***
LTrueTRUE:MTrueTRUE       NA         NA      NA       NA    
Residual standard error: 3.028 on 27 degrees of freedom
Multiple R-squared:  0.8899,    Adjusted R-squared:  0.8817 
F-statistic: 109.1 on 2 and 27 DF,  p-value: 1.162e-13"
summary(lm(Rating ~ 0 + LTrue * MTrue, data = Data)) # weird coef + wrong R2
"                    Estimate Std. Error t value Pr(>|t|)    
LTrueFALSE           -4.5000     1.6583  -2.714   0.0114 *  
LTrueTRUE             5.5000     0.9574   5.745 4.16e-06 ***
MTrueTRUE            20.0000     1.3540  14.771 1.87e-14 ***
LTrueTRUE:MTrueTRUE       NA         NA      NA       NA    
Residual standard error: 3.028 on 27 degrees of freedom
Multiple R-squared:  0.9738,    Adjusted R-squared:  0.9709 
F-statistic: 334.8 on 3 and 27 DF,  p-value: < 2.2e-16"
summary(lm(Rating ~ MTrue + LTrue : MTrue, data = Data)) # correct formula
"                     Estimate Std. Error t value Pr(>|t|)    
(Intercept)            5.5000     0.9574   5.745 4.16e-06 ***
MTrueTRUE             10.0000     1.3540   7.385 6.05e-08 ***
MTrueFALSE:LTrueTRUE       NA         NA      NA       NA    
MTrueTRUE:LTrueTRUE   10.0000     1.3540   7.385 6.05e-08 ***
Residual standard error: 3.028 on 27 degrees of freedom
Multiple R-squared:  0.8899,    Adjusted R-squared:  0.8817 
F-statistic: 109.1 on 2 and 27 DF,  p-value: 1.162e-13"
summary(lm(Rating ~ MTrue + I(LTrue * MTrue), data = Data)) # Removes NA above
"                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)        5.5000     0.9574   5.745 4.16e-06 ***
MTrueTRUE         10.0000     1.3540   7.385 6.05e-08 ***
I(LTrue * MTrue)  10.0000     1.3540   7.385 6.05e-08 ***
Residual standard error: 3.028 on 27 degrees of freedom
Multiple R-squared:  0.8899,    Adjusted R-squared:  0.8817 
F-statistic: 109.1 on 2 and 27 DF,  p-value: 1.162e-13"

summary(lm(Rating ~ 1 + Lny * Mny, data = Data)) # weird coef
"              Estimate Std. Error t value Pr(>|t|)    
(Intercept)     -4.500      1.658  -2.714   0.0114 *  
LnyYes          10.000      1.354   7.385 6.05e-08 ***
MnyYes          20.000      1.354  14.771 1.87e-14 ***
LnyYes:MnyYes       NA         NA      NA       NA    
Residual standard error: 3.028 on 27 degrees of freedom
Multiple R-squared:  0.8899,    Adjusted R-squared:  0.8817 
F-statistic: 109.1 on 2 and 27 DF,  p-value: 1.162e-13"
summary(lm(Rating ~ 0 + Lny * Mny, data = Data)) # weird coef + wrong R2
"              Estimate Std. Error t value Pr(>|t|)    
LnyNo          -4.5000     1.6583  -2.714   0.0114 *  
LnyYes          5.5000     0.9574   5.745 4.16e-06 ***
MnyYes         20.0000     1.3540  14.771 1.87e-14 ***
LnyYes:MnyYes       NA         NA      NA       NA    
Residual standard error: 3.028 on 27 degrees of freedom
Multiple R-squared:  0.9738,    Adjusted R-squared:  0.9709 
F-statistic: 334.8 on 3 and 27 DF,  p-value: < 2.2e-16"
summary(lm(Rating ~ Mny + Lny : Mny, data = Data)) # correct formula
"              Estimate Std. Error t value Pr(>|t|)    
(Intercept)     5.5000     0.9574   5.745 4.16e-06 ***
MnyYes         10.0000     1.3540   7.385 6.05e-08 ***
MnyNo:LnyYes        NA         NA      NA       NA    
MnyYes:LnyYes  10.0000     1.3540   7.385 6.05e-08 ***
Residual standard error: 3.028 on 27 degrees of freedom
Multiple R-squared:  0.8899,    Adjusted R-squared:  0.8817 
F-statistic: 109.1 on 2 and 27 DF,  p-value: 1.162e-13"
summary(lm(Rating ~ Mny + I(Lny : Mny), data = Data)) # NA still there
"                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)         5.5000     0.9574   5.745 4.16e-06 ***
MnyYes             10.0000     1.3540   7.385 6.05e-08 ***
I(Lny:Mny)Yes:No        NA         NA      NA       NA    
I(Lny:Mny)Yes:Yes  10.0000     1.3540   7.385 6.05e-08 ***
Residual standard error: 3.028 on 27 degrees of freedom
Multiple R-squared:  0.8899,    Adjusted R-squared:  0.8817 
F-statistic: 109.1 on 2 and 27 DF,  p-value: 1.162e-13"
summary(lm(Rating ~ Mny + I(Lny == "Yes" & Mny == "Yes"), data = Data)) # no NA
"                                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)                          5.5000     0.9574   5.745 4.16e-06 ***
MnyYes                              10.0000     1.3540   7.385 6.05e-08 ***
I(Lny == Yes & Mny == Yes)TRUE      10.0000     1.3540   7.385 6.05e-08 ***
Residual standard error: 3.028 on 27 degrees of freedom
Multiple R-squared:  0.8899,    Adjusted R-squared:  0.8817 
F-statistic: 109.1 on 2 and 27 DF,  p-value: 1.162e-13"

AIC is used for model comparison and selection. It does not test null hypothesis. To test specific hypotheses of whether one or more coefficients are significant, likelihood ratio test is a good choice. Its implementation in models with random effects have some caveats. Linear mixed models are usually fitted with the restricted maximum likelihood estimator, but this restricted maximum likelihood as the model results cannot be used in likelihood-ratio tests of fixed effects. We need to switch to maximum likelihood, usually by specifying the argument REML = FALSE or method = ML in the estimating function before testing with anova(). To test random effects, in contrast, we are usually interested in seeing if the standard deviation of the random term is zero or larger than zero. This test is on the boundary of the possible range of the parameter and the regular likelihood-ratio test such as anova() is inappropriate. This is because the test statistic under the null hypothesis $H_0: \sigma = 0$ is not $\chi^2(1)$ distributed but a mixture of two or more $\chi^2$ of different degrees of freedom, $\chi^2(0)/2 + \chi^2(1)/2$ in this one-variance case. It requires specific likelihood-ratio test methods for variance components, usually included in the package. See Molenberghs, G., & Verbeke, G. (2007). Likelihood ratio, score, and Wald tests in a constrained parameter space. The American Statistician, 61(1), 22–27. jstor.org/stable/27643833.

library(dplyr)
Data <- data.frame(
  Rating = 1:30, 
  Condition = factor(c(rep("L", 10), rep("M", 10), rep("LM", 10)), levels = c(
    "L", "M", "LM")), 
  L01 = c(rep(1, 10), rep(0, 10), rep(1, 10)),
  M01 = c(rep(0, 10), rep(1, 10), rep(1, 10)), 
  LTrue = c(rep(TRUE, 10), rep(FALSE, 10), rep(TRUE, 10)),
  MTrue = c(rep(FALSE, 10), rep(TRUE, 10), rep(TRUE, 10)), 
  Lny = factor(c(rep("Yes", 10), rep("No", 10), rep("Yes", 10)), levels = c(
    "No", "Yes")),
  Mny = factor(c(rep("No", 10), rep("Yes", 10), rep("Yes", 10)), levels = c(
    "No", "Yes")))
Data |>
  group_by(Condition) |>
  summarise(mean = mean(Rating))
"# A tibble: 3 × 2
  Condition  mean
  <fct>     <dbl>
1 L           5.5
2 M          15.5
3 LM         25.5"

# Factor treatment contrast
summary(lm(Rating ~ Condition, data = Data))
"            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   5.5000     0.9574   5.745 4.16e-06 ***
ConditionM   10.0000     1.3540   7.385 6.05e-08 ***
ConditionLM  20.0000     1.3540  14.771 1.87e-14 ***
Residual standard error: 3.028 on 27 degrees of freedom
Multiple R-squared:  0.8899,    Adjusted R-squared:  0.8817 
F-statistic: 109.1 on 2 and 27 DF,  p-value: 1.162e-13"

# 0/1 encoding
summary(lm(Rating ~ 1 + L01 * M01, data = Data)) # weird coef
"            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   -4.500      1.658  -2.714   0.0114 *  
L01           10.000      1.354   7.385 6.05e-08 ***
M01           20.000      1.354  14.771 1.87e-14 ***
L01:M01           NA         NA      NA       NA
Residual standard error: 3.028 on 27 degrees of freedom
Multiple R-squared:  0.8899,    Adjusted R-squared:  0.8817 
F-statistic: 109.1 on 2 and 27 DF,  p-value: 1.162e-13"
summary(lm(Rating ~ 0 + L01 * M01, data = Data)) # wrong R2 
"        Estimate Std. Error t value Pr(>|t|)    
L01       5.5000     0.9574   5.745 4.16e-06 ***
M01      15.5000     0.9574  16.189 2.00e-15 ***
L01:M01   4.5000     1.6583   2.714   0.0114 *
Residual standard error: 3.028 on 27 degrees of freedom
Multiple R-squared:  0.9738,    Adjusted R-squared:  0.9709 
F-statistic: 334.8 on 3 and 27 DF,  p-value: < 2.2e-16"
summary(lm(Rating ~ M01 + L01 : M01, data = Data)) # correct formula
"            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   5.5000     0.9574   5.745 4.16e-06 ***
M01          10.0000     1.3540   7.385 6.05e-08 ***
M01:L01      10.0000     1.3540   7.385 6.05e-08 ***
Residual standard error: 3.028 on 27 degrees of freedom
Multiple R-squared:  0.8899,    Adjusted R-squared:  0.8817 
F-statistic: 109.1 on 2 and 27 DF,  p-value: 1.162e-13"

# Logical encoding
summary(lm(Rating ~ 1 + LTrue * MTrue, data = Data)) # weird coef
"                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)           -4.500      1.658  -2.714   0.0114 *  
LTrueTRUE             10.000      1.354   7.385 6.05e-08 ***
MTrueTRUE             20.000      1.354  14.771 1.87e-14 ***
LTrueTRUE:MTrueTRUE       NA         NA      NA       NA    
Residual standard error: 3.028 on 27 degrees of freedom
Multiple R-squared:  0.8899,    Adjusted R-squared:  0.8817 
F-statistic: 109.1 on 2 and 27 DF,  p-value: 1.162e-13"
summary(lm(Rating ~ 0 + LTrue * MTrue, data = Data)) # weird coef + wrong R2
"                    Estimate Std. Error t value Pr(>|t|)    
LTrueFALSE           -4.5000     1.6583  -2.714   0.0114 *  
LTrueTRUE             5.5000     0.9574   5.745 4.16e-06 ***
MTrueTRUE            20.0000     1.3540  14.771 1.87e-14 ***
LTrueTRUE:MTrueTRUE       NA         NA      NA       NA    
Residual standard error: 3.028 on 27 degrees of freedom
Multiple R-squared:  0.9738,    Adjusted R-squared:  0.9709 
F-statistic: 334.8 on 3 and 27 DF,  p-value: < 2.2e-16"
summary(lm(Rating ~ MTrue + LTrue : MTrue, data = Data)) # correct formula
"                     Estimate Std. Error t value Pr(>|t|)    
(Intercept)            5.5000     0.9574   5.745 4.16e-06 ***
MTrueTRUE             10.0000     1.3540   7.385 6.05e-08 ***
MTrueFALSE:LTrueTRUE       NA         NA      NA       NA    
MTrueTRUE:LTrueTRUE   10.0000     1.3540   7.385 6.05e-08 ***
Residual standard error: 3.028 on 27 degrees of freedom
Multiple R-squared:  0.8899,    Adjusted R-squared:  0.8817 
F-statistic: 109.1 on 2 and 27 DF,  p-value: 1.162e-13"
summary(lm(Rating ~ MTrue + I(LTrue * MTrue), data = Data)) # Removes NA above
"                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)        5.5000     0.9574   5.745 4.16e-06 ***
MTrueTRUE         10.0000     1.3540   7.385 6.05e-08 ***
I(LTrue * MTrue)  10.0000     1.3540   7.385 6.05e-08 ***
Residual standard error: 3.028 on 27 degrees of freedom
Multiple R-squared:  0.8899,    Adjusted R-squared:  0.8817 
F-statistic: 109.1 on 2 and 27 DF,  p-value: 1.162e-13"

# Binary factor encoding
summary(lm(Rating ~ 1 + Lny * Mny, data = Data)) # weird coef
"              Estimate Std. Error t value Pr(>|t|)    
(Intercept)     -4.500      1.658  -2.714   0.0114 *  
LnyYes          10.000      1.354   7.385 6.05e-08 ***
MnyYes          20.000      1.354  14.771 1.87e-14 ***
LnyYes:MnyYes       NA         NA      NA       NA    
Residual standard error: 3.028 on 27 degrees of freedom
Multiple R-squared:  0.8899,    Adjusted R-squared:  0.8817 
F-statistic: 109.1 on 2 and 27 DF,  p-value: 1.162e-13"
summary(lm(Rating ~ 0 + Lny * Mny, data = Data)) # weird coef + wrong R2
"              Estimate Std. Error t value Pr(>|t|)    
LnyNo          -4.5000     1.6583  -2.714   0.0114 *  
LnyYes          5.5000     0.9574   5.745 4.16e-06 ***
MnyYes         20.0000     1.3540  14.771 1.87e-14 ***
LnyYes:MnyYes       NA         NA      NA       NA    
Residual standard error: 3.028 on 27 degrees of freedom
Multiple R-squared:  0.9738,    Adjusted R-squared:  0.9709 
F-statistic: 334.8 on 3 and 27 DF,  p-value: < 2.2e-16"
summary(lm(Rating ~ Mny + Lny : Mny, data = Data)) # correct formula
"              Estimate Std. Error t value Pr(>|t|)    
(Intercept)     5.5000     0.9574   5.745 4.16e-06 ***
MnyYes         10.0000     1.3540   7.385 6.05e-08 ***
MnyNo:LnyYes        NA         NA      NA       NA    
MnyYes:LnyYes  10.0000     1.3540   7.385 6.05e-08 ***
Residual standard error: 3.028 on 27 degrees of freedom
Multiple R-squared:  0.8899,    Adjusted R-squared:  0.8817 
F-statistic: 109.1 on 2 and 27 DF,  p-value: 1.162e-13"
summary(lm(Rating ~ Mny + I(Lny : Mny), data = Data)) # NA still there
"                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)         5.5000     0.9574   5.745 4.16e-06 ***
MnyYes             10.0000     1.3540   7.385 6.05e-08 ***
I(Lny:Mny)Yes:No        NA         NA      NA       NA    
I(Lny:Mny)Yes:Yes  10.0000     1.3540   7.385 6.05e-08 ***
Residual standard error: 3.028 on 27 degrees of freedom
Multiple R-squared:  0.8899,    Adjusted R-squared:  0.8817 
F-statistic: 109.1 on 2 and 27 DF,  p-value: 1.162e-13"
summary(lm(Rating ~ Mny + I(Lny == "Yes" & Mny == "Yes"), data = Data)) # no NA
"                                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)                          5.5000     0.9574   5.745 4.16e-06 ***
MnyYes                              10.0000     1.3540   7.385 6.05e-08 ***
I(Lny == Yes & Mny == Yes)TRUE      10.0000     1.3540   7.385 6.05e-08 ***
Residual standard error: 3.028 on 27 degrees of freedom
Multiple R-squared:  0.8899,    Adjusted R-squared:  0.8817 
F-statistic: 109.1 on 2 and 27 DF,  p-value: 1.162e-13"

AIC is used for model comparison and selection. It does not test null hypothesis. To test specific hypotheses of whether one or more coefficients are significant, likelihood ratio test is a good choice. Its implementation in models with random effects have some caveats. Linear mixed models are usually fitted with the restricted maximum likelihood estimator, but this restricted maximum likelihood as the model results cannot be used in likelihood-ratio tests of fixed effects. We need to switch to maximum likelihood, usually by specifying the argument REML = FALSE or method = ML in the estimating function before testing with anova(). To test random effects, in contrast, we are usually interested in seeing if the standard deviation of the random term is zero or larger than zero. This test is on the boundary of the possible range of the parameter and the regular likelihood-ratio test such as anova() is inappropriate. This is because the test statistic under the null hypothesis $H_0: \sigma = 0$ is not $\chi^2(1)$ distributed but a mixture of two or more $\chi^2$ of different degrees of freedom, $\chi^2(0)/2 + \chi^2(1)/2$ in this one-variance case. It requires specific likelihood-ratio test methods for variance components, usually included in the package. See Molenberghs, G., & Verbeke, G. (2007). Likelihood ratio, score, and Wald tests in a constrained parameter space. The American Statistician, 61(1), 22–27. https://www.jstor.org/stable/27643833.

added 150 characters in body
Source Link
DrJerryTAO
DrJerryTAO
  • 2.3k
  • 6
  • 22

AIC is used for model comparison and selection. It does not test null hypothesis. To test specific hypotheses of whether one or more coefficients are significant, likelihood ratio test is a good choice. Its implementation in models with random effects have some caveats. Linear mixed models are usually fitted with the restricted maximum likelihood estimator, but this restricted maximum likelihood as the model results cannot be used in likelihood-ratio tests of fixed effects. We need to switch to maximum likelihood, usually by specifying the argument REML = FALSE or method = ML in the estimating function. To test random effects, in contrast, we are usually interested in seeing if the standard deviation of the random term is zero or larger than zero. This test is on the boundary and should not be done by the regular likelihood-ratio test such as through before testing with anova(). To test random effects, in contrast, we are usually interested in seeing if the standard deviation of the random term is zero or larger than zero. This test is on the boundary of the possible range of the parameter and the regular likelihood-ratio test such as anova() is inappropriate. This is because the test statistic under the null hypothesis $H_0: \sigma = 0$ is not $\chi^2(1)$ distributed but a mixture of two or more $\chi^2$ of different degrees of freedom, $\chi^2(0)/2 + \chi^2(1)/2$ in this one-variance case. It requires specific likelihood-ratio test methods for variance components, usually included in the package. See Molenberghs, G., & Verbeke, G. (2007). Likelihood ratio, score, and Wald tests in a constrained parameter space. The American Statistician, 61(1), 22–27. jstor.org/stable/27643833.

AIC is used for model comparison and selection. It does not test null hypothesis. To test specific hypotheses of whether one or more coefficients are significant, likelihood ratio test is a good choice. Its implementation in models with random effects have some caveats. Linear mixed models are usually fitted with the restricted maximum likelihood estimator, but this restricted maximum likelihood as the model results cannot be used in likelihood-ratio tests of fixed effects. We need to switch to maximum likelihood, usually by specifying the argument REML = FALSE or method = ML in the estimating function. To test random effects, in contrast, we are usually interested in seeing if the standard deviation of the random term is zero or larger than zero. This test is on the boundary and should not be done by the regular likelihood-ratio test such as through anova(). It requires specific likelihood-ratio test methods for variance components, usually included in the package.

AIC is used for model comparison and selection. It does not test null hypothesis. To test specific hypotheses of whether one or more coefficients are significant, likelihood ratio test is a good choice. Its implementation in models with random effects have some caveats. Linear mixed models are usually fitted with the restricted maximum likelihood estimator, but this restricted maximum likelihood as the model results cannot be used in likelihood-ratio tests of fixed effects. We need to switch to maximum likelihood, usually by specifying the argument REML = FALSE or method = ML in the estimating function before testing with anova(). To test random effects, in contrast, we are usually interested in seeing if the standard deviation of the random term is zero or larger than zero. This test is on the boundary of the possible range of the parameter and the regular likelihood-ratio test such as anova() is inappropriate. This is because the test statistic under the null hypothesis $H_0: \sigma = 0$ is not $\chi^2(1)$ distributed but a mixture of two or more $\chi^2$ of different degrees of freedom, $\chi^2(0)/2 + \chi^2(1)/2$ in this one-variance case. It requires specific likelihood-ratio test methods for variance components, usually included in the package. See Molenberghs, G., & Verbeke, G. (2007). Likelihood ratio, score, and Wald tests in a constrained parameter space. The American Statistician, 61(1), 22–27. jstor.org/stable/27643833.

added 150 characters in body
Source Link
DrJerryTAO
DrJerryTAO
  • 2.3k
  • 6
  • 22

Notice that using a three-level factor gives coefficients as contrasts to L. Separating the factor into two binary indicators, however, often introduces in redundant terms if there is a missing group. Removing the intercept to combat multicollinearity often jeopardizejeopardizes some model summary statistics like $R^2$. Instead, specifying one main effect and one interaction M + L : M seem to offer the easiest interpretation among three groups where the coefficient of M gives its additionalhigher effect than L and the coefficient of L : M interaction gives the difference between ML and M. To find the difference between ML and L + M, use the difference between the coefficient of L : M interaction and the intercept.

Notice that using a three-level factor gives coefficients as contrasts to L. Separating the factor into two binary indicators, however, often introduces in redundant terms if there is a missing group. Removing the intercept to combat multicollinearity often jeopardize some model summary statistics like $R^2$. Instead, specifying one main effect and one interaction M + L : M seem to offer the easiest interpretation where the coefficient of M gives its additional effect than L and the coefficient of L : M interaction gives the difference between ML and M.

Notice that using a three-level factor gives coefficients as contrasts to L. Separating the factor into two binary indicators, however, often introduces in redundant terms if there is a missing group. Removing the intercept to combat multicollinearity often jeopardizes some model summary statistics like $R^2$. Instead, specifying one main effect and one interaction M + L : M seem to offer the easiest interpretation among three groups where the coefficient of M gives its higher effect than L and the coefficient of L : M interaction gives the difference between ML and M. To find the difference between ML and L + M, use the difference between the coefficient of L : M interaction and the intercept.

Source Link
DrJerryTAO
DrJerryTAO
  • 2.3k
  • 6
  • 22