random effects in dispformula

[mebrooks]

This isn't quite ready to merge, but I wanted some feedback. The code compiles, but vignettes and tests fail. I think part of the problem comes from the files saved in
file_ok <- gt_load("test_data/models.rda", update_gauss_disp = TRUE). How are these updated and should that happen in the Makefile?

Here are a list of potential issues

• In glmmTMB.R, line 1169-1174 does something to ziformula. Not sure if it should be done to dispformula.
• Why is disp treated as other in confint.glmmTMB()?
• I did not create an equivalent of ziPredictCode.
• I didn’t do anything with emmeans.

Some changes were to make our handling of "zi" and "disp" more consistent because some objects had "d" in the name instead of "disp".
etad -> etadisp
Zd -> Zdisp
Xd -> Xdisp
XdS -> XdispS
doffset -> dispoffset
sparseXd -> sparseXdisp

[bbolker]

cd glmmTMB/inst/test_data
touch make_ex.R  ## to force remaking
make

working on this now ...

[mebrooks]

I'm not sure how to solve this problem. The point of some tests is to test back compatibility ('test-sparseX.R:22:5', 'test-methods.R:161:5'), but just reading in the oldfit gives an error.

readRDS(system.file("test_data", "oldfit.rds",
+                               package="glmmTMB",
+                               mustWork=TRUE))
Formula:          log.r ~ epoch + source + (1 | outbreak.year)
Dispersion:             ~1 + offset(log(sdvals^2))
Data: d
     AIC      BIC   logLik df.resid 
 9.91377 14.54931  1.04311       10 
Error in .Call("MakeDoubleFunObject", data, parameters, reportenv, PACKAGE = DLL) : 
  Incorrect number of arguments (3), expecting 4 for 'MakeDoubleFunObject'

[bbolker]

If the internal structure of glmmTMB objects has changed we may have to modify/update up2date ... ? (Or if we're lucky we just have to run up2date on the object?)

[mebrooks]

In up2date, if a model didn't have Zdisp in its data before, what do we want to set it to when we create it? Is it like this line in getXReTrms

Z <- new("dgCMatrix",Dim=c(as.integer(nobs),0L))

[bbolker]

That seems right ... ?

[mebrooks]

This version after the push still gives the error

Error in getParameterOrder(data, parameters, new.env(), DLL = DLL) : 
  Error when reading the variable: 'Zdisp'. Please check data and parameters.
In addition: Warning message:
In getParameterOrder(data, parameters, new.env(), DLL = DLL) :
  Expected sparse matrix of class 'dgTMatrix'.

If I do

oldfit <- readRDS(system.file("test_data", "oldfit.rds",
                           package="glmmTMB",
                            mustWork=TRUE))

and then go line by line in up2date

> with(ee,
+                        TMB::MakeADFun(data,
+                                       parameters,
+                                       map = map,
+                                       random = random,
+                                       silent = silent,
+                                       DLL = "glmmTMB"))
Error in getParameterOrder(data, parameters, new.env(), DLL = DLL) : 
  Error when reading the variable: 'Zdisp'. Please check data and parameters.
In addition: Warning message:
In getParameterOrder(data, parameters, new.env(), DLL = DLL) :
  Expected sparse matrix of class 'dgTMatrix'.
> ee$data$Zdisp
16 x 0 sparse Matrix of class "dgCMatrix"
     
 [1,]
 [2,]
 [3,]
 [4,]
 [5,]
 [6,]
 [7,]
 [8,]
 [9,]
[10,]
[11,]
[12,]
[13,]
[14,]
[15,]
[16,]
> ee$random
[1]  5  6  7  8  9 10 11 12 13

I wonder why it expected dgTMatrix and not dgCMatrix.

[mebrooks]

Using dgTMatrix instead of dgCMatrix fixed that error. I'm not sure why, but I'll keep working on the other errors.

[mebrooks]

This is getting closer to complete, but I need help figuring out the last 3 failing tests.

── Failure (test-methods.R:371:5): profile (no RE) ─────────────────────────────
dim(p0_th) not equal to c(41, 3).
1/2 mismatches
[1] 43 - 41 == 2

> table(p0_th$.par)

     (Intercept)             Days disp~(Intercept) 
              14               14               15 

Error ('test-predict.R:522:5'): weights with attributes are OK
Error in `MakeADFun(data.tmb, parameters, map = mapArg, random = randomArg, 
    profile = NULL, silent = TRUE, DLL = "glmmTMB")`: Only numeric matrices, vectors, arrays, factors, lists or length-1-characters can be interfaced

3. There's a chance that a change I made caused a problem with the estimates in 'test-varstruc.R:44:5' unless this was ill-conditioned from the beginning.

Failure ('test-varstruc.R:44:5'): print ar1 (>1 RE)
cco[12:14] not equal to c(...).
3/3 mismatches
x[1]: "Subject (Intercept) 8e-05 9e-03"
y[1]: "Subject (Intercept) 4e-01 0.6"

x[2]: "Subject.1 row1 4e+03 6e+01 0.87 (ar1)"
y[2]: "Subject.1 row1 4e+03 60.8 0.87 (ar1)"

x[3]: "Residual 8e+01 9e+00"
y[3]: "Residual 8e+01 8.9"

[bbolker]

── Failure (test-methods.R:371:5): profile (no RE) ─────────────────────────────
dim(p0_th) not equal to c(41, 3).
1/2 mismatches
[1] 43 - 41 == 2

> table(p0_th$.par)

     (Intercept)             Days disp~(Intercept) 
              14               14               15 

I'm having a hard time with this one too. When I run it with the current devel version I get the table to be

  (Intercept)          Days d~(Intercept) 
           14            14             6 

(i.e. the number of rows is 34, not 41, and the test fails anyway). Comparing the old (current master branch) and new (dispRE) results:

In this case the old profile is clearly messed up by not having the dispersion parameter corrected to the new scale (i.e. when we changed the definition of dispersion for Gaussian models). I think that some combination of (1) versions of TMB? and (2) scaling of dispersion may account for this difference (i.e., rescaling the dispersion by 2 will change the profile calculation slightly). To be honest, this isn't a great test anyway and I would be fine with updating it to match the current version. (Also wondering why we didn't catch this sooner ... maybe we don't run with NOT_CRAN=true often enough and some tests got skipped?)

Error ('test-predict.R:522:5'): weights with attributes are OK
Error in `MakeADFun(data.tmb, parameters, map = mapArg, random = randomArg, 
    profile = NULL, silent = TRUE, DLL = "glmmTMB")`: Only numeric matrices, vectors, arrays, factors, lists or length-1-characters can be interfaced

We've hit stuff like this before. TMB is fussy about attributes, and I fixed a similar problem recently in 726a705. Not sure how this one slipped by

3. There's a chance that a change I made caused a problem with the estimates in `'test-varstruc.R:44:5' ` unless this was ill-conditioned from the beginning.

Failure ('test-varstruc.R:44:5'): print ar1 (>1 RE)
cco[12:14] not equal to c(...).
3/3 mismatches
x[1]: "Subject (Intercept) 8e-05 9e-03"
y[1]: "Subject (Intercept) 4e-01 0.6"

x[2]: "Subject.1 row1 4e+03 6e+01 0.87 (ar1)"
y[2]: "Subject.1 row1 4e+03 60.8 0.87 (ar1)"

x[3]: "Residual 8e+01 9e+00"
y[3]: "Residual 8e+01 8.9"

Hmm. The new version (first line) is a singular fit for the nugget effect ("Subject (Intercept)") where it was small but non-singular before. I'm out of time for right now, I might try some git bisecting tomorrow to see if I can figure out where this happened ...

[bbolker]

I think problem number 3 dates back to the Gaussian dispersion scale change and wasn't caught before now because we didn't update stored fits. git bisect (based on misc/test-ar1var.R) says it happened at 5939686, but I suspect it might be one commit prior to that ... will work on it more tomorrow.

[bbolker]

example 3 is definitely an unstable fit, and the test definitely broke when we started changing the Gaussian dispersion scale (we should make sure to highlight such possible changes in the NEWS for the Gaussian-dispersion change ...). I think this basically means that both of the remaining test failures are false positives, we just have to decide what to do about them ...