load the data

library(car)
data(peas, package="psych")
str(peas)
## 'data.frame':    700 obs. of  2 variables:
##  $ parent: num  21 21 21 21 21 21 21 21 21 21 ...
##  $ child : num  14.7 14.7 14.7 14.7 14.7 ...

initial plot

plot(child ~ parent, data=peas,
    pch=16, cex.lab=1.25, 
    asp=1, xlim=c(14, 23),
    xlab="Parent seed diameter (0.01 in.)", 
    ylab="Child seed diameter (0.01 in.)")

or, just a semi graphic table

calculate frequencies in the combinations of (parent, child)

library(dplyr)
peas.freq <- peas %>%
    group_by(parent, child) %>%
    summarise( count=n() )
peas.freq
## # A tibble: 52 x 3
## # Groups:   parent [?]
##    parent child count
##     <dbl> <dbl> <int>
##  1    15.  13.8    46
##  2    15.  14.8    14
##  3    15.  15.8     9
##  4    15.  16.8    11
##  5    15.  17.8    14
##  6    15.  18.8     4
##  7    15.  19.8     2
##  8    16.  14.3    34
##  9    16.  15.3    15
## 10    16.  16.3    18
## # ... with 42 more rows

make the plot, adding text labels for the counts

plot(child ~ parent, data=peas,
    cex.lab=1.25, 
    asp=1, xlim=c(14, 23),
    xlab="Parent seed diameter (0.01 in.)", ylab="Child seed diameter (0.01 in.)")
with( peas.freq, 
    text( parent, child, count, col="red", cex=log(count))
    )

try a sunflower plot

sunflowerplot(child ~ parent, data=peas,
    pch=16, cex.lab=1.25, 
    asp=1, 
    xlim=c(14, 23),
    xlab="Parent seed diameter (0.01 in.)", 
    ylab="Child seed diameter (0.01 in.)")

fit the linear regression

peas.mod <- lm(child ~ parent, data=peas)

op <- par(mar=c(5, 4, 1, 1)+.1)
plot(jitter(child) ~ jitter(parent), data=peas,
    pch=16, cex.lab=1.25, 
    asp=1, 
    xlim=c(14, 23),
    xlab="Parent seed diameter (0.01 in.)", 
    ylab="Child seed diameter (0.01 in.)")

# plot the regression line
# NB:  could just use `abline(peas.mod, col="blue", lwd=3)` here.
# However, this draws the line across the whole range of X
# Alternatively,  do this using `predict()``

peas.mod <- lm(child ~ parent, data=peas)
pred <- 14:22
lines(pred, predict(peas.mod, data.frame(parent=pred)), col="blue", lwd=3)

# line of unit slope
mx <- mean(peas$parent)
my <- mean(peas$child)
xp <- 14:22
yp <- my + (xp - mx)

# plot and label the lines
lines(xp, yp, col="darkgray", lwd=3, lty="longdash")
text(23.2, yp[9], "child ~ 1 * parent", cex=1.4, col="darkgray")
text(23, y=18.3, "child ~0.34 * parent", cex=1.4, col="blue")

# reference lines of overall means
abline(v=mx, h=my, col="grey")

# add line of means
means <- aggregate(child ~ parent, data=peas, FUN=mean)
lines (child ~ parent, data=means, type="b", pch="+", cex=2, lwd=7, col="darkgreen")
#text(22, 17.5, "Means", col="darkgreen", cex=1.4)
text(15, 15.3, "means", col="darkgreen", cex=1.4, pos=2)


par(op)

car::scatterplot has many other features

scatterplot(child ~ parent, data=peas, 
    ellipse=TRUE, levels=0.68, robust=FALSE, lwd=3, pch=16,
    col=c("blue", "red", "black"  ),
    jitter=list(x=1, y=1), 
    boxplots=FALSE, spread=FALSE,
    xlab="Parent diameter", 
    ylab="Child diameter")