/usr/lib/R/site-library/TH.data/PSGLMM_MEE/README.txt is in r-cran-th.data 1.0-8-1build1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
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 | **********************************************************************
*
* Reproducibility Information for
*
* Spatio-Phylogenetic Multi-Species Distribution Models
* by A. Kaldhusdal, R. Brandl, J. Mueller, L. Moest and T. Hothorn
*
**********************************************************************
Content
### bird data
data/birdsData.Rdata
### analysis of bird data
analysis_birdData/
mBirdsESM1_2.R mBirdsESM3_2.R
### PSGLMM R code
R/
makeData.R makeEnvir.R psglmm.R psglmmSim.R
### simulation study
simulationStudy/
simModel11bin.R simModel12.R simModel21pois.R simModel23bin.R
simModel11pois.R simModel13bin.R simModel21.R simModel23pois.R
simModel11.R simModel13pois.R simModel22bin.R simModel23.R
simModel12bin.R simModel13.R simModel22pois.R
simModel12pois.R simModel21bin.R simModel22.R
Description
A PSGLMM may be fitted by calling the main function psglmm(). This function
utilises the functions lmer() (for Gaussian responses) or glmer() (for binary or
Poisson responses) from the lme4 package, but with a twist. Instead of having
the covariance structure in the data estimated by the function, psglmm()
requires the user to provide an assumption for this in advance, which is then
incorporated into the estimation procedure. The covariance structure is passed
to the function as a list containing one element per random effect component in
the model.
A formula is passed to the function in the known mixed model (lme4) form, e.g.
y ~ spp - 1 + (spp - 1 | plot).
With spp being a factor expressing to which species an observation belongs and
plot expressing to which plot the observation was made, the above model assumes
a fixed effect for each species and a random effect in the plots for each
species. At this stage psglmm() is merely an extension of (g)lmer() for which
the covariance structure may be assumed to be known.
If several possible covariance structures are plausible, psglmm() also allows us
to test which of these is best fitted to the data. In doing this we pass one
random effect component for each covariance structure to be tested, e.g.:
y ~ spp - 1 + (spp - 1 | plot) + (spp - 1 | plot).
With A_ind being a unity matrix with a size equal to the number of species,
A_phylo being a phylogenetic correlation matrix for the species and B being
a correlation matrix expressing the spatial dependency of the plots, we would
pass the following list of covariance matrices (VCV) to the function (naming
the list elements eases the interpretation of the model output):
list("independence"= kronecker(B, A_ind), "phylogeny"= kronecker(B, A_phylo)).
To perform the model selection, the argument msel must be utilised. The model
selection procedure performs a predictive cross-validation and offers three
different proper scoring rules for this task: LS or DSS for Gaussian responses,
BS or LS for binary responses and DSS for Poisson responses.
psglmm() first calls (g)lFormula() (lme4) to set up the needed model objects.
These are then modified using several helper functions before they are passed to
the lme4 fitting algorithms. These helper functions are briefly described in the
following:
- As lmer() is only allowed to fit less random effects than there are
observations and the above mentioned formulas assume one random effect for
each observation, we apply the code in lines 301-350 (see psglmm.R) to disable
this limitation. Without executing this code a Gaussian model of the above
form cannot be fitted.
- adaptFormula() is used to ensure that the order of the random effect
components in the formula passed on to the lme4 fitting algorithms is indeed
the same as in the formula passed to psglmm() by the user.
- randomTerms() is used to set up the random effect terms using the provided
covariance structures.
- makeVARALL() simply sets up a data.frame for the model results.
- modelSelection() is the function which is called if an automated model
selection is wanted (argument msel set to TRUE in psglmm()). The function in
turn calls the helper functions pred.xval() in which the predictive cross-
validation is performed. It also allows for an LR test to be performed to test
if the variance estimates for the covariance components are > 0 (because of
the issues concerning the testing of variance parameters on the boundary of
their parameter space, the calculated p-value should not be overly
interpreted).
- LR.test() performs the mentioned LR test.
- pred.xval() is the function in which the predictive cross-validation is
actually performed. The function is based on the algorithm outlined by Braun
et. al (2013). The iterated weighted least squares algorithm is hard coded to
abort if no convergence is achieved within 10000 iterations. If this happens,
or if any other numerical issue arises, the model selection is aborted, but
the model fitting is continued. A corresponding warning is displayed in the
model output. Depending on which proper scoring rule is chosen by the user
(argument msel passed to psglmm()), one of the following helper functions is
called:
- BS(), DSS(), LS() incorporate the three proper scoring rules "Brier Score",
"Dawid-Sebastiani Score" and "Logarithmic Score".
The following helper functions are used to simplify the calculation of some
tasks, which were performed often and thus save some space and make the code a
little easier to read:
- nuniq() returns the number of unique elements in a vector.
- equal() checks if all elements of a vector are identical.
- IFELSE() works like ifelse() but allows for the objects to be returned in
either case to be of different classes.
- makePolys() actually didn't find any use in the final version of the paper,
but is practical for setting up the polygons of the plots when displaying
the estimated values from the model.
The results can be reproduced using the following package versions:
lme4 version 1.0-0 (from R-forge), 1.1-6 from CRAN
Matrix version 1.1-4
fields version 6.8
mvtnorm version 0.9-9995
picante version 1.6-2
ggplot2 version 1.0.0
|