Shows what matrices \(\\mathbf{A}, \\mathbf{b}\) look like as the system of linear equations, \(\\mathbf{A x} = \\mathbf{b}\), but written out as a set of equations.
Arguments
- A
either the matrix of coefficients of a system of linear equations, or the matrix
cbind(A,b). The matrix can be numeric or character. Alternatively, can be of class'lm'to print the equations for the design matrix in a linear regression model- b
if supplied, the vector of constants on the right hand side of the equations. When omitted the values
b1, b2, ..., bnwill be used as placeholders- vars
a numeric or character vector of names of the variables. If supplied, the length must be equal to the number of unknowns in the equations. The default is
paste0("x", 1:ncol(A).- simplify
logical; try to simplify the equations?
- reduce
logical; only show the unique linear equations
- fractions
logical; express numbers as rational fractions, using the
fractionsfunction; if you require greater accuracy, you can set thecycles(default 10) and/ormax.denominator(default 2000) arguments tofractionsas a global option, e.g.,options(fractions=list(cycles=100, max.denominator=10^4)).- latex
logical; print equations in a form suitable for LaTeX output?
References
Fox, J. and Friendly, M. (2016). "Visualizing Simultaneous Linear Equations, Geometric Vectors, and Least-Squares Regression with the matlib Package for R". useR Conference, Stanford, CA, June 27 - June 30, 2016.
Examples
A <- matrix(c(2, 1, -1,
-3, -1, 2,
-2, 1, 2), 3, 3, byrow=TRUE)
b <- c(8, -11, -3)
showEqn(A, b)
#> 2*x1 + 1*x2 - 1*x3 = 8
#> -3*x1 - 1*x2 + 2*x3 = -11
#> -2*x1 + 1*x2 + 2*x3 = -3
# show numerically
x <- solve(A, b)
showEqn(A, b, vars=x)
#> 2*2 + 1*3 - 1*-1 = 8
#> -3*2 - 1*3 + 2*-1 = -11
#> -2*2 + 1*3 + 2*-1 = -3
showEqn(A, b, simplify=TRUE)
#> 2*x1 + x2 - 1*x3 = 8
#> -3*x1 - 1*x2 + 2*x3 = -11
#> -2*x1 + x2 + 2*x3 = -3
showEqn(A, b, latex=TRUE)
#> \begin{array}{lllllll}
#> 2 \cdot x_1 &+& 1 \cdot x_2 &-& 1 \cdot x_3 &=& 8 \\
#> -3 \cdot x_1 &-& 1 \cdot x_2 &+& 2 \cdot x_3 &=& -11 \\
#> -2 \cdot x_1 &+& 1 \cdot x_2 &+& 2 \cdot x_3 &=& -3 \\
#> \end{array}
# lower triangle of equation with zeros omitted (for back solving)
A <- matrix(c(2, 1, 2,
-3, -1, 2,
-2, 1, 2), 3, 3, byrow=TRUE)
U <- LU(A)$U
showEqn(U, simplify=TRUE, fractions=TRUE)
#> 2*x1 + x2 + 2*x3 = b1
#> 1/2*x2 + 5*x3 = b2
#> - 16*x3 = b3
showEqn(U, b, simplify=TRUE, fractions=TRUE)
#> 2*x1 + x2 + 2*x3 = 8
#> 1/2*x2 + 5*x3 = -11
#> - 16*x3 = -3
####################
# Linear models Design Matricies
data(mtcars)
ancova <- lm(mpg ~ wt + vs, mtcars)
summary(ancova)
#>
#> Call:
#> lm(formula = mpg ~ wt + vs, data = mtcars)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -3.7071 -2.4415 -0.3129 1.4319 6.0156
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 33.0042 2.3554 14.012 1.92e-14 ***
#> wt -4.4428 0.6134 -7.243 5.63e-08 ***
#> vs 3.1544 1.1907 2.649 0.0129 *
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 2.78 on 29 degrees of freedom
#> Multiple R-squared: 0.801, Adjusted R-squared: 0.7873
#> F-statistic: 58.36 on 2 and 29 DF, p-value: 6.818e-11
#>
showEqn(ancova)
#> 1*x1 + 2.62*x2 + 0*x3 = b1
#> 1*x1 + 2.875*x2 + 0*x3 = b2
#> 1*x1 + 2.32*x2 + 1*x3 = b3
#> 1*x1 + 3.215*x2 + 1*x3 = b4
#> 1*x1 + 3.44*x2 + 0*x3 = b5
#> 1*x1 + 3.46*x2 + 1*x3 = b6
#> 1*x1 + 3.57*x2 + 0*x3 = b7
#> 1*x1 + 3.19*x2 + 1*x3 = b8
#> 1*x1 + 3.15*x2 + 1*x3 = b9
#> 1*x1 + 3.44*x2 + 1*x3 = b10
#> 1*x1 + 3.44*x2 + 1*x3 = b11
#> 1*x1 + 4.07*x2 + 0*x3 = b12
#> 1*x1 + 3.73*x2 + 0*x3 = b13
#> 1*x1 + 3.78*x2 + 0*x3 = b14
#> 1*x1 + 5.25*x2 + 0*x3 = b15
#> 1*x1 + 5.424*x2 + 0*x3 = b16
#> 1*x1 + 5.345*x2 + 0*x3 = b17
#> 1*x1 + 2.2*x2 + 1*x3 = b18
#> 1*x1 + 1.615*x2 + 1*x3 = b19
#> 1*x1 + 1.835*x2 + 1*x3 = b20
#> 1*x1 + 2.465*x2 + 1*x3 = b21
#> 1*x1 + 3.52*x2 + 0*x3 = b22
#> 1*x1 + 3.435*x2 + 0*x3 = b23
#> 1*x1 + 3.84*x2 + 0*x3 = b24
#> 1*x1 + 3.845*x2 + 0*x3 = b25
#> 1*x1 + 1.935*x2 + 1*x3 = b26
#> 1*x1 + 2.14*x2 + 0*x3 = b27
#> 1*x1 + 1.513*x2 + 1*x3 = b28
#> 1*x1 + 3.17*x2 + 0*x3 = b29
#> 1*x1 + 2.77*x2 + 0*x3 = b30
#> 1*x1 + 3.57*x2 + 0*x3 = b31
#> 1*x1 + 2.78*x2 + 1*x3 = b32
showEqn(ancova, simplify=TRUE)
#> x1 + 2.62*x2 = b1
#> x1 + 2.875*x2 = b2
#> x1 + 2.32*x2 + x3 = b3
#> x1 + 3.215*x2 + x3 = b4
#> x1 + 3.44*x2 = b5
#> x1 + 3.46*x2 + x3 = b6
#> x1 + 3.57*x2 = b7
#> x1 + 3.19*x2 + x3 = b8
#> x1 + 3.15*x2 + x3 = b9
#> x1 + 3.44*x2 + x3 = b10
#> x1 + 3.44*x2 + x3 = b11
#> x1 + 4.07*x2 = b12
#> x1 + 3.73*x2 = b13
#> x1 + 3.78*x2 = b14
#> x1 + 5.25*x2 = b15
#> x1 + 5.424*x2 = b16
#> x1 + 5.345*x2 = b17
#> x1 + 2.2*x2 + x3 = b18
#> x1 + 1.615*x2 + x3 = b19
#> x1 + 1.835*x2 + x3 = b20
#> x1 + 2.465*x2 + x3 = b21
#> x1 + 3.52*x2 = b22
#> x1 + 3.435*x2 = b23
#> x1 + 3.84*x2 = b24
#> x1 + 3.845*x2 = b25
#> x1 + 1.935*x2 + x3 = b26
#> x1 + 2.14*x2 = b27
#> x1 + 1.513*x2 + x3 = b28
#> x1 + 3.17*x2 = b29
#> x1 + 2.77*x2 = b30
#> x1 + 3.57*x2 = b31
#> x1 + 2.78*x2 + x3 = b32
showEqn(ancova, vars=round(coef(ancova),2))
#> 1*33 + 2.62*-4.44 + 0*3.15 = b1
#> 1*33 + 2.875*-4.44 + 0*3.15 = b2
#> 1*33 + 2.32*-4.44 + 1*3.15 = b3
#> 1*33 + 3.215*-4.44 + 1*3.15 = b4
#> 1*33 + 3.44*-4.44 + 0*3.15 = b5
#> 1*33 + 3.46*-4.44 + 1*3.15 = b6
#> 1*33 + 3.57*-4.44 + 0*3.15 = b7
#> 1*33 + 3.19*-4.44 + 1*3.15 = b8
#> 1*33 + 3.15*-4.44 + 1*3.15 = b9
#> 1*33 + 3.44*-4.44 + 1*3.15 = b10
#> 1*33 + 3.44*-4.44 + 1*3.15 = b11
#> 1*33 + 4.07*-4.44 + 0*3.15 = b12
#> 1*33 + 3.73*-4.44 + 0*3.15 = b13
#> 1*33 + 3.78*-4.44 + 0*3.15 = b14
#> 1*33 + 5.25*-4.44 + 0*3.15 = b15
#> 1*33 + 5.424*-4.44 + 0*3.15 = b16
#> 1*33 + 5.345*-4.44 + 0*3.15 = b17
#> 1*33 + 2.2*-4.44 + 1*3.15 = b18
#> 1*33 + 1.615*-4.44 + 1*3.15 = b19
#> 1*33 + 1.835*-4.44 + 1*3.15 = b20
#> 1*33 + 2.465*-4.44 + 1*3.15 = b21
#> 1*33 + 3.52*-4.44 + 0*3.15 = b22
#> 1*33 + 3.435*-4.44 + 0*3.15 = b23
#> 1*33 + 3.84*-4.44 + 0*3.15 = b24
#> 1*33 + 3.845*-4.44 + 0*3.15 = b25
#> 1*33 + 1.935*-4.44 + 1*3.15 = b26
#> 1*33 + 2.14*-4.44 + 0*3.15 = b27
#> 1*33 + 1.513*-4.44 + 1*3.15 = b28
#> 1*33 + 3.17*-4.44 + 0*3.15 = b29
#> 1*33 + 2.77*-4.44 + 0*3.15 = b30
#> 1*33 + 3.57*-4.44 + 0*3.15 = b31
#> 1*33 + 2.78*-4.44 + 1*3.15 = b32
showEqn(ancova, vars=round(coef(ancova),2), simplify=TRUE)
#> 33 + 2.62*-4.44 = b1
#> 33 + 2.875*-4.44 = b2
#> 33 + 2.32*-4.44 + 3.15 = b3
#> 33 + 3.215*-4.44 + 3.15 = b4
#> 33 + 3.44*-4.44 = b5
#> 33 + 3.46*-4.44 + 3.15 = b6
#> 33 + 3.57*-4.44 = b7
#> 33 + 3.19*-4.44 + 3.15 = b8
#> 33 + 3.15*-4.44 + 3.15 = b9
#> 33 + 3.44*-4.44 + 3.15 = b10
#> 33 + 3.44*-4.44 + 3.15 = b11
#> 33 + 4.07*-4.44 = b12
#> 33 + 3.73*-4.44 = b13
#> 33 + 3.78*-4.44 = b14
#> 33 + 5.25*-4.44 = b15
#> 33 + 5.424*-4.44 = b16
#> 33 + 5.345*-4.44 = b17
#> 33 + 2.2*-4.44 + 3.15 = b18
#> 33 + 1.615*-4.44 + 3.15 = b19
#> 33 + 1.835*-4.44 + 3.15 = b20
#> 33 + 2.465*-4.44 + 3.15 = b21
#> 33 + 3.52*-4.44 = b22
#> 33 + 3.435*-4.44 = b23
#> 33 + 3.84*-4.44 = b24
#> 33 + 3.845*-4.44 = b25
#> 33 + 1.935*-4.44 + 3.15 = b26
#> 33 + 2.14*-4.44 = b27
#> 33 + 1.513*-4.44 + 3.15 = b28
#> 33 + 3.17*-4.44 = b29
#> 33 + 2.77*-4.44 = b30
#> 33 + 3.57*-4.44 = b31
#> 33 + 2.78*-4.44 + 3.15 = b32
twoway_int <- lm(mpg ~ vs * am, mtcars)
summary(twoway_int)
#>
#> Call:
#> lm(formula = mpg ~ vs * am, data = mtcars)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -6.971 -1.973 0.300 2.036 6.250
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 15.050 1.002 15.017 6.34e-15 ***
#> vs 5.693 1.651 3.448 0.0018 **
#> am 4.700 1.736 2.708 0.0114 *
#> vs:am 2.929 2.541 1.153 0.2589
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 3.472 on 28 degrees of freedom
#> Multiple R-squared: 0.7003, Adjusted R-squared: 0.6682
#> F-statistic: 21.81 on 3 and 28 DF, p-value: 1.735e-07
#>
car::Anova(twoway_int)
#> Anova Table (Type II tests)
#>
#> Response: mpg
#> Sum Sq Df F value Pr(>F)
#> vs 367.41 1 30.4836 6.687e-06 ***
#> am 276.03 1 22.9021 4.984e-05 ***
#> vs:am 16.01 1 1.3283 0.2589
#> Residuals 337.48 28
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
showEqn(twoway_int)
#> 1*x1 + 0*x2 + 1*x3 + 0*x4 = b1
#> 1*x1 + 0*x2 + 1*x3 + 0*x4 = b2
#> 1*x1 + 1*x2 + 1*x3 + 1*x4 = b3
#> 1*x1 + 1*x2 + 0*x3 + 0*x4 = b4
#> 1*x1 + 0*x2 + 0*x3 + 0*x4 = b5
#> 1*x1 + 1*x2 + 0*x3 + 0*x4 = b6
#> 1*x1 + 0*x2 + 0*x3 + 0*x4 = b7
#> 1*x1 + 1*x2 + 0*x3 + 0*x4 = b8
#> 1*x1 + 1*x2 + 0*x3 + 0*x4 = b9
#> 1*x1 + 1*x2 + 0*x3 + 0*x4 = b10
#> 1*x1 + 1*x2 + 0*x3 + 0*x4 = b11
#> 1*x1 + 0*x2 + 0*x3 + 0*x4 = b12
#> 1*x1 + 0*x2 + 0*x3 + 0*x4 = b13
#> 1*x1 + 0*x2 + 0*x3 + 0*x4 = b14
#> 1*x1 + 0*x2 + 0*x3 + 0*x4 = b15
#> 1*x1 + 0*x2 + 0*x3 + 0*x4 = b16
#> 1*x1 + 0*x2 + 0*x3 + 0*x4 = b17
#> 1*x1 + 1*x2 + 1*x3 + 1*x4 = b18
#> 1*x1 + 1*x2 + 1*x3 + 1*x4 = b19
#> 1*x1 + 1*x2 + 1*x3 + 1*x4 = b20
#> 1*x1 + 1*x2 + 0*x3 + 0*x4 = b21
#> 1*x1 + 0*x2 + 0*x3 + 0*x4 = b22
#> 1*x1 + 0*x2 + 0*x3 + 0*x4 = b23
#> 1*x1 + 0*x2 + 0*x3 + 0*x4 = b24
#> 1*x1 + 0*x2 + 0*x3 + 0*x4 = b25
#> 1*x1 + 1*x2 + 1*x3 + 1*x4 = b26
#> 1*x1 + 0*x2 + 1*x3 + 0*x4 = b27
#> 1*x1 + 1*x2 + 1*x3 + 1*x4 = b28
#> 1*x1 + 0*x2 + 1*x3 + 0*x4 = b29
#> 1*x1 + 0*x2 + 1*x3 + 0*x4 = b30
#> 1*x1 + 0*x2 + 1*x3 + 0*x4 = b31
#> 1*x1 + 1*x2 + 1*x3 + 1*x4 = b32
showEqn(twoway_int, reduce=TRUE)
#> 1*x1 + 0*x2 + 1*x3 + 0*x4 = b1
#> 1*x1 + 1*x2 + 1*x3 + 1*x4 = b2
#> 1*x1 + 1*x2 + 0*x3 + 0*x4 = b3
#> 1*x1 + 0*x2 + 0*x3 + 0*x4 = b4
showEqn(twoway_int, reduce=TRUE, simplify=TRUE)
#> x1 + x3 = b1
#> x1 + x2 + x3 + x4 = b2
#> x1 + x2 = b3
#> x1 = b4
# Piece-wise linear regression
x <- c(1:10, 13:22)
y <- numeric(20)
y[1:10] <- 20:11 + rnorm(10, 0, 1.5)
y[11:20] <- seq(11, 15, len=10) + rnorm(10, 0, 1.5)
plot(x, y, pch = 16)
x2 <- as.numeric(x > 10)
mod <- lm(y ~ x + I((x - 10) * x2))
summary(mod)
#>
#> Call:
#> lm(formula = y ~ x + I((x - 10) * x2))
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -2.4470 -0.6140 -0.4304 1.0993 2.4990
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 20.9943 0.9522 22.048 6.05e-14 ***
#> x -0.9731 0.1412 -6.891 2.61e-06 ***
#> I((x - 10) * x2) 1.3166 0.2203 5.977 1.50e-05 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 1.45 on 17 degrees of freedom
#> Multiple R-squared: 0.7444, Adjusted R-squared: 0.7143
#> F-statistic: 24.75 on 2 and 17 DF, p-value: 9.223e-06
#>
lines(x, fitted(mod))
showEqn(mod)
#> 1*x1 + 1*x2 + 0*x3 = b1
#> 1*x1 + 2*x2 + 0*x3 = b2
#> 1*x1 + 3*x2 + 0*x3 = b3
#> 1*x1 + 4*x2 + 0*x3 = b4
#> 1*x1 + 5*x2 + 0*x3 = b5
#> 1*x1 + 6*x2 + 0*x3 = b6
#> 1*x1 + 7*x2 + 0*x3 = b7
#> 1*x1 + 8*x2 + 0*x3 = b8
#> 1*x1 + 9*x2 + 0*x3 = b9
#> 1*x1 + 10*x2 + 0*x3 = b10
#> 1*x1 + 13*x2 + 3*x3 = b11
#> 1*x1 + 14*x2 + 4*x3 = b12
#> 1*x1 + 15*x2 + 5*x3 = b13
#> 1*x1 + 16*x2 + 6*x3 = b14
#> 1*x1 + 17*x2 + 7*x3 = b15
#> 1*x1 + 18*x2 + 8*x3 = b16
#> 1*x1 + 19*x2 + 9*x3 = b17
#> 1*x1 + 20*x2 + 10*x3 = b18
#> 1*x1 + 21*x2 + 11*x3 = b19
#> 1*x1 + 22*x2 + 12*x3 = b20
showEqn(mod, vars=round(coef(mod),2))
#> 1*20.99 + 1*-0.97 + 0*1.32 = b1
#> 1*20.99 + 2*-0.97 + 0*1.32 = b2
#> 1*20.99 + 3*-0.97 + 0*1.32 = b3
#> 1*20.99 + 4*-0.97 + 0*1.32 = b4
#> 1*20.99 + 5*-0.97 + 0*1.32 = b5
#> 1*20.99 + 6*-0.97 + 0*1.32 = b6
#> 1*20.99 + 7*-0.97 + 0*1.32 = b7
#> 1*20.99 + 8*-0.97 + 0*1.32 = b8
#> 1*20.99 + 9*-0.97 + 0*1.32 = b9
#> 1*20.99 + 10*-0.97 + 0*1.32 = b10
#> 1*20.99 + 13*-0.97 + 3*1.32 = b11
#> 1*20.99 + 14*-0.97 + 4*1.32 = b12
#> 1*20.99 + 15*-0.97 + 5*1.32 = b13
#> 1*20.99 + 16*-0.97 + 6*1.32 = b14
#> 1*20.99 + 17*-0.97 + 7*1.32 = b15
#> 1*20.99 + 18*-0.97 + 8*1.32 = b16
#> 1*20.99 + 19*-0.97 + 9*1.32 = b17
#> 1*20.99 + 20*-0.97 + 10*1.32 = b18
#> 1*20.99 + 21*-0.97 + 11*1.32 = b19
#> 1*20.99 + 22*-0.97 + 12*1.32 = b20
showEqn(mod, simplify=TRUE)
#> x1 + x2 = b1
#> x1 + 2*x2 = b2
#> x1 + 3*x2 = b3
#> x1 + 4*x2 = b4
#> x1 + 5*x2 = b5
#> x1 + 6*x2 = b6
#> x1 + 7*x2 = b7
#> x1 + 8*x2 = b8
#> x1 + 9*x2 = b9
#> x1 + 10*x2 = b10
#> x1 + 13*x2 + 3*x3 = b11
#> x1 + 14*x2 + 4*x3 = b12
#> x1 + 15*x2 + 5*x3 = b13
#> x1 + 16*x2 + 6*x3 = b14
#> x1 + 17*x2 + 7*x3 = b15
#> x1 + 18*x2 + 8*x3 = b16
#> x1 + 19*x2 + 9*x3 = b17
#> x1 + 20*x2 + 10*x3 = b18
#> x1 + 21*x2 + 11*x3 = b19
#> x1 + 22*x2 + 12*x3 = b20