Latex-like codes were not processed (Greek letters etc. were not shown) in plot in R

(reproducible code given) I am studying Ugarte2016's "Probability and Statistics with R" 2E. The following code is run in R but Latex-like code is not processed. It seems that the code inside "$...$" is not processed. The code supplied below was from the authors of the book. There seems a problem somehow. What could be the problem?

######### Chapter 12 #############
library(PASWR2); library(ggplot2); library(car); library(scatterplot3d)
library(gridExtra); library(multcomp); library(leaps); library(MASS)
################ Figure 12.1 ###############    
opar <- par(no.readonly = TRUE) # copy of current settings
par(mar=c(2, 14, 2, 1), las = 1)
DF <- data.frame(x = c(1, 4, 9), y = c(1, 4, 9))
plot(y~x, data = DF, xaxt = "n", yaxt = "n", xlim = c(0, 12), ylim = c(-2, 12),  xlab = "", ylab = "", type = "n")
abline(lm(y~x, data = DF), lwd = 2)
axis(side =1, at =c(1, 4, 10), labels = c("$x_1$", "$x_2$", "$x_3$"))
axis(side =2, at =c(1, 4, 10),  labels = c("$E(Y|x_1) = \\beta_0 + \\beta_1x_1$",  "$E(Y|x_1) = \\beta_0 + \\beta_1x_1$", "$E(Y|x_1) = \\beta_0 + \\beta_1x_1$") )
segments(1, -2, 1, 2.5, lty = "dashed")
segments(0, 1, 1 + 0.75, 1, lty = "dashed")
segments(4, -2, 4, 5.5, lty = "dashed")
segments(0, 4, 4 + 0.75, 4, lty = "dashed")
segments(10, -2, 10, 11.5, lty = "dashed")
segments(0, 10, 10 + 0.75, 10, lty = "dashed")
ys <- seq(-1.5, 1.5, length = 200)
xs <- dnorm(ys, 0, 0.5)
lines(xs + 1, ys + 1, type = "l",lwd = 2)
lines(xs + 4, ys + 4, type = "l",lwd = 2)
lines(xs + 10, ys + 10, type = "l",lwd = 2)
text(7.8, 5.5, "$E(Y|x) = \\beta_0 + \\beta_1x$")
arrows(8, 5.7, 7, 7, length = 0.1, lwd = 2)
par(opar)

The code result:

The image in the book:

Solution

Use package latex2exp:

######### Chapter 12 #############
library(PASWR2); library(ggplot2); library(car); library(scatterplot3d)
library(gridExtra); library(multcomp); library(leaps); library(MASS)
library(latex2exp)
################ Figure 12.1 ###############    
opar <- par(no.readonly = TRUE) # copy of current settings
par(mar=c(2, 14, 2, 1), las = 1)
DF <- data.frame(x = c(1, 4, 9), y = c(1, 4, 9))
plot(y~x, data = DF, xaxt = "n", yaxt = "n", xlim = c(0, 12), ylim = c(-2, 12),  xlab = "", ylab = "", type = "n")
abline(lm(y~x, data = DF), lwd = 2)
axis(side =1, at =c(1, 4, 10), labels = TeX(c("$x_1$", "$x_2$", "$x_3$")))
axis(side =2, at =c(1, 4, 10),  labels = TeX(c("$E(Y|x_1) = \\beta_0 + \\beta_1x_1$",  "$E(Y|x_1) = \\beta_0 + \\beta_1x_1$", "$E(Y|x_1) = \\beta_0 + \\beta_1x_1$") ))
segments(1, -2, 1, 2.5, lty = "dashed")
segments(0, 1, 1 + 0.75, 1, lty = "dashed")
segments(4, -2, 4, 5.5, lty = "dashed")
segments(0, 4, 4 + 0.75, 4, lty = "dashed")
segments(10, -2, 10, 11.5, lty = "dashed")
segments(0, 10, 10 + 0.75, 10, lty = "dashed")
ys <- seq(-1.5, 1.5, length = 200)
xs <- dnorm(ys, 0, 0.5)
lines(xs + 1, ys + 1, type = "l",lwd = 2)
lines(xs + 4, ys + 4, type = "l",lwd = 2)
lines(xs + 10, ys + 10, type = "l",lwd = 2)
text(7.8, 5.5, TeX("$E(Y|x) = \\beta_0 + \\beta_1x$"))
arrows(8, 5.7, 7, 7, length = 0.1, lwd = 2)
par(opar)