utils::globalVariables(c("name","age"))
utils::globalVariables(c("ddply"))
utils::globalVariables(c("."))
utils::globalVariables(c("rstudent"))
utils::globalVariables(c("lm"))
utils::globalVariables(c("rstandard"))
utils::globalVariables(c("read.table"))
utils::globalVariables(c("diagsFn"))
readSam<-function(file, reduced=FALSE){
# Function to read a basic fit
ret<-list()
parfile<-as.numeric(scan(paste(file,'.par', sep=''),
what='', n=16, quiet=TRUE)[c(6,11,16)])
ret$nopar<-as.integer(parfile[1])
ret$nlogl<-parfile[2]
ret$maxgrad<-parfile[3]
rep<-scan(paste(file,'.rep', sep=''), quiet=TRUE)
ret$res<-read.table(paste(file,'.res', sep=''),header=FALSE)
ret$stateDim<-rep[1]
ret$years<-rep[-1]
file<-paste(file,'.cor', sep='')
lin<-readLines(file)
ret$npar<-length(lin)-2
ret$logDetHess<-as.numeric(strsplit(lin[1], '=')[[1]][2])
sublin<-lapply(strsplit(lin[1:ret$npar+2], ' '),function(x)x[x!=''])
ret$names<-unlist(lapply(sublin,function(x)x[2]))
ret$est<-as.numeric(unlist(lapply(sublin,function(x)x[3])))
ret$std<-as.numeric(unlist(lapply(sublin,function(x)x[4])))
ret$cor<-matrix(NA, ret$npar, ret$npar)
for(i in 1:ret$npar){
ret$cor[1:i,i]<-as.numeric(unlist(lapply(sublin[i],
function(x)x[5:(4+i)])))
ret$cor[i,1:i]<-ret$cor[1:i,i]
}
ret$cov<-ret$cor*(ret$std%o%ret$std)
mslh<-function(name){
idx<-which(ret$names==name)
x<-cbind(ret$est[idx], ret$std[idx], ret$est[idx]-2*ret$std[idx],
ret$est[idx]+2*ret$std[idx])
colnames(x)<-c('est', 'std', 'low', 'hig')
return(x)
}
ret$ssb<-mslh('ssb')
ret$fbar<-mslh('fbar')
ret$tsb<-mslh('tsb')
ret$logssb<-mslh('logssb')
ret$logfbar<-mslh('logfbar')
ret$logtsb<-mslh('logtsb')
ret$logscale<-mslh('logScale')
ret$logFpar<-mslh('logFpar')
ret$logCatch<-mslh('logCatch')
x<-mslh('U')
ret$stateEst<-matrix(x[,1],ncol=ret$stateDim, byrow=TRUE)
ret$stateStd<-matrix(x[,2],ncol=ret$stateDim, byrow=TRUE)
ret$stateLow<-matrix(x[,3],ncol=ret$stateDim, byrow=TRUE)
ret$stateHig<-matrix(x[,4],ncol=ret$stateDim, byrow=TRUE)
ret$R<-cbind(exp(ret$stateEst[,1]), NA, exp(ret$stateLow[,1]),
exp(ret$stateHig[,1]))
if(reduced){
ret<-ret[which(!names(ret)%in%c('cov','cor'))]
}
return(ret)
}
utils::globalVariables(c("index"))
.diagUsam=function(x){
sam =readSam(x)
res =sam$res
names(res)=c("year","name","age","obs","hat","residual")
res=transform(res,roriginal=residual,
residual =residual*exp(sam$est[sam$names=="logSdLogObs"])[name])
res=ddply(res,.(name,age), function(x) transform(x, rstudent =rstudent( lm(obs~0+hat)),
rstandard=rstandard(lm(obs~0+hat))))
res=ddply(res,.(name,age), diagsFn)
names(res)[seq(length(names(res)))[names(res)=="obs"]]="index"
res}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.