Sums of the form $\sum\limits_{a_1=1}^n\sum\limits_{a_2=1}^{a_1}\cdots \sum\limits_{a_r=1}^{a_{r-1}}{[\gcd(a_1,a_2,\dots ,a_r)=1]}$ Background:
Consider sums of the form
$$\sum_{a_1=1}^n\sum_{a_2=1}^{a_1}\cdots \sum_{a_r=1}^{a_{r-1}}{[\gcd(a_1,a_2,\dots ,a_r)=1]},$$
with $[\gcd(a_1,a_2,\dots ,a_r)=1]$ being equal to $1$ if the condition holds true and $0$ otherwise. After a bit of experimentation it appears that
$$\sum_{a_1=1}^n\sum_{a_2=1}^{a_1}\cdots \sum_{a_r=1}^{a_{r-1}}{[\gcd(a_1,a_2,\dots ,a_r)=1]}=\sum_{m\leq n}\sum_{d\mid m}\mu\left(\frac{m}{d}\right){d+r-2\choose r-1}. \tag 1$$
This method is really fast and I would love to know how to simplify this further. On the other hand we have the well known result
$$\sum_{i=1}^n\sum_{j=1}^{i}{[\gcd(i,j)=1]}=\sum_{k=1}^n\varphi(k)=\frac{3n^2}{\pi^2}+O\left(n\log n\right) \tag 2$$
which we could use to get something like
$$\sum_{a_1=1}^n\sum_{a_2=1}^{a_1}\cdots \sum_{a_r=1}^{a_{r-1}}{[\gcd(a_1,a_2,\dots ,a_r)=1]}= \frac{n^r}{r!\zeta(r)}+O(?). \tag 3$$

Questions:
1) How can we further simplify $(1)$?

2) What's the error term in $(3)$? I'm not sure what I'm doing wrong but I can get it lower than $n^{r-1}\log n^{r-1}$, which seems ridiculously high.

3) Is there anything to be found in the literature about sums of the above form?

Addendum (2019-04-24):
From A002088 we also have
$$\sum_{k=1}^n\varphi(k)=\frac{1}{2}\sum_{k=1}^n\mu(k)\left\lfloor\frac{n}{k}\right\rfloor\left\lfloor\frac{n}{k}+1\right\rfloor \tag 4$$
which we can use to get
$$\sum_{a_1=1}^n\sum_{a_2=1}^{a_1}\cdots \sum_{a_r=1}^{a_{r-1}}{[\gcd(a_1,a_2,\dots ,a_r)=1]}=\frac{1}{r!}\sum_{k=1}^n\mu(k)\prod_{j=0}^{r-1}\left\lfloor\frac{n}{j}+j\right\rfloor. \tag 5$$
Is this a better approach? Still, I don't know how to simplify the RHS sum-product.
 A: As you note in the addendum you added, the sum can be re-arranged as:
$$\sum_{a_{1}=1}^{n}\sum_{a_{2}=1}^{a_{1}}\cdots\sum_{a_{r}=1}^{a_{r-1}}\sum_{d|a_{1},\dots,a_{r}}\mu(d)=\sum_{\substack{a_{r}\leq\cdots\leq a_{1}\\
a_{1}\leq n
}
}^{n}\sum_{d|a_{1},\dots,a_{r}}\mu(d)
$$
$$=\sum_{d\leq n}\mu(d)\sum_{\substack{a_{r}\leq\cdots\leq a_{1}\\
a_{1}\leq\frac{n}{d}
}
}^{n}1$$
$$=\sum_{d\leq n}\mu(d)\binom{\lfloor\frac{n}{d}\rfloor+r-1}{r}.$$
By expanding the polynomial/binomial coefficient, and then expanding the floor function as $\lfloor x\rfloor=x-\{x\}$, one can show that the original sum is equal to:
$$\frac{n^r}{r!\zeta(r)}+O(n^{r-1})$$ for $r>2$. Note that when $r=2$, the error term is larger, and is provably $$\Omega_{\pm}\left(n\sqrt{\log \log n}\right).$$ This is a 1987 result of Montgomery, see this answer for more details: Asymptotic formula for $\sum_{n\leq x}\mu(n)[x/n]^2$ and the Totient summatory function $\sum_{n\leq x} \phi(n)$ 
For $r>2$, I do not believe that the error term can be improved from $O(n^{r-1})$. The issue is that you have to deal with a sum of the form 
$$n^{r-1}\sum_{d=1}^\infty \frac{\mu(d)}{d^{r-1}} \left\{\frac{n}{d}\right\}.$$ The fractional part term has an average of $$\frac{1}{n}\sum_{d\leq n} \left\{\frac{n}{d}\right\} \sim 1-\gamma$$ where $\gamma$ is the Euler-Mascheroni constant, and one could hope that $$\sum_{d=1}^\infty\frac{\mu(d)}{d^{r-1}} \left\{\frac{n}{d}\right\}=^?\frac{1-\gamma}{\zeta(r-1)}+o(1).$$ If this could be proven, then there would be a term of size $n^{r-1}$ appearing, and by including this we could improve the main term. However, this equality is not true. One way to see that is by calculating the mean, and then the standard deviation of $$R(x)=\sum_{d\leq x}\frac{\mu(d)}{d^{r-1}} \left\{\frac{x}{d}\right\}.$$ The mean is $1-\gamma$, and the standard deviation is a constant $>0$. This implies that for the original sum, there is some kind of oscillating term of size $n^{r-1}$ that needs to be understood, determined, and included as part of the main term, if a better error term is desired.
