How to interpret GLMM results?

My question is related with my previous post Extract variance of the fixed effect in a glmm.
However, in this case I change the model that the GLMM follow.
It follows a log family and as there are many zeros in my dataset,so I used a zero inflation method.

I would like to get the variation (variance component) in incidence (inc.) within each habitat while being mindful of random factors such as season and site

This is my data set:

## Incidence:

Incidence <- data.frame(Inc. = c(0.4400, 0.5102, 0.2979, 0.2667, 0.0000, 0.0000,
                                 0.0200, 0.0213, 0.0000, 0.0238, 0.0256, 0.0000,
                                 0.0000, 0.1538, 0.0417, 0.0000, 0.0734, 0.0000,
                                 0.0000, 0.0000, 0.1293, 0.0072, 0.0000, 0.0078,
                                 0.0000, 0.0000, 0.0000, 0.0068, 0.0000, 0.0000,
                                 0.0068), 
                        Habitat = c("Crop", "Crop", "Crop", "Crop", "Edge", "Edge", 
                                    "Edge", "Edge", "Edge", "Edge", "Edge", "Edge", 
                                    "Edge", "Edge", "Edge", "Oakwood", "Oakwood", 
                                    "Oakwood", "Oakwood", "Oakwood", "Oakwood", 
                                    "Oakwood", "Oakwood", "Wasteland", "Wasteland", 
                                    "Wasteland", "Wasteland", "Wasteland", "Wasteland", 
                                    "Wasteland", "Wasteland"), 
                        Season = c("Summer", "Summer", "Summer", "Summer", "Autumn", 
                                   "Autumn", "Autumn", "Autumn", "Spring", "Spring", 
                                   "Spring", "Spring", "Summer", "Summer", "Summer", 
                                   "Autumn", "Autumn", "Autumn", "Autumn", "Spring",
                                   "Spring", "Spring", "Spring", "Autumn", "Autumn", 
                                   "Autumn", "Autumn", "Spring", "Spring", "Spring", 
                                   "Spring"), 
                        Site = c("M1", "M2", "M3", "M4", "L1", "L2", "L3", "L4", 
                                 "L1", "L2", "L3", "L4", "L1", "L2", "L3", "Q1", 
                                 "Q2", "Q3", "Q4", "Q1", "Q2", "Q3", "Q4", "E1", 
                                 "E2", "E3", "E4", "E1", "E2", "E3", "E4"))

With the aim to get the variation I check previously with a shapiro wilk test how is the distribution of my dataset by Rstudio.

shapiro.test(x = Incidence$Inc.):

       Shapiro-Wilk normality test
       data:  Incidence$Incidence
       W = 0.56708, p-value = 2.092e-08

Moreover I got the homocedasticity with a levene test:

leveneTest(y = Incidence$Inc., group = Incidence$Habitat, center = "median")

     Levene's Test for Homogeneity of Variance (center = "median")
           Df F value   Pr(>F)   
     group  3  6.3481 0.002129 **
     27                    
     ---
     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Afterward I check how is the distribution using:

Input_2 <- Incidence$Inc.
library(rriskDistributions)
Prueba <- fit.cont(as.vector(t(Input_2)))

and I got a log distribution

Then I performed a glmm of this dataset in R:

GlM_habitats <- glmmTMB(Inc.~ Habitat + (1|Season)+ (1|Site),
                        data = Incidence,
                        ziformula = ~1,
                        family = poisson(link = "log")) 
 
#Warning message:
#In glmmTMB(Inc.~ Habitat + (1 | Season) + (1 | Site), data = Incidence,  :
#non-integer counts in a poisson model

Anova(GlM_habitats)
Analysis of Deviance Table (Type II Wald chisquare tests)

Response: Incidence
         Chisq Df Pr(>Chisq)
Habitat 3.0632  3      0.382 

summary(GlM_habitats)
   Family: poisson  ( log )
   Formula:          Inc.~ Habitat + (1 | Season) + (1 | Site)
   Zero inflation:             ~1
   Data: Incidence

    AIC      BIC   logLik deviance df.resid 
    23.5     33.5     -4.7      9.5       24 

 Random effects:

 Conditional model:
  Groups Name        Variance  Std.Dev.
  Season (Intercept) 5.656e-13 7.52e-07
  Site   (Intercept) 1.176e-13 3.43e-07
 Number of obs: 31, groups:  Season, 3; Site, 16

 Conditional model:
                  Estimate Std. Error z value Pr(>|z|)
 (Intercept)       -0.9710     0.8125  -1.195    0.232
 HabitatEdge       -2.6780     2.0382  -1.314    0.189
 HabitatOakwood    -2.6696     2.3290  -1.146    0.252
 HabitatWasteland  -4.9528     6.8841  -0.720    0.472

 Zero-inflation model:
             Estimate Std. Error z value Pr(>|z|)
 (Intercept)    -24.1    43216.9  -0.001        1

Then as in previous post they anwered to me I tried to extract the variance of fixed effect:

 # Variance of random effects: 
 vc <- lme4::VarCorr(GlM_habitats)
 print(vc,comp=c("Variance","Std.Dev."),digits=2)

 Conditional model:
  Groups Name        Variance Std.Dev.
  Season (Intercept) 5.7e-13  7.5e-07 
  Site   (Intercept) 1.2e-13  3.4e-07 

  # Variance-Covariance Matrix of fixed effects: 
  vc_fixed <- as.matrix(vcov(GlM_habitats))
  # Variance of fixed effects: 
  var_fixed <- diag(vc_fixed); var_fixed

  [[1]]
                   (Intercept) HabitatEdge HabitatOakwood HabitatWasteland
  (Intercept)         0.660153   -0.660153      -0.660153        -0.660153
  HabitatEdge        -0.660153    4.154245       0.660153         0.660153
  HabitatOakwood     -0.660153    0.660153       5.424338         0.660153
  HabitatWasteland   -0.660153    0.660153       0.660153        47.390362



  # Standard errors of fixed effects: 
  se_fixed <- sqrt(var_fixed); se_fixed

 

When I perform this analysis I got this

 Error in sqrt(var_fixed) : non-numeric argument to mathematical function

  

I would like to know how to interpret this result and to know if they have been performed OK. I can not believe that Season and Site have very low variance and the ANOVA results give a p value that is not significant. Moreover, I do not know why the Standard errors of fixed effects do not work.

What Am I doing wrong?

> m0 <- lmer(Inc.~ Habitat + (1|Season)+ (1|Site),
+            data = Incidence)
> summary(m0)
Linear mixed model fit by REML ['lmerMod']
Formula: Inc. ~ Habitat + (1 | Season) + (1 | Site)
   Data: Incidence

REML criterion at convergence: -78.9

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.45229 -0.30319 -0.01575  0.20558  2.53994 

Random effects:
 Groups   Name        Variance  Std.Dev.
 Site     (Intercept) 0.0031294 0.05594 
 Season   (Intercept) 0.0005702 0.02388 
 Residual             0.0008246 0.02872 
Number of obs: 31, groups:  Site, 16; Season, 3

Fixed effects:
                 Estimate Std. Error t value
(Intercept)       0.35450    0.03607   9.827
HabitatEdge      -0.32669    0.04475  -7.301
HabitatOakwood   -0.31616    0.04637  -6.818
HabitatWasteland -0.33973    0.04637  -7.326

hist(residuals(m0))