Why is my method/(R) code for calculating the proportion of values outside x standard deviations incorrect? Imagine I need to calculate the proportion of samples greater than one standard deviation from the mean.
For an ideal, normally distributed population in theory in R:
pnorm(1, mean=0, sd = 1) - pnorm(-1, mean=0, sd = 1)

N.B. I add mean = and sd = for illustration here, as these are the default values in R
By the same logic, for my dataset y:
library(rafalib)
y.sd <- popsd(y)  
y.mean <- mean(y)  

pnorm((y.mean + y.sd), mean = y.mean, sd = y.sd) - pnorm((y.mean - y.sd), mean = y.mean, sd = y.sd)

However, this gives me an incorrect answer, and instead I am supposed to use:
z <- ( y - mean(y) ) / popsd(y)
mean( abs(z) <=1 )

While I understand what the second block of code accomplishes, I don't understand why my code/method does not work?
 A: In your code, you are looking for proportion of samples within 1 standard deviation, not greater.
So this code:
z <- ( y - mean(y) ) / popsd(y)
mean( abs(z) <=1 )

popsd() actually returns the biased sample standard deviation and uses that to estimate the proportion... so it's a bit odd, but regardless, it checks among your samples, how many of have a difference of less than 1 sd from the mean, both estimated from the sample. Straight forward.
This works regardless of whether your sample size is small or large, follows a normal distribution or not.
Your approach goes one step further, you actually assume the data follows a normal distribution and it will always give one answer which is:
pnorm(1) - pnorm(-1)
[1] 0.6826895

So in situations where it does not follow a normal distribution, it will not work. Or for example, in your sample, the size is small or the distribution is slightly biased, you will see the theoretical value differs from what you see in the sample:
set.seed(777)
y = rnorm(50,2,5)
y.mean <- mean(y)  
y.sd = popsd(y)

z <- ( y - y.mean ) / y.sd
mean( abs(z) <=1 )
[1] 0.7

So if the question or the task you have at hand is to simply tell how much of your samples fall within 1 sd of its mean, I guess you need to actually calculate it from the sample.
A: First try a small sample with $n = 50$ observations from
$\mathsf{Norm}(\mu = 100, \sigma=15).$
set.seed(624)
x = rnorm(50, 100, 15)

Find the sample average a and SD s.
a = mean(x);  s = sd(x);  a;  s
[1] 98.67048
[1] 17.71619

Then you want the proportion of the $n = 50$ observations more than
one SD from the mean, that's the number below $\bar X - S$ or above
$\bar X + S.$ In R, you can make a logical vector with TRUEs for
the observations outside the interval $\bar X \pm S$ and FALSEs for
those inside. The mean of a logical vector is the proportion of its
TRUEs. (In combining logical vectors, R uses | for or.)
For clarity, begin by sorting the sample vector x in order from smallest to largest.
x = sort(x)
> x < a - s
 [1]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE
[11] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[21] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[31] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[41] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
x > a + s
 [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[11] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[21] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[31] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[41] FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
(x < a-s) | (x > a + s)
 [1]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE
[11] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[21] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[31] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[41] FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
mean((x < a-s) | (x > a + s))
[1] 0.36

There are
$9+9=18$ TRUEs out of 50 so that's 36%.
stripchart(x, pch="|")
  abline(v = a-s, col="red");  abline(v = a+s)


In any normal population the probability outside of $\mu \pm \sigma$ is
about $0.317 = 31.7\%.$ In a sample as small as $n = 50$ you can expect to
get somewhere near 31.7% of observations outside $\bar X \pm S,$ but not
exactly.
2 * pnorm(-1)
[1] 0.3173105
1 - diff(pnorm(c(-1,1)))
[1] 0.3173105

For a huge sample of size $n = 100\,000,$ you can expect to come
much closer to 31.7%, as follows:
set.seed(2020)
y = rnorm(10^5, 100, 15)
a = mean(y);  s = sd(y);  a;  s
[1] 99.9584
[1] 14.95503
mean((y < a-s)|(y > a+s))
[1] 0.31763

hist(y, prob=T, br=50, col="skyblue2")
curve(dnorm(x,100,15), add=T, lwd=2)  
      # 'x' argument mandatory for 'curve'
abline(v=c(a-s,a+s), col="red")


