Why doesn't Var(x) equal mean(x^2) - mean(x)^2 in Rstudio? I'm probably being very dumb, but I'm trying to figure out the variance by hand using R, and I am using a dummy dataset to figure out the formula. My dataset is
x = (1,1,2,3,3,3,3,3)
Using var(x) I get 0.8392857, however when using mean(x^2)-mean(x)^2 I get  0.734375. In my mind, these values should be the same, could someone explain to me why they are not?

The var function uses an $n-1$ denominator to calculate the variance.
For your decomposition, the variance is defined in terms of the expected value of the squared deviations from the mean, which translates into an $n$ denominator.
Thus, you are using the wrong formula to calculate variance when you call the var function. The formula you give applies to the population values of mean, variance, and the uncentered second moment, not to estimates of those values. To calculate population moments from a finite set of numbers like you have, you apply the expected value formula to the numbers. That is, you add up some values and then divide by the number of values.
If you set x = rep(c(1,1,2,3,3,3,3,3), 1000) to repeat your $x$-values a thousand times (thus increasing the $n$ by a factor of a thousand), the difference between $n$ and $n-1$ will be very small, and your numbers will agree out to many decimal places.
If you write your own function to calculate the $n$-denominator variance, then your values will agree.
x <- c(1,1,2,3,3,3,3,3)
var_n <- function(x){
    return(
        mean((x - mean(x))^2)
    )
}
(mean(x^2) - mean(x)^2) - var_n(x) # I get 0
(mean(x^2) - mean(x)^2) - var(x)   # I get -0.1049107, same as your 0.734375 - 0.8392857

I do not know an R function offhand that calculates an $n$-denominator variance. For Python users, numpy.var uses an $n$-denominator variance by default and has to have ddof = 1 passed as an argument to calculate variance with a denominator of $n-1$.
