Hi ,
I want to recreate this SAS plot in R:
On this plot there are 95% CI for fitted regression line and 95% Prediction Intervals for income variable - educ_12 variable was centered.
This is my code:
dat1 <- structure(list(educ = c(20, 20, 20, 20, 20, 20, 20, 20, 20, 20,
20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 19,
19, 19, 19, 19, 19, 19, 19, 19, 19, 19, 19, 19, 19, 19, 19, 19,
18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18,
18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18,
18, 18, 18, 18, 18, 18, 18, 17, 17, 17, 17, 17, 17, 17, 17, 17,
17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 17,
17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 16, 16, 16,
16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16,
16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16,
16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16,
16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16,
16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16,
16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16,
16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 15, 15, 15,
15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15,
15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15,
15, 15, 15, 15, 15, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14,
14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14,
14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14,
14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14,
14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14,
14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14,
14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 13, 13, 13,
13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13,
13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13,
13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13,
13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13,
13, 13, 13, 13, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12,
12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12,
12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12,
12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12,
12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12,
12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12,
12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12,
12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12,
12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12,
12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12,
12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12,
12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12,
12, 12, 12, 12, 12, 12, 12, 12, 12, 11, 11, 11, 11, 11, 11, 11,
11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11,
11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11,
11, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 9, 9, 9, 9, 9, 9,
9, 9, 9, 9, 9, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 6, 6, 6,
6, 6, 6, 6, 6, 6, 6, 6, 5, 4, 4, 4, 4, 2), educ_12 = c(8, 8,
8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
8, 8, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 6, 6,
6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,
6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 5, 5, 5, 5, 5,
5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 4, 4, 4, 4, 4, 4, 4, 4, 4,
4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 3, 3,
3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 2, 2, 2, 2,
2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, -1, -1, -1, -1, -1, -1,
-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
-1, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -3, -3, -3, -3, -3,
-3, -3, -3, -3, -3, -3, -4, -4, -4, -4, -4, -4, -4, -4, -4, -4,
-4, -4, -4, -6, -6, -6, -6, -6, -6, -6, -6, -6, -6, -6, -7, -8,
-8, -8, -8, -10), income = c(26.95, 15.925, 58.8, 33.075, 33.075,
26.95, 22.05, 58.8, 15.925, 5.5125, 10.4125, 18.375, 49, 7.9625,
49, 18.375, 3.185, 3.675, 33.075, 33.075, 40.425, 22.05, 15.925,
26.95, 9.1875, 22.05, 40.425, 18.375, 22.05, 15.925, 26.95, 40.425,
26.95, 22.05, 22.05, 33.075, 22.05, 0.98, 33.075, 13.475, 22.05,
33.075, 22.05, 58.8, 40.425, 22.05, 26.95, 49, 26.95, 2.205,
49, 22.05, 13.475, 49, 33.075, 15.925, 33.075, 33.075, 26.95,
40.425, 49, 26.95, 22.05, 26.95, 33.075, 40.425, 40.425, 33.075,
22.05, 15.925, 22.05, 26.95, 22.05, 58.8, 58.8, 22.05, 6.7375,
26.95, 10.4125, 7.9625, 26.95, 15.925, 22.05, 7.9625, 40.425,
33.075, 49, 40.425, 26.95, 49, 0.245, 40.425, 40.425, 58.8, 33.075,
10.4125, 22.05, 18.375, 0.98, 5.5125, 13.475, 33.075, 33.075,
33.075, 26.95, 49, 26.95, 13.475, 15.925, 11.6375, 0.98, 9.1875,
22.05, 11.6375, 9.1875, 33.075, 18.375, 26.95, 10.4125, 3.185,
49, 49, 33.075, 40.425, 33.075, 49, 33.075, 3.675, 40.425, 40.425,
26.95, 33.075, 58.8, 2.205, 5.5125, 10.4125, 26.95, 33.075, 22.05,
49, 26.95, 49, 49, 33.075, 26.95, 1.715, 22.05, 18.375, 9.1875,
58.8, 13.475, 33.075, 15.925, 68.6, 26.95, 15.925, 7.9625, 6.7375,
9.1875, 22.05, 40.425, 18.375, 26.95, 2.695, 40.425, 58.8, 18.375,
0.245, 2.695, 11.6375, 2.695, 58.8, 33.075, 18.375, 9.1875, 10.4125,
10.4125, 49, 1.715, 33.075, 58.8, 49, 26.95, 49, 68.6, 6.7375,
3.675, 18.375, 15.925, 3.185, 18.375, 9.1875, 4.41, 22.05, 68.6,
5.5125, 18.375, 40.425, 33.075, 9.1875, 7.9625, 18.375, 58.8,
5.5125, 7.9625, 26.95, 49, 26.95, 22.05, 15.925, 11.6375, 5.5125,
22.05, 7.9625, 49, 33.075, 0.98, 3.675, 49, 33.075, 7.9625, 15.925,
13.475, 40.425, 6.7375, 26.95, 6.7375, 49, 11.6375, 22.05, 6.7375,
9.1875, 58.8, 10.4125, 49, 22.05, 22.05, 26.95, 40.425, 11.6375,
13.475, 22.05, 1.715, 15.925, 5.5125, 2.695, 10.4125, 5.5125,
0.98, 10.4125, 18.375, 13.475, 10.4125, 40.425, 15.925, 22.05,
0.245, 2.695, 11.6375, 7.9625, 15.925, 58.8, 5.5125, 6.7375,
13.475, 13.475, 10.4125, 10.4125, 13.475, 9.1875, 10.4125, 33.075,
26.95, 22.05, 22.05, 7.9625, 13.475, 18.375, 3.675, 18.375, 6.7375,
6.7375, 15.925, 22.05, 5.5125, 40.425, 18.375, 13.475, 6.7375,
26.95, 0.245, 33.075, 10.4125, 49, 22.05, 2.205, 6.7375, 33.075,
13.475, 9.1875, 18.375, 22.05, 0.245, 13.475, 11.6375, 33.075,
26.95, 11.6375, 22.05, 9.1875, 26.95, 40.425, 18.375, 9.1875,
40.425, 5.5125, 0.245, 6.7375, 0.245, 18.375, 11.6375, 13.475,
13.475, 6.7375, 2.205, 1.715, 7.9625, 2.205, 1.715, 13.475, 10.4125,
22.05, 10.4125, 5.5125, 15.925, 11.6375, 7.9625, 18.375, 18.375,
26.95, 18.375, 13.475, 5.5125, 22.05, 33.075, 33.075, 11.6375,
22.05, 22.05, 13.475, 3.675, 6.7375, 0.98, 18.375, 13.475, 18.375,
26.95, 5.5125, 10.4125, 18.375, 13.475, 13.475, 0.245, 22.05,
11.6375, 0.98, 22.05, 33.075, 6.7375, 9.1875, 11.6375, 6.7375,
0.98, 22.05, 3.185, 22.05, 4.41, 5.5125, 22.05, 15.925, 4.41,
2.695, 40.425, 26.95, 6.7375, 22.05, 11.6375, 6.7375, 33.075,
6.7375, 2.695, 0.245, 33.075, 7.9625, 2.205, 22.05, 18.375, 9.1875,
15.925, 13.475, 3.185, 18.375, 10.4125, 15.925, 13.475, 22.05,
1.715, 3.185, 10.4125, 15.925, 15.925, 26.95, 3.675, 5.5125,
26.95, 13.475, 13.475, 6.7375, 6.7375, 22.05, 22.05, 5.5125,
3.185, 5.5125, 3.675, 22.05, 18.375, 18.375, 6.7375, 7.9625,
18.375, 15.925, 4.41, 5.5125, 11.6375, 18.375, 6.7375, 33.075,
22.05, 11.6375, 4.41, 22.05, 13.475, 6.7375, 26.95, 26.95, 26.95,
22.05, 22.05, 22.05, 6.7375, 33.075, 33.075, 22.05, 15.925, 3.185,
26.95, 1.715, 18.375, 10.4125, 4.41, 11.6375, 22.05, 33.075,
22.05, 15.925, 13.475, 0.245, 5.5125, 22.05, 26.95, 10.4125,
10.4125, 2.695, 40.425, 18.375, 40.425, 18.375, 26.95, 2.205,
0.245, 18.375, 11.6375, 13.475, 10.4125, 13.475, 9.1875, 26.95,
5.5125, 22.05, 5.5125, 13.475, 0.98, 22.05, 13.475, 15.925, 15.925,
22.05, 1.715, 3.185, 11.6375, 1.715, 18.375, 5.5125, 13.475,
22.05, 40.425, 15.925, 40.425, 0.245, 3.675, 18.375, 5.5125,
11.6375, 33.075, 18.375, 6.7375, 13.475, 22.05, 40.425, 7.9625,
1.715, 22.05, 40.425, 10.4125, 13.475, 0.245, 13.475, 5.5125,
2.695, 18.375, 9.1875, 10.4125, 26.95, 11.6375, 1.715, 3.675,
10.4125, 26.95, 40.425, 0.245, 3.185, 15.925, 6.7375, 13.475,
6.7375, 13.475, 3.675, 5.5125, 6.7375, 11.6375, 0.245, 3.675,
11.6375, 10.4125, 26.95, 11.6375, 11.6375, 10.4125, 0.245, 3.675,
13.475, 22.05, 15.925, 33.075, 5.5125, 9.1875, 11.6375, 22.05,
18.375, 13.475, 22.05, 5.5125, 1.715, 15.925, 10.4125, 0.98,
6.7375, 18.375, 18.375, 0.245, 18.375, 0.98, 15.925, 13.475,
13.475, 22.05, 9.1875, 6.7375, 6.7375, 15.925, 10.4125, 33.075,
0.98, 22.05, 13.475, 3.185, 10.4125, 11.6375, 13.475, 22.05,
7.9625, 15.925, 10.4125, 6.7375, 22.05, 4.41, 40.425, 1.715,
40.425, 2.695, 6.7375, 13.475, 5.5125, 3.185, 15.925, 40.425,
33.075, 6.7375, 13.475, 9.1875, 15.925, 9.1875, 1.715, 7.9625,
13.475, 22.05, 22.05, 5.5125, 15.925, 4.41, 7.9625, 15.925, 6.7375,
0.98, 9.1875, 5.5125, 33.075, 22.05, 18.375, 15.925, 15.925,
9.1875, 13.475, 11.6375, 15.925, 6.7375, 9.1875, 1.715, 3.675,
4.41, 22.05, 22.05, 0.98, 1.715, 15.925, 5.5125, 3.185, 6.7375,
9.1875, 0.98, 9.1875, 10.4125, 26.95, 0.98, 5.5125, 2.695, 5.5125,
15.925, 22.05, 9.1875, 5.5125, 2.205, 6.7375, 18.375, 6.7375,
2.695, 0.98, 7.9625, 4.41, 2.205, 11.6375, 5.5125, 0.98, 15.925,
18.375, 13.475, 3.185, 15.925, 1.715, 0.245, 15.925, 1.715, 1.715,
4.41, 10.4125, 1.715, 26.95, 18.375, 7.9625, 5.5125, 6.7375,
11.6375, 11.6375, 11.6375, 5.5125, 7.9625, 33.075, 6.7375, 9.1875,
18.375, 7.9625, 1.715, 22.05, 4.41, 0.245, 2.205, 9.1875, 2.695,
6.7375, 7.9625, 5.5125, 3.675, 15.925, 1.715, 13.475, 10.4125,
9.1875, 13.475, 7.9625, 10.4125, 22.05, 1.715, 15.925)), row.names = c(NA,
-734L), class = c("tbl_df", "tbl", "data.frame"))
model <- lm(income ~ educ_12, data = dat1)
predictions <- data.frame(
educ_12 = dat1$educ_12,
fit = predict(model, interval = "confidence", level = 0.95),
pred = predict(model, interval = "prediction", level = 0.95)
)
colnames(predictions) <- c("educ_12", "fit", "lwr_conf", "upr_conf", "lwr_pred", "upr_pred")
ggplot(dat1, aes(x = educ_12, y = income)) +
geom_point() + # Plot the original data points
geom_line(data = predictions, aes(x = educ_12, y = fit), color = "blue") + # Fitted regression line
geom_ribbon(data = predictions, aes(x = educ_12, ymin = lwr_conf, ymax = upr_conf), alpha = 0.2, fill = "blue", inherit.aes = FALSE) + # 95% confidence interval
geom_ribbon(data = predictions, aes(x = educ_12, ymin = lwr_pred, ymax = upr_pred), alpha = 0.1, fill = "red", inherit.aes = FALSE) + # 95% prediction interval
labs(
title = "Income vs. Education (12-year equivalent)",
x = "Years of Education (12-year equivalent)",
y = "Income"
) +
theme_minimal()
it gives me this:
What do I do wrong ?