Podem haver outras abordagens, mas considere essa de início.

da <- with(subset(sub_data, tratamento=="A"),
           list(y=imoveis, n=moveis+imoveis))

library(bbmle)

fitA <- mle2(y~dbinom(prob=p, size=n),
             start=list(p=0.01),
             data=da,
             optimizer="nlminb",
             lower=c(p=10e-5))
confint(fitA)

Eu tive que fazer uma restrição de borda à esquerda para que aconteça a convergência (lower=10e-5). O valor estimado é justamente o valor que usei na borda. Tentei valores menore mas ocorre falha (under/overflow numérico), então mantive esse. Isso não é um problema, porque sabendo que todos os y são zero, a estimativa amostral de p = y/n é zero (fora do espaço, claro). De qualquer maneira, essa rotina é para calcular o IC e o confint retorna o limite superior. Existe algumas preocupações sobre esses IC calculados para estimativas de borda mas use em complemento ao fato da estimativa amostral ter sido zero e que esse IC de verossimilhança é o melhor que se pode fazer para reportar alguma incerteza sobre a estimativa.

​À
​ disposição.
Walmes.​