Can I work around a small number of samples in a logistic regression? [duplicate]

I am deeply frustrated, I have a dataframe with a small number of samples (N=13) and am trying to do a logistic regression with binomial variables. I know, that especially given that there is an even split in my outcome I can't really do this, but there is no way of me increasing my N and this is all I have to work with. I know it has been previously discussed, but I am looking to see if I show the following it might make some sense. If I show that a) there is deviance between the null and the model (using chi-squared) b) the summary of the model is suggestive of separation (which it is) and whether or not I can then do c) and still use the predict() function to determine whether or not each N was being correctly assigned, which it appears is happening. What are your thoughts on this and does it maybe make sense to explain all of this knowing that it is difficult to work with such small samples?

Dataframe:

coxy <- structure(list(first_cyto = c(0, 1, 1, 0, 0, 1, 1, 1, 1, 1, 0,
1, 1), second_cyto = c(1, 1, 0, 0, 0, 1, 0, 1, 1, 0, 0, 1, 1),
third_cyto = c(0, 0, 1, 0, 1, 1, 1, 1, 0, 0, 0, 1, 0), fourth_cyto = c(0,
0, 0, 1, 0, 1, 1, 0, 1, 0, 0, 0, 0), fifth_cyto = c(0, 0,
0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0), sixth_cyto = c(0, 0, 0,
0, 0, 0, 1, 1, 0, 0, 0, 0, 0), seventh_cyto = c(0, 0, 0,
0, 0, 1, 1, 1, 0, 1, 0, 1, 1), eighth_cyto = c(0, 0, 0, 0,
0, 0, 1, 1, 0, 0, 0, 1, 0), ninth_cyto = c(0, 1, 1, 0, 0,
0, 1, 1, 0, 0, 0, 1, 1), tenth_cyto = c(0, 0, 1, 0, 1, 0,
1, 1, 0, 0, 0, 1, 0), eleventh_cyto = c(0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 1, 0, 1), time = c(4, 26, 28, 53, 23, 23, 2,
3, 22, 3, 3, 33, 3), status = c(0, 0, 0, 0, 0, 0, 1, 1, 1,
1, 1, 1, 1)), row.names = c("160088-T",
"S5", "S7", "S21", "S27", "S31", "160098-T",
"160084-T",
"160093-T",
"160095-T",
"S9", "S13", "S17"), class = "data.frame")


where I have made the following models from:

model.full = glm(status ~ first_cyto + second_cyto + third_cyto + fourth_cyto + fifth_cyto + sixth_cyto + seventh_cyto + ninth_cyto + tenth_cyto + eleventh_cyto,
data=coxy,
)


and the null:

model.null = glm(status ~ 1,
data=coxy,
)


And I have calculated the deviance between the null hypothesis model and the actual model using a Chi-squared test:

anova(model.full,
model.null,
test="Chisq")


Which suggested deviance between the two models, albeit not significant:

Analysis of Deviance Table

Model 1: status ~ first_cyto + second_cyto + third_cyto + fourth_cyto +
fifth_cyto + sixth_cyto + seventh_cyto + ninth_cyto + tenth_cyto +
eleventh_cyto
Model 2: status ~ 1
Resid. Df Resid. Dev  Df Deviance Pr(>Chi)
1         2      0.000
2        12     17.945 -10  -17.945   0.0559 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


However when looking at the actual summary of model.full:

> summary(model.full)

Call:
glm(formula = status ~ first_cyto + second_cyto + third_cyto +
fourth_cyto + fifth_cyto + sixth_cyto + seventh_cyto + ninth_cyto +
tenth_cyto + eleventh_cyto, family = binomial(link = "logit"),
data = coxy)

Deviance Residuals:
160088-T          S5          S7         S21         S27         S31
-6.547e-06  -6.547e-06  -6.547e-06  -6.547e-06  -6.547e-06  -6.547e-06
160098-T    160084-T    160093-T    160095-T          S9         S13
6.547e-06   6.547e-06   6.547e-06   6.547e-06   6.547e-06   6.547e-06
S17
6.547e-06

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept)   -2.457e+01  1.155e+05       0        1
first_cyto     4.913e+01  1.610e+05       0        1
second_cyto    7.866e-07  1.105e+05       0        1
third_cyto    -4.913e+01  2.210e+05       0        1
fourth_cyto    7.866e-07  1.204e+05       0        1
fifth_cyto     4.913e+01  2.531e+05       0        1
sixth_cyto     1.512e-09  1.725e+05       0        1
seventh_cyto   3.697e-07  1.610e+05       0        1
ninth_cyto    -4.913e+01  1.538e+05       0        1
tenth_cyto     4.913e+01  3.209e+05       0        1
eleventh_cyto  4.913e+01  1.381e+05       0        1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 1.7945e+01  on 12  degrees of freedom
Residual deviance: 5.5727e-10  on  2  degrees of freedom
AIC: 22

Number of Fisher Scoring iterations: 23


None of the variables are significant, which I was expecting as there just aren't enough samples (N=13). But according to a previous post, there is still separation as there is a high Fisher score (23) and the positive coefficients are high with large standard errors. Does it still make sense to then use predict() function on the data? And complete a confusion matrix?

coxy$predy = predict(model.full, type="response")  where it appears coxy$predy is correctly assigning each into the correct outcome group?

(dput of coxy$predy) coxy$predy <- c(8.94455269520935e-12, 2.93370332953973e-11, 2.22044604925031e-16,
1.27189451620771e-11, 2.22044604925031e-16, 1.36197070333109e-11,
0.999999999997593, 0.999999999996578, 0.999999999976139, 1, 0.999999999986539,
0.999999999997099, 0.999999999976139)


So, my N is small, but there is a deviance between the null and the full model, the estimates and fisher score are high and the prediction of outcome assignment looks right. Can I do this? Does this make sense at all?