Deriving the formula for the Möbius function?

Yes, there is! It is proved in the same way you prove the Möbius inversion formula:

Going directly from the formula you give, it is clear that $$ \mu(1) = 1, \mu(p) = -1 \text{ for $p$ a prime}. $$ In general, for square-free numbers, we can show the desired formula by principle of inclusion-exclusion. For instance, for $n = pq, p, q$ primes, we have $$ \mu(pq) = \sum_{d|pq} \mu(d) - \sum_{d|p} \mu(d) - \sum_{d|q} \mu(d) + \sum_{d|1} \mu(d) = 1 $$ and for three primes, $$ \mu(pqr) = \left (\sum_{d|pqr} - \sum_{d|pq} - \sum_{d|qr} - \sum_{d|pr} + \sum_{d|p} + \sum_{d|q} + \sum_{d|r} - \sum_{d|1} \right )\mu(d) = -1. $$ Now PIE works for the general case as well. If $n = p_1^{e_1} \cdots p_s^{e_s}$, let $m = p_1^{e_1 - 1 } \cdots p_s^{e_s -1}$, i.e. $n$ divided by its square-free part. Then we get that $$\small \mu(n) = \sum_{d|n} \mu(d) - \sum_{p_1|n} \left( \sum_{d|n/p_1} \mu(d) \right ) + \sum_{p_1, p_2 |n, p_1 \neq p_2} \left ( \sum_{d|n/(p_1 p_2)} \mu(d) \right ) -\cdots + (-1)^s \sum_{d|m} \mu(d) = 0 $$ if $m \neq 1$, as desired.

Even though the question demands forgetting about posets, I'll summarize the general approach of Möbius inversion here, since it is no less elementary than a "direct" approach, does not rely on deep insight, and is in fact quite painless.

Suppose $(P,\preceq))$ are a set $P$ and a partial ordering $\preceq$ on $P$ (for instance $P=\mathbb N_{\geq1}$, and $x\preceq y$ means $x\mid y$), and suppose that for each $y\in P$ there are only finitely many $x\in P$ with $x\preceq y$, so that our summations below will make sense. We want a function $\mu(x,y)$ on $P^2$ with the property $$ \sum_{y\preceq z}\mu(x,y)=\delta_{x,z}. $$ for all $x,z\in P$. This equation defines $\mu(x,z)$ once all other $\mu(x,y)$ that occur are known, so the values $\mu(x,y)$ are uniquely defined by this equation, as can be shown for fixed $x$ by strong induction on $y$ along the partial ordering $\preceq$ (or along any total ordering extending it, if one prefers). This shows in particular that $\mu(x,y)=0$ whenever $x\not\preceq y$, and $\mu(x,x)=1$ for all $x$.

In our example we will get $\mu(x,y)=\mu(y/x)$ with $\mu$ now denoting the classical (arithmetic) Möbius function. Indeed for $x=1$ the equation becomes the usual $\sum_{y\mid z}\mu(1,y)=\delta_{1,z}$, and for $x>1$ one can easily check that $\mu(x,y)=\mu(1,y/x)$ for all multiples $y$ of $x$. It is somewhat deranging that the general setting requires a function of two arguments rather than one, but that setting provides no means to express anything like the quotient operation used to reduce to a single argument; nevertheless in concrete situations there is often a similar simple function of $x$ and $y$ in terms of which $\mu(x,y)$ is defined. Note that the equation is asking $\mu(x,y)$ to be the coefficients of the inverse of the matrix $M$, with rows and columns indexed by elements of $P$, given by $M_{x,y}=1$ if $x\preceq y$ and $M_{x,y}=0$ otherwise. Since a left inverse is also right inverse, this point of view allows us to recast the above equation as $$ \sum_{x\preceq y\preceq z}\mu(y,z)=\delta_{x,z}. $$ But enough general theory.

If $\preceq$ is actually a total ordering on $P$, which means (given the finiteness condition) that $P$ is isomorphic to an initial segment of $\mathbb N$ (possibily all of $\mathbb N$), then it is easy to see that one has $\mu(x,y)=\epsilon(y-x)$ where $\epsilon:\mathbb Z\to\{-1,0,1\}$ satisfies $\epsilon(i)=(-1)^i$ if $i\in\{0,1\}$ and $\epsilon(i)=0$ otherwise. This almost trivial computation will allow deducing many important Möbius functions, including the classical Möbius function, using the following

Theorem. If $(P_i,\preceq_i)$ are appropriate posets for $i\in I$, and $(P,\preceq)$ is the product poset (restricted if $I$ is infinite), then the Möbius function for $P$ is the product of the Möbius functions for the posets $P_i$: one has $\mu((x_i)_{i\in I},(y_i)_{i\in I})=\prod_{i\in I}\mu_i(x_i,y_i)$.

Here a product poset is the Cartesian product $\prod_{i\in I}P_i$ with a component-wise partial ordering: $(x_i)_{i\in I}\leq(y_i)_{i\in I}$ if and only if $x_i\leq y_i$ in $P_i$ for all $i\in I$. Also "restricted" means restriction to the subset of the Cartesian product of those $(x_i)_{i\in I}$ such that $x_i$ is a minimal element of $P_i$ for all but finitiely many indices $i$; this restriction is necessary to ensure that there are only finitely many elements below any given element of $P$, so that its Möbius function is well defined. The example that will interest us is with $I$ the set $\mathrm{Pr}$ of prime numbers, and every $(P_i,\preceq_i)$ a copy of $(\mathbb N,\leq)$; then elements of $P$ are collections $(m_p)_{p\in\mathrm{Pr}}$ with $m_p\in\mathbb N$ and $m_p=0$ for all but finitiely many $p$. Such an element can be associated to the product $\prod_{p\in\mathrm{Pr}}p^{m_p}$, which is a bijection of $P$ to the set $\mathbb N_{\geq1}$, under which bijection the partial order $\preceq$ on $P$ corresponds to the divisibility relation on $\mathbb N_{\geq1}$. Thus we see that the theorem describes in particular the classical Möbius function. Indeed it gives $\mu(x,y)=\prod_{p\in\mathrm{Pr}}\epsilon(m_p(y)-m_p(x))$ where $m_p(x)$ denotes the multiplicity of the prime $p$ in the factorisation of $x$; the product equals $\prod_{p\in\mathrm{Pr}}\epsilon(m_p(y/x))$ whenever $x\mid y$, which is the usual definition of $\mu(y/x)$.

Proof. This is simply a formal verification that the proposed Möbius function for $P$ satifies the required equation. Let $x=(x_i)_{i\in I}$ and $z=(z_i)_{i\in I}$ be elements of $P$, and let $J=\{i_1,\ldots,i_n\}$ be a finite set containing all indices for which $z_i$ is not minimal. Then $$ \begin{align} \sum_{y\preceq z}\mu(x,y) &=\sum_{y_{i_1}\preceq z_{i_1}}\cdots\sum_{y_{i_n}\preceq z_{i_n}} \mu_{i_1}(x_{i_1},y_{i_1})\cdots\mu_{i_n}(x_{i_n},y_{i_n}) \\ & =\prod_{j\in J}\,\sum_{y_j\preceq z_j}\mu_j(x_j,y_j) =\prod_{j\in J}\delta_{x_j,z_j} =\delta_{x,z}.\quad\mbox{QED}\\ \end{align} $$

Setting $n=1$ we deduce $\mu(1)=1$ we could continue with n=p prime to get $\mu(1)+\mu(p)=0 \implies \mu(p)=-1$, try the product of two primes to get $\mu(p\cdot q)=1$ and so on...

We could too define the arithmetical function u : u(k)=1 for all k and rewrite your definition as $\mu * u = I$ (of course this is just a shortcut for $\sum_{d|n}\mu(d) u(\frac{n}{d})=I(n)$ with $I$ the function of n defined at the right of your equation).
This and $\mu(1)=1$ is enough to prove the Möbius inversion formula (details in Apostol's excellent introduction to A.N.T. page 30 to 32 : the proofs there only require your definition). You could prove that u and I are multiplicative ($f(m n)=f(m)f(n)$ if $\gcd(m,n)=1$) so that $\mu$ will also be multiplicative (using/proving Theorem 2.16 page 36). Of course this is only the draft of an action plan and I'll let you or others concretize it! :-)