In a hierachical model, we have $$p(x_1, \dots, x_N, z_1, \dots, z_N, \beta) = p(\beta) \prod_{i=1}^N p(z_i | \beta) p(x_i | z_i) $$ In such models, we have $x_i \perp x_j | \beta, i \neq j$. However, in general, we do not have $x_i \perp x_j$, while in general we usually assume that data points $x_{1:N}$ are drawn i.i.d. from some model.
Is this a somewhat contradictory?
The resolution of this issue lies in recognition of the difference between conditional independence versus marginal independence. Conditional independence does not generally imply marginal independence, and indeed, in most cases it implies positive marginal correlation. This occurs because outcomes of one of the random variables in the set gives information on the underlying parameters (the conditioning parameters of the underlying distribution) and this in turn gives information on the other random variables in the set. If you would like a more detailed answer than this, I recommend you read O'Neill (2009). This paper sets out the relationship between conditional and marginal independence in classical and Bayesian models, and it discusses some of the implications of conditional independence. It also gives conditions for induced marginal positive correlation for conditionally IID variables.
In the hierarchical model you have given, the elements $x_1,...,x_N$ are conditionally independent given $\beta$. (You are wrong to assert that they are IID --- this holds only conditionally on the parameter $\beta$ and only assuming that $z_1,...,z_N$ are identically distributed conditional on $\beta$.) Each observed value $x_i$ gives information about the underlying value $z_i$ and through it gives information about $\beta$. Conditional on the underlying parameter value $\beta$ the observations are independent, but once you take that conditioning away they are no longer independent because they each give information on $\beta$, and so they indirectly give information on each other. Unless you are dealing with a pathological case, if you derive the joint density for $x_1,...x_N$ (without conditioning on $\mathbf{z}$ or $\beta$) you will find that there is dependence between the values, and most likely positive correlation. This is a standard result in Bayesian models and a standard outcome for conditionally IID values in general.
Typically in hierarchical models, the variable $Z_i$ is taken to be some categorical variable that groups different observations, not a feature that varies between each and every individual. Taking the SAT example, one could reasonably take $Z_i$ to be the mean score of the school that student $i$ attended, since many students will attend the same school. Taking $Z_i$ as the individual skill level kinda defeats the purpose of using a hierarchical model, since this variable does not group your observations at all.
To make my point clearer, let us take the SAT example. We could use a model where students' scores are normally distributed around their school's mean score, that is:
$$X_i|Z_i\sim\mathcal N(Z_i,\sigma^2)$$
Then, to complete the hierarchical model, we say the schools' mean are normally distributed around a grand mean:
$$Z_i\sim\mathcal N(\mu,\alpha)$$
Now, one has to define priors for $\mu$, $\sigma^2$ and $\alpha$ and presto: the model is now fully specified.
The advantages of this approach include:
If you use $Z_i$ to represent individual characteristics, such as skill level, you lose these advantages. Such a model would be perfectly valid and might be technically a hierarchical model, but a rather atypical one. (Unless you are considering each student takes the test more than once, then it would be a usual hierarchical model, since the variable "students" groups the variable "test score")
Now, suppose students $i$ and $j$ go to the same school (that is, $Z_i=Z_j$). Their SAT scores are dependent: if we learn $i$ got a good score, the probability of $Z_i$ being a good school increases and the expected score for $j$ increases. The only way we can make their scores independent is conditioning in their school's mean score, which you denoted as $x_i\perp x_j|\beta$.
This model does not suppose observations are iid - the observations are neither independent now identically distributed, i they have different values of $Z$. There is no contradiction in this, not all statistical models need iid assumptions. Most do, but not all.
You can find more info about hierarchical models and their implementation in PyMC3 in this link Hope that was helpful!
