Extracting standard error in R for estimated time series
Opened this issue · 7 comments
I am having trouble getting std. error out of the TMB code. I am using ADREPORT(log_x_t); for the estimated time series, but can't seem to access std. error for it in R.
What do you mean exactly?
Typically a SE for a derived quantity can be obtained using summary(SD, "report") where SD is generated by sdreport. If it is not there, then it is likely a bug e.g., where you added the ADREPORT statement but haven't recompiled the dll because it is currently loaded.
Please try playing around a bit and then feel free to send me a minimal example if you're still wanting help!
If you are using the TMBhelper:optimize function you can do the following:
Obj = MakeADFun( data=DataList, parameters=Params, DLL="model", random = random, map = map)
Opt = TMBhelper::Optimize( obj=Obj, newtonsteps=1)
sd <- sdreport(Obj)
summary(sd)
Ok, I've got that now as:
SD_rep<-summary(sdreport(Obj))
producing:
Estimate Std. Error
log_d0 0.000000 3.8297150
alpha 1.000000 0.5509031
rho 1.000000 0.2787609
.
.
.
log_x_t 1.000000 4.6482617
log_x_t 1.000000 6.0601868
log_x_t 1.000000 7.2157557
log_x_t 1.000000 8.1442962
log_x_t 1.000000 8.8701084
How can I extract the last 5 log_x_t values (std. error), so I can save them without hard-coding the index number of where they are in the report?
Another option is to access the fitted values and their variance as
SD_rep$par.random
SD_rep$diag.cov.randomThen take the square root to get the standard error.
kinda clunky for saving, but works
SD_rep <- SD_rep[which(row.names(SD_rep) == "log_x_t"),]
SD_rep_save <- SD_rep[c((nrow(SD_rep) - 4) : nrow(SD_rep)),]
Ended up using:
se.fore[i,]<- SD_rep[(dim(SD_rep)[1]-forecast+1):dim(SD_rep)[1],2]
Jim: would love to learn how to use grep. Maybe after class this week?