We present an $L_q(L_{p})$-theory for the equation $$ \partial_{t}^αu=φ(Δ) u +f, \quad t>0,\, x\in \mathbb{R}^d \quad\, ;\, u(0,\cdot)=u_0. $$ Here $p,q>1$, $α\in (0,1)$, $\partial_{t}^α$ is the Caputo fractional derivative of order $α$, and $φ$ is a Bernstein function satisfying the following: $\exists δ_0\in (0,1]$ and $c>0$ such that \begin{equation} \label{eqn 8.17.1} c \left(\frac{R}{r}\right)^{δ_0}\leq \frac{φ(R)}{φ(r)}, \qquad 0<r<R<\infty. \end{equation} We prove uniqueness and existence results in Sobolev spaces, and obtain maximal regularity results of the solution. In particular, we prove \begin{align*} \| |\partial^α_t u|+|u|+|φ(Δ)u|\|_{L_q([0,T];L_p)}\leq N(\|f\|_{L_q([0,T];L_p)}+ \|u_0\|_{B_{p,q}^{φ,2-2/ αq}}), \end{align*} where $B_{p,q}^{φ,2-2/αq}$ is a modified Besov space on $\mathbb{R}^d$ related to $φ$. Our approach is based on BMO estimate for $p=q$ and vector-valued Calderón-Zygmund theorem for $p\neq q$. The Littlewood-Paley theory is also used to treat the non-zero initial data problem. Our proofs rely on the derivative estimates of the fundamental solution, which are obtained in this article based on the probability theory.