The increase of operating frequencies requires-D electromagnetic (EM) methods, such as the partial element equivalent circuit (PEEC) method, for the analysis and design of high-speed circuits. Very large systems of equations are often produced by ${\bf 3}$-D EM methods and model order reduction (MOR) techniques are used to reduce such a high complexity. When signal waveform rise times decrease and the corresponding frequency content increases, or the geometric dimensions become electrically large, time delays must be included in the modeling. A PEEC formulation, which include delay elements called PEEC method, becomes necessary and leads to systems of neutral delayed differential equations (NDDE). The reduction of large NDDE is still a very challenging research topic, especially for electrically large structures, where delays among coupled elements cannot be neglected or easily approximated by rational basis functions. We propose a novel model order technique for PEEC models that is able to accurately reduce electrically large systems with large delays. It is based on an adaptive multipoint expansion and MOR of equivalent first-order systems. The neutral delayed differential formulation is preserved in the reduced model. Pertinent numerical examples based on PEEC models validate the proposed MOR approach.