hmsc: Hierarchical modelling of species community

Description Usage Arguments Details Value Note Author(s) References See Also Examples

Description

Hierarchical modelling of species community. hmscH is designed to model a community with habitat characteristics only whereas hmscHT also uses species traits to construct the model.

Usage

1
2
3
4
5
6
7
## S3 method for class 'formula'
hmscH(formula, data, start, priors = NULL, niter = 12000, nburn = 2000, nrot = 100, thin=1, likelihoodtype = "logit", details = TRUE, verbose=TRUE, ...)
## S3 method for class 'HMSCdata'
hmscH(data, start, priors = NULL, niter = 12000, nburn = 2000, nrot = 100, thin=1, likelihoodtype = "logit",spAsso = FALSE, details = TRUE, verbose=TRUE, ...)

## S3 method for class 'HMSCdata'
hmscHT(data, start, priors = NULL, niter = 12000, nburn = 2000, nrot = 100, thin=1, likelihoodtype = "logit", details = TRUE, verbose=TRUE, ...)

Arguments

formula

Model formula. The left side presents the response matrix and the right side gives the descriptors with their parameters. (See details)

data

Object of class HMSCdata.

start

Object of class HMSCparam.

priors

Object of class HMSCprior.

niter

A positive integer defining the number of iterations to be carried out. Default is 12000.

nburn

Numeric. A positive integer defining the number of iterations to be used in the burning (first) phase of the algorithm. Default is 2000.

nrot

Numeric. A positive integer defining the number of iterations that needs to be carried out before eigenvector rotation is performed.

thin

Numeric. A positive integer defining thinning. Default is 1, that is all values are considered. (see details)

likelihoodtype

Character string defining how the likelihood should be calculated. Any unambiguous variation of wording used in this argument is accepted. (See details)

spAsso

Logical. Whether species association should be estimated (TRUE) or not (FALSE).

details

Logical. Whether detailed information about the burning, the accept ratio, and the standard deviation should be given. Default is TRUE.

verbose

Logical. Whether a counter printing on the screen how many iterations have been performed so far is used or not. Default is TRUE.

...

Parameters passed to other functions.

Details

For the user the main difference between hmscH and hmscHT is within the data and start arguments.

The formula was designed so that it is necessary to include the parameters and the descriptors. This function is thus capable of carrying out non-linear analyses (although the current implement is far from being very efficient).

When building a model for the first time, the parameters provided in start can be randomly generated using rnorm[stats], as it is shown in the examples. However, nothing prevents the user to predefine the starting values so that the function can more quickly converge to optimal set of parameters.

When thinning, the value given to thin means that all but the k^th values will be considered.

Currently, the following likelihood have been implemented:

- logit: Binomial likelihood using a logit link function

- poisson: Poisson likelihood using a log link function

Value

model

List with the following items:

paramX

A 3-dimensional array where, for each species (first dimension), the regression parameters of X (second dimension) are stored for all iterations (third dimension) after burning.

paramTr

A 3-dimensional array where, for each regression parameter (first dimension), the species traits (second dimension) are stored for all iterations (third dimension) after burning.

means

Matrix of means measuring the response of a typical response variables to the descriptors. Each row represent an iteration result after burning.

sigma

A 3-dimensional array where covariance matrices measuring how the species vary in their responses to the descriptors (first and second dimension) are stored for all iterations (third dimension) after burning.

sdvalue

Matrix of standard deviation used when sampling for the regression parameters (columns) of each response variables (row) after burning.

acceptRatio

Matrix of accept ratio after burning for each species (rows) when sampling for the regression parameters of X (columns). These values should be close to 0.44. If they are too far from 0.44, the model may not be fitted well. The number of burning iterations should be increased.

paramRandom

A 3-dimensional array where, for each species (first dimension), the regression parameters for each levels of the random effect (second dimension) are stored for all iterations (third dimension) after burning.

sdvalueRandom

Matrix of standard deviation used when sampling for the regression parameters for each levels of the random effect (columns) of each response variables (row) after burning.

acceptRatioRandom

Matrix of accept ratio after burning for each species (rows) when sampling for the regression parameters for each levels of the random effect (columns). These values should be close to 0.44. If they are too far from 0.44, the model may not be fitted well. The number of burning iterations should be increased.

RandomVar

Scalar defining the within group variance in Random after burning.

sdvalueRandomCov

Standard deviation used when sampling for RandomCov after burning.

acceptRatioRandomCov

accept ratio after burning for RandomCov after burning. This values should be close to 0.44. If it is too far from 0.44, the model may not be fitted well. The number of burning iterations should be increased.

burning

List with the following items :

paramX

A 3-dimensional array where, for each species (first dimension), the regression parameters (second dimension) are stored for all iterations (third dimension) during burning.

paramTr

A 3-dimensional array where, for each species traits (second dimension), the regression parameter (first dimension) are stored for all iterations (third dimension) during burning.

means

Matrix of means measuring the response of a typical response variables to the descriptors. Each row represent an iteration result during burning.

sigma

A 3-dimensional array where covariance matrices measuring how the species vary in their responses to the descriptors (first and second dimension) are stored for all iterations (third dimension) during burning.

paramRandom

A 3-dimensional array where, for each species (first dimension), the regression parameters for each levels of the random effect (second dimension) are stored for all iterations (third dimension) during burning.

RandomVar

A vector defining the within group variance in Random during burning.

details

This part of the result is only available when the argument details is TRUE List with the following items :

sdvalue

A 3-dimensional array where the first two dimensions represent the standard deviation used when sampling for the regression parameters (columns) of each response variables (row) for the different burning iterations (third dimension).

acceptRatio

A 3-dimensional array where the first two dimensions represent the accept ratio calculated during the burning phase of the model construction when sampling for the regression parameters (columns) of each response variables (row). The accept ratio for all burning iterations are presented as the third dimension.

sdvalueRandom

A 3-dimensional array where the first two dimensions represent the standard deviation used when sampling for the regression parameters of each levels of the random effect (columns) for each response variables (row) for the different burning iterations (third dimension).

acceptRatioRandom

A 3-dimensional array where the first two dimensions represent the accept ratio calculated during the burning phase of the model construction when sampling for the regression parameters of each levels of the random effect (columns) for each response variables (row). The accept ratio for all burning iterations are presented as the third dimension.

sdvalueRandomCov

A vector showing how the standard deviation varies when sampling for RandomCov for the different burning iterations.

acceptRatioRandomCov

A vector showing how the accept ratio varies when sampling for RandomCov for the different burning iterations.

data

An object of class HMSCdata that includes the data used to construct the model. Note that if traits and/or a random effect are considered in the model, data could also have a HMSCdataTrait and/or a HMSCdataRandom class.

Note

These functions are essentially glorified for loops using scamH or scamHT and/or scamRandom, repetitively.

Author(s)

F. Guillaume Blanchet

References

Ovaskainen, O. and Soininen, J. (2011) Making more out of sparse data: hierarchical modeling of species communities. Ecology 92, 289-295.

See Also

scamH

Examples

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
#==================
### Simulating data
#==================
desc<-cbind(1,scale(1:50),scale(1:50)^2)
nspecies<-20
dataBase<-communitySimulH(desc,nsp=nspecies)

#=================================
### Formatting data and parameters
#=================================
formdata<-as.HMSCdata(dataBase$data$Y,desc,Ypattern="sp",interceptX=FALSE,interceptTr=FALSE)

#========================
### Formatting parameters
#========================
startParamX<-matrix(rnorm(nspecies*ncol(desc)),nrow=nspecies,ncol=ncol(desc))
startMean<-rnorm(ncol(desc))
startSigma<-rWishart(1,ncol(desc)+1,diag(ncol(desc)))[,,1]

formparam<-as.HMSCparam(formdata,paramX=startParamX,means=startMean,sigma=startSigma)
formpriors<-as.HMSCprior(formparam,rep(0,length(formparam$means)),0,length(formparam$means),diag(length(formparam$means)))

#=============================================================================================================
### Building model using a linear modelling approach (currently 2000 iterations takes about 50 seconds to run)
#=============================================================================================================
system.time(
modelLinear<-hmscH(formdata,formparam,formpriors,niter=200,nburn=100)
)

#=============================================================================================================
### Building model using a formula modelling approach (currently 2000 iterations takes about 2 minutes to run)
#=============================================================================================================
formule<-sp~p1*x1+p2*x2+p3*x3
system.time(
modelFormula<-hmscH(formule,formdata,formparam,formpriors,niter=200,nburn=100)
)

#===============
### Some results
#===============
### Linear model
summaryLin<-summary(modelLinear)
summaryLin[1:5,1:5,]

### Formula model
summaryForm<-summary(modelFormula)
summaryForm[1:5,1:5,]

#------------------------------
### Draw model from modelLinear
#------------------------------

par(mfrow=c(5,4),mar=c(1,3,2,1))
for(i in 1:20){
	plot(summaryLin[,i,1],type="l",ylim=c(0,1),las=1,lwd=3,xaxt="n",col="red")
	title(paste("Species",i))
	lines(summaryLin[,i,3],col="red",lwd=1)
	lines(summaryLin[,i,5],col="red",lwd=1)
	lines(dataBase$probMat[,i],col="black",lwd=3)
}

#-------------------------------
### Draw model from modelFormula
#-------------------------------
dev.new()
par(mfrow=c(5,4),mar=c(1,3,2,1))
for(i in 1:20){
	plot(summaryForm[,i,1],type="l",ylim=c(0,1),las=1,lwd=3,xaxt="n",col="red")
	title(paste("Species",i))
	lines(summaryForm[,i,3],col="red",lwd=1)
	lines(summaryForm[,i,5],col="red",lwd=1)
	lines(dataBase$probMat[,i],col="black",lwd=3)
}

#------------------
### Mixing - Linear
#------------------
#___________________
#### Community Sigma
#___________________
par(mfrow=c(3,3),mar=c(1,4,3,1))
for(i in 1:3){
	for(j in 1:3){
		mixing(modelLinear,toplot="sigma",param=c(i,j),col="red",xlab="Sites",ylab="Sigma",ylim=range(c(modelLinear$model$sigma[i,j,],dataBase$param$sigma[i,j])),las=1)
		abline(h=dataBase$param$sigma[i,j],col="red",lwd=3)
		title(paste("Community Sigma",i,"-",j))
	}
}

#__________________
### Community means
#__________________
par(mfrow=c(1,1),mar=c(5,4,4,2))
mixing(modelLinear,toplot="means",param=1,ylim=range(rbind(modelLinear$model$means,dataBase$param$means)),xlab="Sites",ylab="Means",las=1)
title("Community means")
mixing(modelLinear,toplot="means",param=2,add=TRUE,col="blue",xlab="Sites",ylab="Means",las=1)
mixing(modelLinear,toplot="means",param=3,add=TRUE,col="red",xlab="Sites",ylab="Means",las=1)
abline(h=dataBase$param$means,col=c("black","blue","red"),lwd=3)

#_____________________________
### Parameter for each species 
#_____________________________
par(mfrow=c(5,4),mar=c(1,4,2,1))
for(i in 1:20){
	mixing(modelLinear,toplot="paramX",param=c(i,1),ylim=range(cbind(modelLinear$model$paramX[i,,],dataBase$param$paramX[i,])),xlab="Sites",ylab="paramX",las=1,xaxt="n")
	title(paste("Species",i))
	mixing(modelLinear,toplot="paramX",param=c(i,2),add=TRUE,col="blue")
	mixing(modelLinear,toplot="paramX",param=c(i,3),add=TRUE,col="red")
	abline(h=dataBase$param$paramX[i,],col=c("black","blue","red"),lwd=3)
}

#------------------------------------------
### Check the distribution of the posterior
#------------------------------------------
#__________________
### Community means
#__________________
pairs(modelLinear$model$means,asp=1,pch=19,cex=0.2)

#-------------------
### Mixing - Formula
#-------------------
#___________________
#### Community Sigma
#___________________
par(mfrow=c(3,3),mar=c(1,3,3,1))
for(i in 1:3){
	for(j in 1:3){
		mixing(modelFormula,toplot="sigma",param=c(i,j),col="red",xlab="Sites",ylab="Sigma",ylim=range(c(modelFormula$model$sigma[i,j,],dataBase$param$sigma[i,j])),las=1)
		abline(h=dataBase$param$sigma[i,j],col="red",lwd=3)
		title(paste("Community Sigma",i,"-",j))
	}
}

#__________________
### Community means
#__________________
par(mfrow=c(1,1),mar=c(5,4,4,2))
mixing(modelFormula,toplot="means",param=1,ylim=range(rbind(modelFormula$model$means,dataBase$param$means)),xlab="Sites",ylab="Means",las=1)
title("Community means")
mixing(modelFormula,toplot="means",param=2,add=TRUE,col="blue",xlab="Sites",ylab="Means",las=1)
mixing(modelFormula,toplot="means",param=3,add=TRUE,col="red",xlab="Sites",ylab="Means",las=1)
abline(h=dataBase$param$means,col=c("black","blue","red"),lwd=3)

#_____________________________
### Parameter for each species 
#_____________________________
par(mfrow=c(5,4),mar=c(1,3,2,1))
for(i in 1:20){
	mixing(modelFormula,toplot="paramX",param=c(i,1),ylim=range(cbind(modelFormula$model$paramX[i,,],dataBase$param$param[i,])),xlab="Sites",ylab="param",las=1,xaxt="n")
	title(paste("Species",i))
	mixing(modelFormula,toplot="paramX",param=c(i,2),add=TRUE,col="blue")
	mixing(modelFormula,toplot="paramX",param=c(i,3),add=TRUE,col="red")
	abline(h=dataBase$param$param[i,],col=c("black","blue","red"),lwd=3)
}

#------------------------------------------
### Check the distribution of the posterior
#------------------------------------------
#__________________
### Community means
#__________________
pairs(modelFormula$model$means,asp=1,pch=19,cex=0.2)

########################
### Example using Traits
########################
#==============================
### Simulating data with traits
#==============================
desc <- cbind(1,scale(1:50),scale(1:50)^2)
trait <- rbind(1,1:20/5)
dataTrait<-communitySimulHT(desc,trait)

#=================================
### Formatting data and parameters
#=================================
formdata<-as.HMSCdata(dataTrait$data$Y,desc,trait,Ypattern="sp",interceptX=FALSE,interceptTr=FALSE)

#========================
### Formatting parameters
#========================
startParamX<-matrix(rnorm(nrow(trait)*ncol(desc)),nrow=ncol(dataTrait$data$Y),ncol=ncol(desc))
startParamTr<-matrix(rnorm(ncol(trait)*ncol(desc)),nrow=ncol(desc),ncol=nrow(trait))
startSigma<-rWishart(1,ncol(desc)+1,diag(ncol(desc)))[,,1]

formparam<-as.HMSCparam(formdata,paramX=startParamX,paramTr=startParamTr,sigma=startSigma)
formpriors<-as.HMSCprior(formparam,kappa0=0,nu0=length(formparam$means),Lambda0=diag(nrow(formparam$sigma)))

#=============================================================================================================
### Building model using a linear modelling approach (currently 2000 iterations takes about 50 seconds to run)
#=============================================================================================================
system.time(
modelTraits<-hmscHT(formdata,formparam,formpriors,niter=200,nburn=100)
)

### Traits model
summaryTraits<-summary(modelTraits)
summaryTraits[1:5,1:5,]

#------------------
### Mixing - Traits
#------------------
#___________________
#### Community Sigma
#___________________
par(mfrow=c(3,3),mar=c(1,3,3,1))
for(i in 1:3){
	for(j in 1:3){
		mixing(modelTraits,toplot="sigma",param=c(i,j),col="red",xlab="Sites",ylab="Sigma",ylim=range(c(modelTraits$model$sigma[i,j,],dataTrait$param$sigma[i,j])),las=1)
		abline(h=dataTrait$param$sigma[i,j],col="red",lwd=3)
		title(paste("Community Sigma",i,"-",j))
	}
}

#__________________
### Community means
#__________________
par(mfrow=c(5,4),mar=c(1,4,2,1))
for(i in 1:20){
	mixing(modelTraits,toplot="means",param=c(1,i),ylim=range(cbind(modelTraits$model$means[,i,],dataTrait$param$means[,i])),xlab="Sites",ylab="means",las=1,xaxt="n")
	title(paste("Species",i))
	mixing(modelTraits,toplot="means",param=c(2,i),add=TRUE,col="blue")
	mixing(modelTraits,toplot="means",param=c(3,i),add=TRUE,col="red")
	abline(h=dataTrait$param$means[,i],col=c("black","blue","red"),lwd=3)
}

#_____________________________
### Parameter for each species 
#_____________________________
par(mfrow=c(5,4),mar=c(1,4,2,1))
for(i in 1:20){
	mixing(modelTraits,toplot="paramX",param=c(i,1),ylim=range(cbind(modelTraits$model$paramX[i,,],dataTrait$param$paramX[i,])),xlab="Sites",ylab="paramX",las=1,xaxt="n")
	title(paste("Species",i))
	mixing(modelTraits,toplot="paramX",param=c(i,2),add=TRUE,col="blue")
	mixing(modelTraits,toplot="paramX",param=c(i,3),add=TRUE,col="red")
	abline(h=dataTrait$param$paramX[i,],col=c("black","blue","red"),lwd=3)
}

#____________________________
### Parameter for each traits 
#____________________________
par(mfrow=c(3,1),mar=c(1,4,2,1))
for(i in 1:3){
	mixing(modelTraits,toplot="paramTr",param=c(i,1),ylim=range(cbind(modelTraits$model$paramTr[i,,],dataTrait$param$paramTr[i,])),xlab="Sites",ylab="paramTr",las=1,xaxt="n")
	title(paste("Habitat",i))
	mixing(modelTraits,toplot="paramTr",param=c(i,2),add=TRUE,col="blue")
	abline(h=dataTrait$param$paramTr[i,],col=c("black","blue"),lwd=3)
}

#------------------------------
### Draw model from modelTraits
#------------------------------
par(mfrow=c(5,4),mar=c(1,3,2,1))
for(i in 1:20){
	plot(summaryTraits[,i,1],type="l",ylim=c(0,1),las=1,lwd=3,xaxt="n",col="red")
	title(paste("Species",i))
	lines(summaryTraits[,i,3],col="red",lwd=1)
	lines(summaryTraits[,i,5],col="red",lwd=1)
	lines(dataTrait$probMat[,i],col="black",lwd=3)
}

#################################
### Example using a Random effect
#################################
#==============================
### Simulating data with traits
#==============================
desc <- cbind(1,scale(1:50),scale(1:50)^2)
arbi<-as.factor(c(1,rep(1:3,each=16),3))
dataRandom<-communitySimulH(desc,Random=arbi,nsp=20)

#=================================
### Formatting data and parameters
#=================================
formdata<-as.HMSCdata(dataRandom$data$Y,desc,Random=arbi,Ypattern="sp",interceptX=FALSE,scaleX=FALSE)

#========================
### Formatting parameters
#========================
startParamX<-matrix(rnorm(ncol(dataRandom$data$Y)*ncol(desc)),nrow=ncol(dataRandom$data$Y),ncol=ncol(desc))
startSigma<-rWishart(1,ncol(desc)+1,diag(ncol(desc)))[,,1]
startRandom<-matrix(rnorm(ncol(dataRandom$data$Y)*nlevels(arbi)),nrow=ncol(dataRandom$data$Y),ncol=nlevels(arbi))
startMean<-rnorm(ncol(desc))

formparam<-as.HMSCparam(formdata,paramX=startParamX,paramRandom=startRandom,means=startMean,sigma=startSigma,RandomVar=dataRandom$param$RandomCommSp+1,RandomCommSp=dataRandom$param$RandomCommSp)

#============================================================================
### Building model (currently 2000 iterations takes about 2.5 minutes to run)
#============================================================================
system.time(
modelRandom<-hmscH(formdata,formparam,niter=200,nburn=100)
)

### Traits model
summaryRandom<-summary(modelRandom)
summaryRandom[1:5,1:5,]

#------------------
### Mixing - Random
#------------------
#___________________
#### Community Sigma
#___________________
par(mfrow=c(3,3),mar=c(1,3,3,1))
for(i in 1:3){
	for(j in 1:3){
		mixing(modelRandom,toplot="sigma",param=c(i,j),col="red",xlab="Sites",ylab="Sigma",ylim=range(c(modelRandom$model$sigma[i,j,],dataRandom$param$sigma[i,j])),las=1)
		abline(h=dataRandom$param$sigma[i,j],col="red",lwd=3)
		title(paste("Community Sigma",i,"-",j))
	}
}

#__________________
### Community means
#__________________
par(mfrow=c(3,1),mar=c(1,4,2,1))
for(i in 1:3){
	mixing(modelRandom,toplot="means",param=i,col="red",ylim=range(cbind(modelRandom$model$means[,i],dataRandom$param$means[i])),xlab="Sites",ylab="means",las=1,xaxt="n")
	abline(h=dataRandom$param$means[i],col="red",lwd=3)
}

#__________________________________
### Parameter of X for each species 
#__________________________________
par(mfrow=c(5,4),mar=c(1,4,2,1))
for(i in 1:20){
	mixing(modelRandom,toplot="paramX",param=c(i,1),ylim=range(cbind(modelRandom$model$paramX[i,,],dataRandom$param$paramX[i,])),xlab="Sites",ylab="paramX",las=1,xaxt="n")
	title(paste("Species",i))
	mixing(modelRandom,toplot="paramX",param=c(i,2),add=TRUE,col="blue")
	mixing(modelRandom,toplot="paramX",param=c(i,3),add=TRUE,col="red")
	abline(h=dataRandom$param$paramX[i,],col=c("black","blue","red"),lwd=3)
}

#_________________________________________________
### Parameter for each levels of the Random effect
#_________________________________________________
par(mfrow=c(5,4),mar=c(1,4,2,1))
for(i in 1:20){
	mixing(modelRandom,toplot="paramRandom",param=c(i,1),ylim=range(cbind(modelRandom$model$paramRandom[i,,],dataRandom$param$paramRandom[i,])),xlab="Sites",ylab="paramRandom",las=1,xaxt="n")
	title(paste("Species",i))
	mixing(modelRandom,toplot="paramRandom",param=c(i,2),add=TRUE,col="blue")
	mixing(modelRandom,toplot="paramRandom",param=c(i,3),add=TRUE,col="red")
	abline(h=dataRandom$param$paramRandom[i,],col=c("black","blue","red"),lwd=3)
}

#________________________________
### Variance of the Random effect
#________________________________
par(mfrow=c(1,1),mar=c(1,4,2,1))
mixing(modelRandom,toplot="RandomVar",param=i,col="red",ylim=range(c(modelRandom$model$RandomVar,dataRandom$param$RandomVar)),xlab="Sites",ylab="RandomVar",las=1,xaxt="n")
abline(h=dataRandom$param$RandomVar,col="red",lwd=3)

#------------------------------
### Draw model from modelRandom
#------------------------------
par(mfrow=c(5,4),mar=c(1,3,2,1))
for(i in 1:20){
	plot(summaryRandom[,i,1],type="l",ylim=c(0,1),las=1,lwd=3,xaxt="n",col="red")
	title(paste("Species",i))
	lines(summaryRandom[,i,3],col="red",lwd=1)
	lines(summaryRandom[,i,5],col="red",lwd=1)
	lines(dataRandom$probMat[,i],col="black",lwd=3)
}

##########################################
### Example using Traits and Random effect
##########################################
#==============================
### Simulating data with traits
#==============================
desc <- cbind(1,scale(1:50),scale(1:50)^2)
arbi<-as.factor(c(1,rep(1:3,each=16),3))
trait <- rbind(1,1:20/5)
dataTraitsRandom<-communitySimulHT(desc,trait,Random=arbi)

#=================================
### Formatting data and parameters
#=================================
formdata<-as.HMSCdata(dataTraitsRandom$data$Y,desc,trait,Random=arbi,Ypattern="sp",interceptX=FALSE,interceptTr=FALSE)

#========================
### Formatting parameters
#========================
startParamX<-matrix(rnorm(nrow(trait)*ncol(desc)),nrow=ncol(dataTraitsRandom$data$Y),ncol=ncol(desc))
startParamTr<-matrix(rnorm(ncol(trait)*ncol(desc)),nrow=ncol(desc),ncol=nrow(trait))
startSigma<-rWishart(1,ncol(desc)+1,diag(ncol(desc)))[,,1]
startRandom<-matrix(rnorm(ncol(dataTraitsRandom$data$Y)*nlevels(arbi)),nrow=ncol(dataTraitsRandom$data$Y),ncol=nlevels(arbi))

formparam<-as.HMSCparam(formdata,paramX=startParamX,paramRandom=startRandom,paramTr=startParamTr,sigma=startSigma,RandomVar=dataTraitsRandom$param$RandomVar,RandomCommSp=dataTraitsRandom$param$RandomCommSp)

#=============================================================================================================
### Building model using a linear modelling approach (currently 2000 iterations takes about 50 seconds to run)
#=============================================================================================================
system.time(
modelTraitsRandom<-hmscHT(formdata,formparam,niter=200,nburn=100)
)

### Traits model
summaryTraitsRandom<-summary(modelTraitsRandom)
summaryTraitsRandom[1:5,1:5,]

#-------------------------
### Mixing - Traits Random
#-------------------------
#___________________
#### Community Sigma
#___________________
par(mfrow=c(3,3),mar=c(1,3,3,1))
for(i in 1:3){
	for(j in 1:3){
		mixing(modelTraitsRandom,toplot="sigma",param=c(i,j),col="red",xlab="Sites",ylab="Sigma",ylim=range(c(modelTraitsRandom$model$sigma[i,j,],dataTraitsRandom$param$sigma[i,j])),las=1)
		abline(h=dataTraitsRandom$param$sigma[i,j],col="red",lwd=3)
		title(paste("Community Sigma",i,"-",j))
	}
}

#__________________
### Community means
#__________________
par(mfrow=c(5,4),mar=c(1,4,2,1))
for(i in 1:20){
	mixing(modelTraitsRandom,toplot="means",param=c(1,i),ylim=range(cbind(modelTraitsRandom$model$means[,i,],dataTraitsRandom$param$means[,i])),xlab="Sites",ylab="means",las=1,xaxt="n")
	title(paste("Species",i))
	mixing(modelTraitsRandom,toplot="means",param=c(2,i),add=TRUE,col="blue")
	mixing(modelTraitsRandom,toplot="means",param=c(3,i),add=TRUE,col="red")
	abline(h=dataTraitsRandom$param$means[,i],col=c("black","blue","red"),lwd=3)
}

#__________________________________
### Parameter of X for each species 
#__________________________________
par(mfrow=c(5,4),mar=c(1,4,2,1))
for(i in 1:20){
	mixing(modelTraitsRandom,toplot="paramX",param=c(i,1),ylim=range(cbind(modelTraitsRandom$model$paramX[i,,],dataTraitsRandom$param$paramX[i,])),xlab="Sites",ylab="paramX",las=1,xaxt="n")
	title(paste("Species",i))
	mixing(modelTraitsRandom,toplot="paramX",param=c(i,2),add=TRUE,col="blue")
	mixing(modelTraitsRandom,toplot="paramX",param=c(i,3),add=TRUE,col="red")
	abline(h=dataTraitsRandom$param$paramX[i,],col=c("black","blue","red"),lwd=3)
}

#____________________________
### Parameter for each traits 
#____________________________
par(mfrow=c(3,1),mar=c(1,4,2,1))
for(i in 1:3){
	mixing(modelTraitsRandom,toplot="paramTr",param=c(i,1),ylim=range(cbind(modelTraitsRandom$model$paramTr[i,,],dataTraitsRandom$param$paramTr[i,])),xlab="Sites",ylab="paramTr",las=1,xaxt="n")
	title(paste("Habitat",i))
	mixing(modelTraitsRandom,toplot="paramTr",param=c(i,2),add=TRUE,col="blue")
	abline(h=dataTraitsRandom$param$paramTr[i,],col=c("black","blue"),lwd=3)
}

#------------------------------------
### Draw model from modelTraitsRandom
#------------------------------------
par(mfrow=c(5,4),mar=c(1,3,2,1))
for(i in 1:20){
	plot(summaryTraitsRandom[,i,1],type="l",ylim=c(0,1),las=1,lwd=3,xaxt="n",col="red")
	title(paste("Species",i))
	lines(summaryTraitsRandom[,i,3],col="red",lwd=1)
	lines(summaryTraitsRandom[,i,5],col="red",lwd=1)
	lines(dataTraitsRandom$probMat[,i],col="black",lwd=3)
}

HMSC documentation built on May 31, 2017, 5:14 a.m.