Generando números aleatorios de un modelo lineal en R

Simulando 22 números con el siguiente modelo lineal: y = β₀ + β₁x + ε. En donde ε ~N(0,2²). Asumiendo x~N(0,1²),  y β₁ = 2

# Generando números para "x", "e", "y"
set.seed(4)
x <- rnorm(22)
e <- rnorm(22, 0, 2)
y <- 0.5 + 2*x + e

# Mostrando números generados
xy  = cbind(y, x)
xy

            y        x
 [1,]  3.5488  0.21675
 [2,]  1.9915 -0.54249
 [3,]  3.4681  0.89114
 [4,]  1.1261  0.59598
 [5,]  6.2830  1.63562
 [6,]  3.6982  0.68928
 [7,] -3.9185 -1.28125
 [8,]  2.5541 -0.21314
 [9,]  4.6000  1.89654
[10,]  6.1576  1.77686
[11,]  0.1248  0.56660
[12,] -2.4329  0.01572
[13,]  2.9884  0.38306
[14,] -0.3993 -0.04514
[15,]  0.1139  0.03435
[16,]  2.7062  0.16903
[17,]  1.8983  1.16503
[18,] -0.8635 -0.04420
[19,]  2.9867 -0.10037
[20,]  0.2962 -0.28344
[21,]  6.1667  1.54081
[22,] -2.5458  0.16517

# Estadística descriptiva de números generados
summary(xy)

       y                x          
 Min.   :-3.919   Min.   :-1.2812  
 1st Qu.: 0.117   1st Qu.:-0.0449  
 Median : 2.273   Median : 0.1929  
 Mean   : 1.843   Mean   : 0.4196  
 3rd Qu.: 3.529   3rd Qu.: 0.8407  
 Max.   : 6.283   Max.   : 1.8965  

# Trazando gráfico
plot(x,y)

hist(y)