Nothing
hmscHT.HMSCdata <-
function(data,start,priors=NULL,niter=12000,nburn=2000,nrot=100,thin=1,likelihoodtype="logit",details=TRUE,verbose=TRUE,...){
#### F. Guillaume Blanchet - July 2013
##########################################################################################
call<-match.call()
### Basic objects
nparamX<-ncol(start$paramX)
nsp<-nrow(start$paramX)
nparamTr<-ncol(start$paramTr)
nexp<-nrow(start$paramTr)
if(is.null(start$means)){
start$means<-matrix(NA,ncol=nsp,nrow=nexp)
for(i in 1:nsp){
start$means[,i]<-tcrossprod(data$Tr[,i],start$paramTr)
}
}
### General checks
if(niter < nburn){
stop("'niter' should be equal or larger than 'burning'")
}
if(is.null(priors)){
if(!any(class(data)=="HMSCdataRandom")){
priors<-list(means0=rep(0,nparamX),kappa0=0,nu0=rep(0,nparamX),Lambda0=diag(nparamX))
}else{
priors<-list(means0=rep(0,nparamX),kappa0=0,nu0=rep(0,nparamX),Lambda0=diag(nparamX),priorRandom=c(-1,1))
nparamRandom<-ncol(start$paramRandom)
}
}
#=====================================================================================================
#### Construct the model using Single Component Adaptative Metropolis (SCAM) with eigenvector rotation
#=====================================================================================================
#--------------------
#### Starting objects
#--------------------
#### Rotation objects (to use when no rotation is carried out)
if(!any(class(data)=="HMSCdataRandom")){
rotRes<-vector("list",length=4)
names(rotRes)<-c("values","vectors","wParamX","wOuterParamX")
rotRes[[1]]<-matrix(1,nrow=nsp,ncol=nparamX)
rotRes[[2]]<-array(diag(nparamX),dim=c(nparamX,nparamX,nsp))
rotRes[[3]]<-matrix(0,nrow=nsp,ncol=nparamX)
rotRes[[4]]<-array(0,dim=c(nparamX,nparamX,nsp))
}else{
rotRes<-vector("list",length=9)
names(rotRes)<-c("values","vectors","wParamX","wOuterParamX","valuesRandom","vectorsRandom","wParamRandom","wOuterParamRandom","wRandomCov")
rotRes[[1]]<-matrix(1,nrow=nsp,ncol=nparamX)
rotRes[[2]]<-array(diag(nparamX),dim=c(nparamX,nparamX,nsp))
rotRes[[3]]<-matrix(0,nrow=nsp,ncol=nparamX)
rotRes[[4]]<-array(0,dim=c(nparamX,nparamX,nsp))
rotRes[[5]]<-matrix(1,nrow=nsp,ncol=nparamRandom)
rotRes[[6]]<-array(diag(nparamRandom),dim=c(nparamRandom,nparamRandom,nsp))
rotRes[[7]]<-matrix(0,nrow=nsp,ncol=nparamRandom)
rotRes[[8]]<-array(0,dim=c(nparamRandom,nparamRandom,nsp))
rotRes[[9]]<-1
}
if(!any(class(data)=="HMSCdataRandom")){
accept<-array(0,dim=c(nsp,nparamX,2))
sdvalue<-matrix(1,nrow=nsp,ncol=nparamX)
}else{
acceptparamX<-array(0,dim=c(nsp,nparamX,2))
acceptparamRandom<-array(0,dim=c(nsp,nparamRandom,2))
acceptRandomCov<-c(0,0)
accept<-list(accept=acceptparamX,acceptRandom=acceptparamRandom,acceptRandomCov=acceptRandomCov)
sdvalueParamX<-matrix(1,nrow=nsp,ncol=nparamX)
sdvalueParamRandom<-matrix(1,nrow=nsp,ncol=nparamRandom)
sdvalueRandomCov<-1
sdvalue<-list(sdvalue=sdvalueParamX,sdvalueRandom=sdvalueParamRandom,sdvalueRandomCov=sdvalueRandomCov)
}
scamRes<-list(start=start,rotation=rotRes,sdvalue=sdvalue,accept=accept)
#### Store burning results
if(!any(class(data)=="HMSCdataRandom")){
burningObj<-vector("list",length=4)
names(burningObj)<-c("paramX","paramTr","means","sigma")
burningObj[[1]]<-array(dim=c(nsp,nparamX,nburn))
burningObj[[2]]<-array(dim=c(nparamX,nparamTr,nburn))
burningObj[[3]]<-array(dim=c(nparamX,nsp,nburn))
burningObj[[4]]<-array(dim=c(nparamX,nparamX,nburn))
}else{
burningObj<-vector("list",length=6)
names(burningObj)<-c("paramX","paramTr","means","sigma","paramRandom","RandomVar")
burningObj[[1]]<-array(dim=c(nsp,nparamX,nburn))
burningObj[[2]]<-array(dim=c(nparamX,nparamTr,nburn))
burningObj[[3]]<-array(dim=c(nparamX,nsp,nburn))
burningObj[[4]]<-array(dim=c(nparamX,nparamX,nburn))
burningObj[[5]]<-array(dim=c(nsp,nparamRandom,nburn))
burningObj[[6]]<-vector("numeric",length=nburn)
}
#### Store detailed information
if(details){
if(!any(class(data)=="HMSCdataRandom")){
detailsObj<-vector("list",length=2)
names(detailsObj)<-c("sdvalue","acceptRatio")
detailsObj[[1]]<-array(dim=c(nsp,nparamX,nburn))
detailsObj[[1]][,,1]<-sdvalue
detailsObj[[2]]<-array(dim=c(nsp,nparamX,nburn))
detailsObj[[2]][,,1]<-matrix(0,nrow=nsp,ncol=nparamX)
}else{
detailsObj<-vector("list",length=6)
names(detailsObj)<-c("sdvalue","acceptRatio","acceptRatioRandomCov","acceptRatioRandom","sdvalueRandom","sdvalueRandomCov")
detailsObj[[1]]<-array(dim=c(nsp,nparamX,nburn))
detailsObj[[1]][,,1]<-sdvalueParamX
detailsObj[[2]]<-array(dim=c(nsp,nparamX,nburn))
detailsObj[[2]][,,1]<-matrix(0,nrow=nsp,ncol=nparamX)
detailsObj[[3]]<-numeric()
detailsObj[[3]][1]<-0
detailsObj[[4]]<-array(dim=c(nsp,nparamRandom,nburn))
detailsObj[[4]][,,1]<-matrix(0,nrow=nsp,ncol=nparamRandom)
detailsObj[[5]]<-array(dim=c(nsp,nparamRandom,nburn))
detailsObj[[5]][,,1]<-sdvalueParamRandom
detailsObj[[6]]<-numeric()
detailsObj[[6]][1]<-0
}
}
#### Store results after burning
if(!any(class(data)=="HMSCdataRandom")){
modelres<-vector("list",length=6)
names(modelres)<-c("paramX","paramTr","means","sigma","sdvalue","acceptRatio")
nvalues<-floor((niter-nburn)/thin)
modelres[[1]]<-array(dim=c(nsp,nparamX,nvalues))
modelres[[2]]<-array(dim=c(nparamX,nparamTr,nvalues))
modelres[[3]]<-array(dim=c(nparamX,nsp,nvalues))
modelres[[4]]<-array(dim=c(nparamX,nparamX,nvalues))
}else{
modelres<-vector("list",length=11)
names(modelres)<-c("paramX","paramTr","means","sigma","sdvalue","acceptRatio","paramRandom","sdvalueRandomCov","sdvalueRandom","acceptRatioRandomCov","acceptRatioRandom")
nvalues<-floor((niter-nburn)/thin)
modelres$paramX<-array(dim=c(nsp,nparamX,nvalues))
modelres$paramTr<-array(dim=c(nparamX,nparamTr,nvalues))
modelres$means<-array(dim=c(nparamX,nsp,nvalues))
modelres$sigma<-array(dim=c(nparamX,nparamX,nvalues))
modelres$paramRandom<-array(dim=c(nsp,nparamRandom,nvalues))
}
### Generate print counter
if(verbose){
printCounter<-round(seq(niter/5,niter,length=5))
}
#==================================
#### SCAM with eigenvector rotation
#==================================
itsave<-1
if(!any(class(data)=="HMSCdataRandom")){
for(i in 1:niter){
scamRes<-scamHT(data,scamRes$start,priors,scamRes$rotation$values,scamRes$rotation$vectors,scamRes$rotation$wParamX,scamRes$rotation$wOuterParamX,scamRes$sdvalue,scamRes$accept,likelihoodtype=likelihoodtype,burning=(i < nburn),rotation=(i > nrot & i < nburn),iter=i) # Pourrait changer
#-------------------------
#### Store burning results
#-------------------------
if(i < nburn){
burningObj[[1]][,,i]<-scamRes$start$paramX
burningObj[[2]][,,i]<-scamRes$start$paramTr
burningObj[[3]][,,i]<-scamRes$start$means
burningObj[[4]][,,i]<-scamRes$start$sigma
if(details){
detailsObj[[1]][,,i+1]<-scamRes$sdvalue
detailsObj[[2]][,,i+1]<-scamRes$accept[,,1]/scamRes$accept[,,2]
}
}else{
it<-i-nburn
if(it==thin*itsave){
modelres[[1]][,,itsave]<-scamRes$start$paramX
modelres[[2]][,,itsave]<-scamRes$start$paramTr
modelres[[3]][,,itsave]<-scamRes$start$means
modelres[[4]][,,itsave]<-scamRes$start$sigma
itsave<-itsave+1
}
}
### Print number of iterations completed
if(verbose){
if(any(i==printCounter)){
print(i)
}
}
}
}else{
for(i in 1:niter){
#---------------
### Fixed effect
#---------------
scamResHT<-scamHT(data,scamRes$start,priors,scamRes$rotation$values,scamRes$rotation$vectors,scamRes$rotation$wParamX,scamRes$rotation$wOuterParamX,scamRes$sdvalue$sdvalue,scamRes$accept$accept,likelihoodtype=likelihoodtype,burning=(i < nburn),rotation=(i > nrot & i < nburn),iter=i) # Pourrait changer
### Assign the new estimated pieces of scamResH to scamRes
scamRes$start$paramX<-scamResHT$start$paramX
scamRes$start$sigma<-scamResHT$start$sigma
scamRes$start$means<-scamResHT$start$means
scamRes$start$paramTr<-scamResHT$start$paramTr
scamRes$rotation$values<-scamResHT$rotation$values
scamRes$rotation$vectors<-scamResHT$rotation$vectors
scamRes$rotation$wParamX<-scamResHT$rotation$wParamX
scamRes$rotation$wOuterParamX<-scamResHT$rotation$wOuterParamX
scamRes$sdvalue$sdvalue<-scamResHT$sdvalue
scamRes$accept$accept<-scamResHT$accept
#----------------
### Random effect
#----------------
scamResRandom<-scamRandom(data,start=scamRes$start,priors=priors,rotationVal=scamRes$rotation$valuesRandom,rotationVect=scamRes$rotation$vectorsRandom,wParamRandom=scamRes$rotation$wParamRandom,wOuterParamRandom=scamRes$rotation$wOuterParamRandom,wRandomCov=scamRes$rotation$wRandomCov,sdvalue=scamRes$sdvalue,accept=scamRes$accept,likelihoodtype=likelihoodtype,burning=(i < nburn),rotation=(i > nrot & i < nburn),iter=i)
### Assign the new estimated pieces of scamResH to scamRes
scamRes$start$paramRandom<-scamResRandom$start$paramRandom
scamRes$start$RandomVar<-scamResRandom$start$RandomVar
scamRes$rotation$valuesRandom<-scamResRandom$rotation$valuesRandom
scamRes$rotation$vectorsRandom<-scamResRandom$rotation$vectorsRandom
scamRes$rotation$wParamRandom<-scamResRandom$rotation$wParamRandom
scamRes$rotation$wOuterParamRandom<-scamResRandom$rotation$wOuterParamRandom
scamRes$rotation$wRandomCov<-scamResRandom$rotation$wRandomCov
scamRes$sdvalue<-scamResRandom$sdvalue
scamRes$accept<-scamResRandom$accept
#-------------------------
#### Store burning results
#-------------------------
if(i < nburn){
burningObj$paramX[,,i]<-scamRes$start$paramX
burningObj$paramTr[,,i]<-scamRes$start$paramTr
burningObj$means[,,i]<-scamRes$start$means
burningObj$sigma[,,i]<-scamRes$start$sigma
burningObj$paramRandom[,,i]<-scamRes$start$paramRandom
burningObj$RandomVar[i]<-scamRes$start$RandomVar
if(details){
detailsObj$sdvalue[,,i+1]<-scamRes$sdvalue$sdvalue
detailsObj$acceptRatio[,,i+1]<-scamRes$accept$accept[,,1]/scamRes$accept$accept[,,2]
detailsObj$sdvalueRandomCov[i+1]<-scamRes$sdvalue$sdvalueRandomCov[1]
detailsObj$acceptRatioRandomCov[i+1]<-scamRes$accept$acceptRandomCov[1]/scamRes$accept$acceptRandomCov[2]
detailsObj$sdvalueRandom[,,i+1]<-scamRes$sdvalue$sdvalueRandom
detailsObj$acceptRatioRandom[,,i+1]<-scamRes$accept$acceptRandom[,,1]/scamRes$accept$acceptRandom[,,2]
}
}else{
it<-i-nburn
if(it==thin*itsave){
modelres$paramX[,,itsave]<-scamRes$start$paramX
modelres$paramTr[,,itsave]<-scamRes$start$paramTr
modelres$means[,,itsave]<-scamRes$start$means
modelres$sigma[,,itsave]<-scamRes$start$sigma
modelres$paramRandom[,,itsave]<-scamRes$start$paramRandom
modelres$RandomVar[itsave]<-scamRes$start$RandomVar
itsave<-itsave+1
}
}
### Print number of iterations completed
if(verbose){
if(any(i==printCounter)){
print(i)
}
}
}
}
if(!any(class(data)=="HMSCdataRandom")){
modelres$sdvalue<-scamRes$sdvalue
modelres$acceptRatio<-scamRes$accept[,,1]/scamRes$accept[,,2]
}else{
modelres$sdvalue<-scamRes$sdvalue$sdvalue
modelres$acceptRatio<-scamRes$accept$accept[,,1]/scamRes$accept$accept[,,2]
modelres$sdvalueRandomCov<-scamRes$sdvalue$sdvalueRandomCov
modelres$acceptRatioRandomCov<-scamRes$accept$acceptRandomCov[1]/scamRes$accept$acceptRandomCov[2]
modelres$sdvalueRandom<-scamRes$sdvalue$sdvalueRandom
modelres$acceptRatioRandom<-scamRes$accept$acceptRandom[,,1]/scamRes$accept$acceptRandom[,,2]
}
#===================
#### Name the pieces
#===================
### paramX
dimnames(modelres$paramX)[[3]]<-paste("iter",1:dim(modelres[[1]])[3],sep="")
dimnames(modelres$paramX)[[1]]<-rownames(start$paramX)
dimnames(modelres$paramX)[[2]]<-colnames(start$paramX)
### paramTr
dimnames(modelres$paramTr)[[3]]<-paste("iter",1:dim(modelres[[1]])[3],sep="")
dimnames(modelres$paramTr)[[1]]<-rownames(start$paramTr)
dimnames(modelres$paramTr)[[2]]<-colnames(start$paramTr)
### means
dimnames(modelres$means)[[3]]<-paste("iter",1:dim(modelres[[1]])[3],sep="")
dimnames(modelres$means)[[1]]<-colnames(start$paramX)
dimnames(modelres$means)[[2]]<-rownames(start$paramX)
### sigma
dimnames(modelres$sigma)[[3]]<-paste("iter",1:dim(modelres[[1]])[3],sep="")
dimnames(modelres$sigma)[[1]]<-colnames(start$paramX)
dimnames(modelres$sigma)[[2]]<-colnames(start$paramX)
### paramRandom
if(any(class(data)=="HMSCdataRandom")){
dimnames(modelres$paramRandom)[[1]]<-rownames(start$paramRandom)
dimnames(modelres$paramRandom)[[2]]<-colnames(start$paramRandom)
}
#==========
### Results
#==========
if(details){
res<-list(call=call,model=modelres,burning=burningObj,details=detailsObj,data=data)
}else{
res<-list(call=call,model=modelres,burning=burningObj,data=data)
}
attr(res,"likelihood")<-likelihoodtype
if(!any(class(data)=="HMSCdataRandom")){
class(res)<-c("hmsc","HMSCdata","HMSCdataTrait")
}else{
class(res)<-c("hmsc","HMSCdata","HMSCdataTrait","HMSCdataRandom")
}
return(res)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.