This is true and this is not about numbers $j^j$ at all. Consider a more genral matrix $(c_j-i^j)_{1\leqslant i,j\leqslant n}$. Denote $f_j(x)=c_j-x^j$, and find the numbers $\beta$, $\alpha_1,\ldots,\alpha_{n-1}$ such that $f_n(x)+\sum_{j=1}^{n-1}\alpha_j f_j(x)=\beta$ for all $x\in \{1,2,\ldots,n\}$. Then the matrix $(f_j(i))_{i,j}$ has the same determinant as if we replace the column with $f_n$ by a column of $\beta$'s, that is seen from operations with columns. Now by operations with columns we may replace each $f_j$ to $-x^j$, so we get almost a Vandermonde matrix (up to $n-1$ changes of signs in the columns and the cyclic shift, and $\beta$'s instead of 1's), so its determinant equals to the determinant of Vandermonde matrix for the numbers $1,\ldots,n$ times $(-1)^{n-1}(-1)^{n-1}\beta=\beta$. Thus, it appears to find $\beta$. For this, write $$ f_n+\sum_{j=1}^{n-1}\alpha_jf_j(x)=\beta-(x-1)(x-2)\ldots(x-n)= \beta-\sum_{k=0}^n s(n+1,k+1)x^k. $$ Equalising the coefficients of $x^k$ we get $\alpha_k=s(n+1,k+1)$ for all $k=1,\ldots,n$ (where $\alpha_n=1$ by agreement), thus $\beta-s(n+1,1)=\sum_{k=1}^n c_k s(n+1,k+1)$, $\beta=\sum_{k=0}^n c_k s(n+1,k+1)$, where $c_0=1$ by agreement.
