Multiple Nested Random Effects affecting a Mixed Model

Question

I have a quite complex psychophysiological data dependant of different nested data in a repeated measures experiment.

The first nested structure comes from the data collection were there are several blocks containing several trials with 3 repetitions per trial. The second comes from grouping as we are measuring the adjustment in a particular task between a pair of participants in this repetitions of the same trial.

I'd like to do various measurements within pair of subjects but I always end up having the same constraint with the nested data.

The result of not having a good structured data makes the models I try to fit have enormous amounts of df wich consequently make me having irrelevant models (that explain only a ridiculous amount of variability).

Logic tells me it should be nested by random effects twice:

(1|CoupleID/SubjectID) <- Regarding the nesting of the subjects
(1|Block/Trial/TrialRep) <- Regarding at what point in the experiment was the response data collected

Therefore, this would be some of the models I would be interested in to be served as an example, so I can relate a particular response to a behavioural measurement:

1. lmer(electrode_response~TrialRep)
2. lmer(electrode_response~divergence_withinparticipants)
3. lmer(electrode_response~participant_adjustment)

The random effects term I am trying looks like this:

lmer(electrode_response~TrialRep+(TrialRep |SubjectID)+(1|CoupleID/SubjectID)+(1|Block/Trial/TrialRep)

Which a priori could look like it gives a good explanation:

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: electrode_response ~ TrialRep + (TrialRep | SubjectID) + (1 | CoupleID/SubjectID) + (1 | Block/Trial/TrialRep)
   Data: finaldb
Control: lmerControl(optimizer = "optimx", optCtrl = list(method = "nlminb"))

REML criterion at convergence: 80726

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-5.3489 -0.5861  0.0203  0.6204  6.7150 

Random effects:
 Groups                 Name        Variance  Std.Dev. Corr 
 TrialRep:(Trial:Block) (Intercept)   0.82127  0.9062       
 Trial:Block            (Intercept)   0.05282  0.2298       
 SubjectID:CoupleID     (Intercept)   9.21819  3.0361       
 SubjectID              (Intercept)  19.82985  4.4531       
                          TrialRep    3.16028  1.7777  -1.00
 CoupleID               (Intercept)   2.03290  1.4258       
 Block                  (Intercept)   0.08479  0.2912       
 Residual                           100.53917 10.0269       
Number of obs: 10800, groups:  TrialRep:(Trial:Block), 300; Trial:Block, 100; SubjectID:CoupleID, 36; SubjectID, 36; CoupleID, 18; Block, 4

Fixed effects:
            Estimate Std. Error      df t value Pr(>|t|)    
(Intercept)   8.5483     0.9857 24.1529   8.672 6.98e-09 ***
TrialRep     -1.8157     0.3254 37.6691  -5.581 2.20e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
         (Intr)

As you can see I'm using an optimizer, but even with this I often, in the data collected in some electrodes, can find the subsequent warnings, even though I use the same nesting structure in the data:

boundary (singular) fit: see ?isSingular
Warning message:
In optwrap(optimizer, devfun, getStart(start, rho$lower, rho$pp),  :
  convergence code 1 from optimx

Or

Warning messages:
1: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
  unable to evaluate scaled gradient
2: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
  Model failed to converge: degenerate  Hessian with 1 negative eigenvalues
3: Model failed to converge with 1 negative eigenvalue: -7.8e-01 

Wich I fail to comprehend. Also, the Rsq I get with r.squaredGLMM() of MuMIn package tells me not a lot is explained by my model.

> r.squaredGLMM(electrode_response)
            R2m       R2c
[1,] 0.01769641 0.1905829

I would really appreciate some feedback on this matter, I have gone out of ideas.