Unstructured error covariance matrix in a multilevel growth model

Yep, I can confirm this method produces the desired variance/covariance method, but it loses the desired linear model. Here are the results:

# load 
library(tidyverse)

# download and wrangle the data
data1 <- read.table("https://stats.idre.ucla.edu/stat/r/examples/alda/data/opposites_pp.txt", header = TRUE, sep = ",")
data2 <- subset(data1, select = -c(wave))
data3 <- spread(data2, time, opp)
colnames(data3) <- c("id", "cog", "ccog", "t0", "t1", "t2", "t3")

# fit the model
fit.unstructured <-
  brm(data = data3, 
      family = gaussian,
      bf(mvbind(t0, t1, t2, t3) ~ ccog) + 
        set_rescor(TRUE))
 Family: MV(gaussian, gaussian, gaussian, gaussian) 
  Links: mu = identity; sigma = identity
         mu = identity; sigma = identity
         mu = identity; sigma = identity
         mu = identity; sigma = identity 
Formula: t0 ~ ccog 
         t1 ~ ccog 
         t2 ~ ccog 
         t3 ~ ccog 
   Data: data3 (Number of observations: 35) 
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup samples = 4000

Population-Level Effects: 
             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
t0_Intercept   164.50      5.99   152.54   176.02 1.00     2626     2525
t1_Intercept   192.15      5.32   181.51   201.90 1.00     2315     2526
t2_Intercept   216.84      5.50   206.03   227.29 1.00     2380     2851
t3_Intercept   246.27      5.66   235.30   257.53 1.00     2787     2737
t0_ccog         -0.01      0.50    -0.99     1.02 1.00     2704     3082
t1_ccog          0.17      0.45    -0.73     1.07 1.00     2314     2921
t2_ccog          0.72      0.47    -0.22     1.63 1.00     2148     2562
t3_ccog          1.25      0.48     0.30     2.16 1.00     2675     2492

Family Specific Parameters: 
         Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma_t0    36.26      4.31    29.05    45.71 1.00     3077     2876
sigma_t1    32.46      3.64    26.04    40.25 1.00     2526     2603
sigma_t2    33.66      3.82    26.90    41.97 1.00     2596     2893
sigma_t3    34.47      4.09    27.53    43.35 1.00     3142     2892

Residual Correlations: 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
rescor(t0,t1)     0.74      0.08     0.56     0.86 1.00     2638     3053
rescor(t0,t2)     0.65      0.10     0.43     0.80 1.00     2606     2792
rescor(t1,t2)     0.79      0.07     0.63     0.89 1.00     2898     3058
rescor(t0,t3)     0.33      0.14     0.02     0.58 1.00     2493     2706
rescor(t1,t3)     0.63      0.10     0.40     0.79 1.00     3344     3004
rescor(t2,t3)     0.72      0.08     0.52     0.85 1.00     3499     2813

Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

Who knows how to fit this?

\begin{align*} \text{opp}_{ij} & \sim \operatorname{Normal}(\mu_{ij}, \mathbf r) \\ \mu_{ij} & = \gamma_{00} + \gamma_{10} \text{time}_{ij} + \gamma_{01} (\text{cog}_i - \overline{\text{cog}}) + \gamma_{11} (\text{cog}_i - \overline{\text{cog}}) \times \text{time}_{ij} \\ \mathbf r & = \begin{bmatrix} \mathbf{\Sigma_r} & \mathbf 0 & \mathbf 0 & \dots & \mathbf 0 \\ \mathbf 0 & \mathbf{\Sigma_r} & \mathbf 0 & \dots & \mathbf 0 \\ \mathbf 0 & \mathbf 0 & \mathbf{\Sigma_r} & \dots & \mathbf 0 \\ \vdots & \vdots & \vdots & \ddots & \mathbf 0 \\ \mathbf 0 & \mathbf 0 & \mathbf 0 & \mathbf 0 & \mathbf{\Sigma_r} \end{bmatrix}, \text{where} \\ \mathbf{\Sigma_r} & = \begin{bmatrix} \sigma_1^2 & \sigma_{12} & \sigma_{13} & \sigma_{14} \\ \sigma_{21} & \sigma_2^2 & \sigma_{23} & \sigma_{24} \\ \sigma_{31} & \sigma_{32} & \sigma_3^2 & \sigma_{34} \\ \sigma_{41} & \sigma_{42} & \sigma_{43} & \sigma_4^2 \end{bmatrix}. \end{align*}