"""Fase 1b, parte B: agrupamentos espaciais (estatística de varredura de Kulldorff, Poisson). Plano registrado em PLANO-fase1b.md (antes de rodar). Resumo: - janelas circulares centradas no centroide de cada município, com os vizinhos mais próximos, até 100 km de raio e até 1 milhão de habitantes (Censo 2022); - óbitos esperados em relação ao próprio estado (municipios_rmp.csv), de modo que os agrupamentos são excessos DENTRO do estado; só janelas com observado > esperado; - valor de p por Monte Carlo (999 réplicas, multinomial condicionada ao total de óbitos); - agrupamento mais provável e secundários sem sobreposição com p < 0,05 (os secundários são comparados com a mesma distribuição do máximo, o que é conservador, como no SaTScan); - esôfago é a análise principal; os outros 12 grupos são exploratórios (sem correção entre eles); - os agrupamentos de esôfago são testados em 2015–2019 (RMP agrupada, replica se > 1 e p < 0,05). Saídas: resultados/varredura_agrupamentos.csv, varredura_municipios.csv, varredura_validacao_esofago.csv Uso: python varredura.py (depois de analisar.py e validar.py) """ import numpy as np import pandas as pd from scipy import stats from analisar import carregar, centroides, comparar from config import ANOS_VALIDACAO, GRUPOS, PROCESSADOS, RESULTADOS from validar import ic_exato, populacao_2017 RAIO_KM = 100 POP_MAX = 1_000_000 REPLICAS = 999 ALFA = 0.05 SEMENTE = 20261005 def janelas(lat, lon, pop): """Matriz (centros × K) com os índices dos vizinhos em ordem de distância (-1 = vazio).""" la, lo = np.radians(lat), np.radians(lon) listas, raios = [], [] for i in range(len(lat)): h = np.sin((la - la[i]) / 2) ** 2 + np.cos(la[i]) * np.cos(la) * np.sin((lo - lo[i]) / 2) ** 2 d = 6371 * 2 * np.arcsin(np.sqrt(h)) ordem = np.argsort(d, kind="stable") ordem = ordem[d[ordem] <= RAIO_KM] ordem = ordem[np.cumsum(pop[ordem]) <= POP_MAX] listas.append(ordem) raios.append(d[ordem]) k = max(len(o) for o in listas) viz = np.full((len(lat), k), -1, dtype=np.int64) dist = np.full((len(lat), k), np.nan) for i, (o, r) in enumerate(zip(listas, raios)): viz[i, :len(o)] = o dist[i, :len(o)] = r return viz, dist def llr(c, e, total): """Log-verossimilhança de Kulldorff (Poisson); 0 se não há excesso na janela.""" with np.errstate(divide="ignore", invalid="ignore"): dentro = np.where(c > 0, c * np.log(c / e), 0.0) fora = np.where(total - c > 0, (total - c) * np.log((total - c) / (total - e)), 0.0) v = dentro + fora return np.where((c > e) & (e > 0), v, 0.0) def acumular(valores, viz): v = np.where(viz >= 0, valores[np.clip(viz, 0, None)], 0.0) return np.cumsum(v, axis=1) def varrer(obs, esp, viz, rng): total = obs.sum() esp = esp * total / esp.sum() # esperado sob a hipótese nula, somando o mesmo total valido = viz >= 0 ce, ee = acumular(obs, viz), acumular(esp, viz) l = np.where(valido, llr(ce, ee, total), 0.0) probs = esp / esp.sum() maximos = np.empty(REPLICAS) for r in range(REPLICAS): sim = rng.multinomial(int(total), probs).astype(float) cs = acumular(sim, viz) maximos[r] = np.where(valido, llr(cs, ee, total), 0.0).max() return l, ce, ee, maximos, total def agrupamentos(l, ce, ee, viz, dist, maximos, total): """Mais provável e secundários sem sobreposição, em ordem de LLR.""" ordem = np.argsort(l, axis=None)[::-1] usados = set() saida = [] for idx in ordem: i, j = np.unravel_index(idx, l.shape) v = l[i, j] if v <= 0: break p = (1 + (maximos >= v).sum()) / (REPLICAS + 1) if p >= ALFA: break membros = viz[i, :j + 1] if usados.intersection(membros.tolist()): continue usados.update(membros.tolist()) c, e = ce[i, j], ee[i, j] saida.append({"centro": i, "membros": membros, "raio_km": dist[i, j], "observados": c, "esperados": e, "rmp": c / e, "risco_relativo": (c / e) / ((total - c) / (total - e)), "llr": v, "p": p}) return saida def esperados_validacao(grupo): """Observados e esperados (em relação ao estado) em 2015–2019, por município.""" pop22, _, _ = carregar() pop, _ = populacao_2017(pop22) ob = pd.read_csv(PROCESSADOS / f"obitos_agregados_{ANOS_VALIDACAO[0]}_{ANOS_VALIDACAO[-1]}.csv") ob["uf"] = ob["cod6"] // 10000 ob = ob[ob["cod6"].isin(set(pop["cod6"])) & (ob["grupo"] == grupo)] idx = pd.Index(sorted(pop["cod6"].unique()), name="cod6") obs = ob.groupby("cod6")["obitos"].sum().reindex(idx).fillna(0) uf = comparar(obs, pop, ob, "uf") return pd.DataFrame({"obs_val": obs, "esp_val": uf["esperados"]}) def main(): res = pd.read_csv(RESULTADOS / "municipios_rmp.csv") cent = centroides() muns = res.drop_duplicates("cod_ibge").sort_values("cod_ibge").reset_index(drop=True) muns = muns[muns["cod_ibge"].isin(cent)].reset_index(drop=True) lat = np.array([cent[c][0] for c in muns["cod_ibge"]]) lon = np.array([cent[c][1] for c in muns["cod_ibge"]]) viz, dist = janelas(lat, lon, muns["populacao"].to_numpy(float)) print(f"{len(muns)} municípios; até {viz.shape[1]} municípios por janela; " f"{int((viz >= 0).sum()):,} janelas") rng = np.random.default_rng(SEMENTE) linhas, membros_linhas = [], [] for grupo in GRUPOS: d = res[res["grupo"] == grupo].set_index("cod_ibge").reindex(muns["cod_ibge"]) obs = d["observados"].fillna(0).to_numpy(float) esp = d["esperados_uf"].fillna(0).to_numpy(float) l, ce, ee, maximos, total = varrer(obs, esp, viz, rng) achados = agrupamentos(l, ce, ee, viz, dist, maximos, total) print(f"{GRUPOS[grupo]['nome']:42s} agrupamentos com p < 0,05: {len(achados)}") for n, a in enumerate(achados, 1): m = muns.loc[a["membros"]] dm = d.loc[m["cod_ibge"]] centro = muns.loc[a["centro"]] sinais = dm[dm["sinal"].fillna(False).astype(bool)] linhas.append({ "grupo": grupo, "nome_grupo": GRUPOS[grupo]["nome"], "ordem": n, "principal_ou_exploratorio": "principal" if grupo == "esofago" else "exploratório", "centro": centro["municipio"], "uf_centro": centro["uf_sigla"], "cod_ibge_centro": centro["cod_ibge"], "raio_km": round(a["raio_km"], 1), "municipios": len(m), "ufs": ",".join(sorted(m["uf_sigla"].unique())), "populacao": int(m["populacao"].sum()), "observados": int(a["observados"]), "esperados": a["esperados"], "rmp": a["rmp"], "risco_relativo": a["risco_relativo"], "llr": a["llr"], "p_monte_carlo": a["p"], "inclui_capital": bool(dm["capital"].any()), "sinais_dentro": len(sinais), "sinais_prioritarios_dentro": int((sinais["contexto"] == "prioritário").sum()), }) for cod in m["cod_ibge"]: membros_linhas.append({"grupo": grupo, "ordem": n, "cod_ibge": cod, "municipio": d.loc[cod, "municipio"], "uf_sigla": d.loc[cod, "uf_sigla"], "observados": int(d.loc[cod, "observados"]), "esperados_uf": d.loc[cod, "esperados_uf"], "sinal": bool(d.loc[cod, "sinal"]), "contexto": d.loc[cod, "contexto"]}) ag = pd.DataFrame(linhas) mb = pd.DataFrame(membros_linhas) ag.to_csv(RESULTADOS / "varredura_agrupamentos.csv", index=False, float_format="%.4g") mb.to_csv(RESULTADOS / "varredura_municipios.csv", index=False, float_format="%.4g") # Validação dos agrupamentos de esôfago em 2015–2019 val = esperados_validacao("esofago") linhas = [] for _, a in ag[ag["grupo"] == "esofago"].iterrows(): cods = mb[(mb["grupo"] == "esofago") & (mb["ordem"] == a["ordem"])]["cod_ibge"] // 10 o, e = val.loc[cods, "obs_val"].sum(), val.loc[cods, "esp_val"].sum() inf, sup = ic_exato(o, e) p = stats.poisson.sf(o - 1, e) linhas.append({"ordem": a["ordem"], "centro": a["centro"], "uf_centro": a["uf_centro"], "municipios": a["municipios"], "rmp_2020_24": a["rmp"], "obs_2015_19": int(o), "esp_2015_19": e, "rmp_2015_19": o / e, "ic95_inf": inf, "ic95_sup": sup, "p": p, "replica": bool(o / e > 1 and p < ALFA)}) vv = pd.DataFrame(linhas) vv.to_csv(RESULTADOS / "varredura_validacao_esofago.csv", index=False, float_format="%.4g") pd.set_option("display.width", 250, "display.max_columns", 30) print() print(ag[["nome_grupo", "ordem", "centro", "uf_centro", "raio_km", "municipios", "observados", "esperados", "rmp", "llr", "p_monte_carlo", "inclui_capital", "sinais_dentro", "sinais_prioritarios_dentro"]] .round(3).to_string(index=False)) print() print(vv.round(3).to_string(index=False)) return ag, mb, vv if __name__ == "__main__": main()