Skip to main content
added 96 characters in body
Source Link
dipetkov
  • 10.7k
  • 2
  • 20
  • 55

Once you fit the model you can look at any number of comparisons between (functions of) the coefficients. For lm, glm, gls and geepack::geese linear models, you don't have to define the contrasts by hand; use the contrastcontrast package instead.

# Fit the model with the default treatment contrast. (Any contrast will do.)
fit <- lm(y ~ type, data = dat)

library("contrast")

# Example 1: Compare the control to each of the treatments.
print(
  contrast(
    fit,
    list(type = c("A", "B", "C")),
    list(type = "Contr")
  ),
  X = TRUE
)
#> lm model parameter contrast
#> 
#>  Contrast     S.E.    Lower    Upper     t df Pr(>|t|)
#>  18.61365 1.277391 16.02299 21.20432 14.57 36        0
#>  21.01564 1.277391 18.42497 23.60631 16.45 36        0
#>  15.11610 1.277391 12.52543 17.70677 11.83 36        0
#> 
#> Contrast coefficients:
#>  (Intercept) typeA typeB typeC
#>            0     1     0     0
#>            0     0     1     0
#>            0     0     0     1

# Example 2: Compare the control to the average of the treatments.
print(
  contrast(fit,
    list(type = c("A", "B", "C")),
    list(type = "Contr"),
    type = "average"
  ),
  X = TRUE
)
#> lm model parameter contrast
#> 
#>   Contrast     S.E.    Lower    Upper    t df Pr(>|t|)
#> 1 18.24847 1.042985 16.13319 20.36374 17.5 36        0
#> 
#> Contrast coefficients:
#>   (Intercept)     typeA     typeB     typeC
#> 1           0 0.3333333 0.3333333 0.3333333

# Example 3: Compare the control to a weighted average of the treatments.
print(
  contrast(fit,
    list(type = c("A", "B", "C")),
    list(type = "Contr"),
    weights = c(1 / 2, 1 / 4, 1 / 4),
    type = "average"
  ),
  X = TRUE
)
#> lm model parameter contrast
#> 
#>   Contrast     S.E.    Lower    Upper     t df Pr(>|t|)
#> 1 18.33976 1.059157 16.19169 20.48783 17.32 36        0
#> 
#> Contrast coefficients:
#>   (Intercept) typeA typeB typeC
#> 1           0   0.5  0.25  0.25

Once you fit the model you can look at any number of comparisons between (functions of) the coefficients. For lm, glm, gls and geepack::geese linear models, you don't have to define the contrasts by hand; use the contrast package instead.

# Fit the model with the default treatment contrast. (Any contrast will do.)
fit <- lm(y ~ type, data = dat)

library("contrast")

# Compare the control to each of the treatments.
print(
  contrast(
    fit,
    list(type = c("A", "B", "C")),
    list(type = "Contr")
  ),
  X = TRUE
)
#> lm model parameter contrast
#> 
#>  Contrast     S.E.    Lower    Upper     t df Pr(>|t|)
#>  18.61365 1.277391 16.02299 21.20432 14.57 36        0
#>  21.01564 1.277391 18.42497 23.60631 16.45 36        0
#>  15.11610 1.277391 12.52543 17.70677 11.83 36        0
#> 
#> Contrast coefficients:
#>  (Intercept) typeA typeB typeC
#>            0     1     0     0
#>            0     0     1     0
#>            0     0     0     1

# Compare the control to the average of the treatments.
print(
  contrast(fit,
    list(type = c("A", "B", "C")),
    list(type = "Contr"),
    type = "average"
  ),
  X = TRUE
)
#> lm model parameter contrast
#> 
#>   Contrast     S.E.    Lower    Upper    t df Pr(>|t|)
#> 1 18.24847 1.042985 16.13319 20.36374 17.5 36        0
#> 
#> Contrast coefficients:
#>   (Intercept)     typeA     typeB     typeC
#> 1           0 0.3333333 0.3333333 0.3333333

# Compare the control to a weighted average of the treatments.
print(
  contrast(fit,
    list(type = c("A", "B", "C")),
    list(type = "Contr"),
    weights = c(1 / 2, 1 / 4, 1 / 4),
    type = "average"
  ),
  X = TRUE
)
#> lm model parameter contrast
#> 
#>   Contrast     S.E.    Lower    Upper     t df Pr(>|t|)
#> 1 18.33976 1.059157 16.19169 20.48783 17.32 36        0
#> 
#> Contrast coefficients:
#>   (Intercept) typeA typeB typeC
#> 1           0   0.5  0.25  0.25

Once you fit the model you can look at any number of comparisons between (functions of) the coefficients. For lm, glm, gls and geepack::geese linear models, you don't have to define the contrasts by hand; use the contrast package instead.

# Fit the model with the default treatment contrast. (Any contrast will do.)
fit <- lm(y ~ type, data = dat)

library("contrast")

# Example 1: Compare the control to each of the treatments.
print(
  contrast(
    fit,
    list(type = c("A", "B", "C")),
    list(type = "Contr")
  ),
  X = TRUE
)
#> lm model parameter contrast
#> 
#>  Contrast     S.E.    Lower    Upper     t df Pr(>|t|)
#>  18.61365 1.277391 16.02299 21.20432 14.57 36        0
#>  21.01564 1.277391 18.42497 23.60631 16.45 36        0
#>  15.11610 1.277391 12.52543 17.70677 11.83 36        0
#> 
#> Contrast coefficients:
#>  (Intercept) typeA typeB typeC
#>            0     1     0     0
#>            0     0     1     0
#>            0     0     0     1

# Example 2: Compare the control to the average of the treatments.
print(
  contrast(fit,
    list(type = c("A", "B", "C")),
    list(type = "Contr"),
    type = "average"
  ),
  X = TRUE
)
#> lm model parameter contrast
#> 
#>   Contrast     S.E.    Lower    Upper    t df Pr(>|t|)
#> 1 18.24847 1.042985 16.13319 20.36374 17.5 36        0
#> 
#> Contrast coefficients:
#>   (Intercept)     typeA     typeB     typeC
#> 1           0 0.3333333 0.3333333 0.3333333

# Example 3: Compare the control to a weighted average of the treatments.
print(
  contrast(fit,
    list(type = c("A", "B", "C")),
    list(type = "Contr"),
    weights = c(1 / 2, 1 / 4, 1 / 4),
    type = "average"
  ),
  X = TRUE
)
#> lm model parameter contrast
#> 
#>   Contrast     S.E.    Lower    Upper     t df Pr(>|t|)
#> 1 18.33976 1.059157 16.19169 20.48783 17.32 36        0
#> 
#> Contrast coefficients:
#>   (Intercept) typeA typeB typeC
#> 1           0   0.5  0.25  0.25
added 616 characters in body
Source Link
dipetkov
  • 10.7k
  • 2
  • 20
  • 55
# Fit the model with the default treatment contrast. (Any contrast will do.)
fit <- lm(y ~ type, data = dat)

library("contrast")

# Compare the control to each of the treatments.
print(
  contrast(
    fit,
    list(type = c("A", "B", "C")),
    list(type = "Contr")
  ),
  X = TRUE
)
#> lm model parameter contrast
#> 
#>  Contrast     S.E.    Lower    Upper     t df Pr(>|t|)
#>  18.61365 1.277391 16.02299 21.20432 14.57 36        0
#>  21.01564 1.277391 18.42497 23.60631 16.45 36        0
#>  15.11610 1.277391 12.52543 17.70677 11.83 36        0
#> 
#> Contrast coefficients:
#>  (Intercept) typeA typeB typeC
#>            0     1     0     0
#>            0     0     1     0
#>            0     0     0     1

# Compare the control to the average of the treatments.
print(
  contrast(fit,
    list(type = c("A", "B", "C")),
    list(type = "Contr"),
    type = "average"
  ),
  X = TRUE
)
#> lm model parameter contrast
#> 
#>   Contrast     S.E.    Lower    Upper    t df Pr(>|t|)
#> 1 18.24847 1.042985 16.13319 20.36374 17.5 36        0
#> 
#> Contrast coefficients:
#>   (Intercept)     typeA     typeB     typeC
#> 1           0 0.3333333 0.3333333 0.3333333

# Compare the control to a weighted average of the treatments.
print(
  contrast(fit,
    list(type = c("A", "B", "C")),
    list(type = "Contr"),
    weights = c(1 / 2, 1 / 4, 1 / 4),
    type = "average"
  ),
  X = TRUE
)
#> lm model parameter contrast
#> 
#>   Contrast     S.E.    Lower    Upper     t df Pr(>|t|)
#> 1 18.33976 1.059157 16.19169 20.48783 17.32 36        0
#> 
#> Contrast coefficients:
#>   (Intercept) typeA typeB typeC
#> 1           0   0.5  0.25  0.25
```

As the examples show, a contrast is a linear combination of the factor levels; it can be represented as a row vector. If you can extract the coefficients and their covariance matrix from the fitted model, you can find the contrast and its standard error by hand.

# Compute by hand the contrast that compares the control 
# to the average of the treatments.
X <- matrix(c(0, 1 / 3, 1 / 3, 1 / 3), nrow = 1)

# Find the value of the contrast
X %*% fit$coefficients
#>          [,1]
#> [1,] 18.24847
# ... and its standard error.
sqrt(X %*% vcov(fit) %*% t(X))
#>          [,1]
#> [1,] 1.042985
# Fit the model with the default treatment contrast. (Any contrast will do.)
fit <- lm(y ~ type, data = dat)

library("contrast")

# Compare the control to each of the treatments.
print(
  contrast(
    fit,
    list(type = c("A", "B", "C")),
    list(type = "Contr")
  ),
  X = TRUE
)
#> lm model parameter contrast
#> 
#>  Contrast     S.E.    Lower    Upper     t df Pr(>|t|)
#>  18.61365 1.277391 16.02299 21.20432 14.57 36        0
#>  21.01564 1.277391 18.42497 23.60631 16.45 36        0
#>  15.11610 1.277391 12.52543 17.70677 11.83 36        0
#> 
#> Contrast coefficients:
#>  (Intercept) typeA typeB typeC
#>            0     1     0     0
#>            0     0     1     0
#>            0     0     0     1

# Compare the control to the average of the treatments.
print(
  contrast(fit,
    list(type = c("A", "B", "C")),
    list(type = "Contr"),
    type = "average"
  ),
  X = TRUE
)
#> lm model parameter contrast
#> 
#>   Contrast     S.E.    Lower    Upper    t df Pr(>|t|)
#> 1 18.24847 1.042985 16.13319 20.36374 17.5 36        0
#> 
#> Contrast coefficients:
#>   (Intercept)     typeA     typeB     typeC
#> 1           0 0.3333333 0.3333333 0.3333333

# Compare the control to a weighted average of the treatments.
print(
  contrast(fit,
    list(type = c("A", "B", "C")),
    list(type = "Contr"),
    weights = c(1 / 2, 1 / 4, 1 / 4),
    type = "average"
  ),
  X = TRUE
)
#> lm model parameter contrast
#> 
#>   Contrast     S.E.    Lower    Upper     t df Pr(>|t|)
#> 1 18.33976 1.059157 16.19169 20.48783 17.32 36        0
#> 
#> Contrast coefficients:
#>   (Intercept) typeA typeB typeC
#> 1           0   0.5  0.25  0.25
```
# Fit the model with the default treatment contrast. (Any contrast will do.)
fit <- lm(y ~ type, data = dat)

library("contrast")

# Compare the control to each of the treatments.
print(
  contrast(
    fit,
    list(type = c("A", "B", "C")),
    list(type = "Contr")
  ),
  X = TRUE
)
#> lm model parameter contrast
#> 
#>  Contrast     S.E.    Lower    Upper     t df Pr(>|t|)
#>  18.61365 1.277391 16.02299 21.20432 14.57 36        0
#>  21.01564 1.277391 18.42497 23.60631 16.45 36        0
#>  15.11610 1.277391 12.52543 17.70677 11.83 36        0
#> 
#> Contrast coefficients:
#>  (Intercept) typeA typeB typeC
#>            0     1     0     0
#>            0     0     1     0
#>            0     0     0     1

# Compare the control to the average of the treatments.
print(
  contrast(fit,
    list(type = c("A", "B", "C")),
    list(type = "Contr"),
    type = "average"
  ),
  X = TRUE
)
#> lm model parameter contrast
#> 
#>   Contrast     S.E.    Lower    Upper    t df Pr(>|t|)
#> 1 18.24847 1.042985 16.13319 20.36374 17.5 36        0
#> 
#> Contrast coefficients:
#>   (Intercept)     typeA     typeB     typeC
#> 1           0 0.3333333 0.3333333 0.3333333

# Compare the control to a weighted average of the treatments.
print(
  contrast(fit,
    list(type = c("A", "B", "C")),
    list(type = "Contr"),
    weights = c(1 / 2, 1 / 4, 1 / 4),
    type = "average"
  ),
  X = TRUE
)
#> lm model parameter contrast
#> 
#>   Contrast     S.E.    Lower    Upper     t df Pr(>|t|)
#> 1 18.33976 1.059157 16.19169 20.48783 17.32 36        0
#> 
#> Contrast coefficients:
#>   (Intercept) typeA typeB typeC
#> 1           0   0.5  0.25  0.25

As the examples show, a contrast is a linear combination of the factor levels; it can be represented as a row vector. If you can extract the coefficients and their covariance matrix from the fitted model, you can find the contrast and its standard error by hand.

# Compute by hand the contrast that compares the control 
# to the average of the treatments.
X <- matrix(c(0, 1 / 3, 1 / 3, 1 / 3), nrow = 1)

# Find the value of the contrast
X %*% fit$coefficients
#>          [,1]
#> [1,] 18.24847
# ... and its standard error.
sqrt(X %*% vcov(fit) %*% t(X))
#>          [,1]
#> [1,] 1.042985
added 56 characters in body
Source Link
dipetkov
  • 10.7k
  • 2
  • 20
  • 55

You are trying to define two contrasts with one matrix. To fit thea model you have to choose one contrast, so that lmthe model understands how to interpret the coefficients and estimate them.

Once you fit the model you can look at any number of comparisons between (functions of) the coefficients. AndFor lm, glm, gls and geepack::geese linear models, you don't have to define the contrasts by hand; use the contrastcontrast package instead.

You are trying to define two contrasts with one matrix. To fit the model you have to choose one contrast, so that lm understands how to interpret the coefficients and estimate them.

Once you fit the model you can look at any number of comparisons. And you don't have to define the contrasts by hand; use the contrast package instead.

You are trying to define two contrasts with one matrix. To fit a model you have to choose one contrast, so that the model understands how to interpret the coefficients and estimate them.

Once you fit the model you can look at any number of comparisons between (functions of) the coefficients. For lm, glm, gls and geepack::geese linear models, you don't have to define the contrasts by hand; use the contrast package instead.

Source Link
dipetkov
  • 10.7k
  • 2
  • 20
  • 55
Loading