brainchart / Lifespan

Code used in the lifespan brainchart project
https://www.nature.com/articles/s41586-022-04554-y
Other
62 stars 14 forks source link

Errors running the tutorial R - 4.3.1 #5

Open richardbeare opened 5 months ago

richardbeare commented 5 months ago

Hi, I'm working through the tutorial and encounter the following. I must be doing something stupid, but can't see it.

# 
source("220.simulation-omega-setup.r")

source("100.common-variables.r")
source("101.common-functions.r")

source("300.variables.r")
source("301.functions.r")

## 310-script
PATHS <- Create.Folders( "omega-Wand__.n0000" )
HOLDER <- Load.Subset.Wrapper( Tag="omega-Wand__.n0000", LSubset=TRUE )
HOLDER$MODEL <- readRDS( file.path( PATHS$MODEL, "base200.GGalt.fp.rds" ) )
EXTRACT <- Extract.Wrapper( HOLDER, Store.Full=TRUE ) ## store full fitting object to use as initial point of bfpNA() re-fit [expect 12 iterations]
summary(EXTRACT$FIT.FULL) ## standard summary of a GAMLSS fit
getSmo(EXTRACT$FIT.FULL,what="mu")
getSmo(EXTRACT$FIT.FULL,what="mu")$power

HOLDER$MODEL <- Make.bfpNA.model.from.extract( EXTRACT$param )
saveRDS( HOLDER$MODEL, file.path( PATHS$MODEL, "base200.GGalt.bfpNA.rds" ) )
EXTRACT.bfp <- Extract.Wrapper( HOLDER, Fit.Full=FALSE, start.from=EXTRACT$FIT ) ## helpful to start.from, improves convergence speed [expect 5 iterations]
Save.Extracted( EXTRACT.bfp, PATHS, "base200.GGalt.bfpNA.rds", Save.Full=FALSE )

## 310-script
EXTRACT.bfp$param$BIC ## compare BIC on all fitted models

file.copy(from=file.path(PATHS$FIT.EXTRACT,"base200.GGalt.bfpNA.rds"),to=file.path(PATHS$PATH,"MODEL.rds"))
file.copy(from=file.path(PATHS$FIT.EXTRACT,"base200.GGalt.bfpNA.rds"),to=file.path(PATHS$PATH,"FIT.EXTRACT.rds"))

HOLDER <- Load.Subset.Wrapper( Tag="omega-Wand__.n0000", LSubset=TRUE, LModel=TRUE, LFit=TRUE )
Error in get(OBJ$MODEL$family) : invalid first argument

This line is accessing MODEL$family, but seems it should access MODEL$params$family. Changing that line gets things along a bit further, but runs into problems at the boot stage:

BOOT <- list()
BOOT[[1]] <- Boot.Function(n=1,Base.Seed=12345,Holder=HOLDER)
BOOT[[1]]
$diagnostics
[1] "try-error"

$iter
[1] -1

$base
[1] 12345

$offset
[1] 1

Setting a breakpoint inside gamlssWrapper and examining Model shows that there's no Model$mu - there is a Model$param$mu, which is a larger object, rather than a formula.

What have I missed?

richardbeare commented 5 months ago

I think I've sorted out the problem - the copy command in the README was duplicated so that MODEL.rds and FIT.EXTRACT.rds had the same content. See PR #6