library(PLNmodels)library(ggplot2)library(corrplot)

data(trichoptera)trichoptera<-prepare_data(trichoptera$Abundance, trichoptera$Covariate)

PLN: geometrical view

Covariates and offsets

This model generalizes naturally to a formulation closer to amultivariate generalized linear model, where the main effect is due to alinear combination of\(d\) covariates\(\mathbf{x}_i\) (including a vector ofintercepts). We also let the possibility to add some offsets for the\(p\) variables in in each sample, thatis\(\mathbf{o}_i\). Hence, theprevious model generalizes to

\[\begin{equation} \mathbf{Y}_i | \mathbf{Z}_i \sim\mathcal{P}\left(\exp\{\mathbf{Z}_i\}\right), \qquad \mathbf{Z}_i \sim\mathcal{N}({\mathbf{o}_i +\mathbf{x}_i^\top\mathbf{B}},\boldsymbol\Sigma), \\\end{equation}\] where\(\mathbf{B}\) is a\(d\times p\) matrix of regressionparameters. When all individuals\(i=1,\dots,n\) are stacked together, thedata matrices available to feed the model are

• the\(n\times p\) matrix of counts\(\mathbf{Y}\)
• the\(n\times d\) matrix of design\(\mathbf{X}\)
• the\(n\times p\) matrix of offsets\(\mathbf{O}\)

Inference in PLN then focuses on the regression parameters\(\mathbf{B}\) and on the covariance matrix\(\boldsymbol\Sigma\).

Optimization by Variational inference

Technically speaking, we adopt inPLNmodels avariational strategy to approximate the log-likelihood function andoptimize the consecutive variational surrogate of the log-likelihoodwith a gradient-ascent-based approach. To this end, we rely on the CCSAalgorithm ofSvanberg (2002) implemented in the C++ library(Johnson 2011), whichwe link to the package.

Adjusting a fit

In order to become familiar with the functionPLN andits outputs, let us first fit a simple PLN model with just an interceptfor each species:

myPLN<-PLN(Abundance~1, trichoptera)

## ##  Initialization...##  Adjusting a full covariance PLN model with nlopt optimizer##  Post-treatments...##  DONE!

Note the use of theformula object to specify the model:the vector\(\boldsymbol\mu\) of maineffects in the mathematical formulation (one per column species) isspecified in the call with the term~ 1 in theright-hand-side of the formula.Abundance is a variable inthe data frametrichoptera corresponding to a matrix of 17columns and theresponse in the model, occurring on theleft-hand-side of the formula.

ThePLNfit object

myPLN is anR6 object with classPLNfit, which comes with a couple of methods, as recalledwhen printing/showing such an object in theR console:

myPLN

## A multivariate Poisson Lognormal fit with full covariance model.## ==================================================================##  nb_param    loglik       BIC       ICL##       170 -1130.047 -1460.852 -2229.083## ==================================================================## * Useful fields##     $model_par, $latent, $latent_pos, $var_par, $optim_par##     $loglik, $BIC, $ICL, $loglik_vec, $nb_param, $criteria## * Useful S3 methods##     print(), coef(), sigma(), vcov(), fitted()##     predict(), predict_cond(), standard_error()

See also?PLNfit for more comprehensive information.

Field access

Accessing public fields of aPLNfit object can be donejust like with a traditional list,e.g.,

c(myPLN$loglik, myPLN$BIC, myPLN$ICL)

## [1] -1130.047 -1460.852 -2229.083

myPLN$criteria

##   nb_param    loglik       BIC       ICL## 1      170 -1130.047 -1460.852 -2229.083

GLM-like interface

We provide a set of S3-methods forPLNfit that mimic thestandard (G)LM-like interface ofR::stats, which we presentnow.

One can access the fitted value of the counts (Abundance–\(\hat{\mathbf{Y}}\)) and check thatthe algorithm basically learned correctly from the data¹:

data.frame(fitted   =as.vector(fitted(myPLN)),observed =as.vector(trichoptera$Abundance))%>%ggplot(aes(x = observed,y = fitted))+geom_point(size = .5,alpha =.25 )+scale_x_log10()+scale_y_log10()+theme_bw()+annotation_logticks()

fitted value vs. observation

The residual correlation matrix better displays as an imagematrix:

myPLN%>%sigma()%>%cov2cor()%>%corrplot()

Observation weights

It is also possible to use observation weights like in standard(G)LMs:

myPLN_weighted<-PLN(    Abundance~1,data    = trichoptera,weights =runif(nrow(trichoptera)),control =PLN_param(trace =0)  )data.frame(unweighted =as.vector(fitted(myPLN)),weighted   =as.vector(fitted(myPLN_weighted)))%>%ggplot(aes(x = unweighted,y = weighted))+geom_point(size = .5,alpha =.25 )+scale_x_log10()+scale_y_log10()+theme_bw()+annotation_logticks()

myPLN_offsets<-PLN(Abundance~1+offset(log(Offset)),data = trichoptera,control =PLN_param(trace =0))

rbind(  myPLN$criteria,  myPLN_offsets$criteria)%>% knitr::kable()

myPLN_wind<-PLN(Abundance~1+ Wind+offset(log(Offset)),data = trichoptera)

rbind(  myPLN_offsets$criteria,  myPLN_wind$criteria)%>% knitr::kable()

myPLN_spherical<-PLN(    Abundance~1+offset(log(Offset)),data = trichoptera,control =PLN_param(covariance ="spherical",trace =0)  )

myPLN_diagonal<-PLN(    Abundance~1+offset(log(Offset)),data = trichoptera,control =PLN_param(covariance ="diagonal",trace =0)  )

myPLN_default<-PLN(Abundance~1,data = trichoptera, )

myPLN_full<-PLN(Abundance~1,data = trichoptera,control =PLN_param(covariance ="full"))

rbind(  myPLN_offsets$criteria,  myPLN_diagonal$criteria,  myPLN_spherical$criteria)%>%as.data.frame(row.names =c("full","diagonal","spherical"))%>%  knitr::kable()

myPLN_final<-PLN(    Abundance~1+ Wind+offset(log(Offset)),data    = trichoptera,control =PLN_param(covariance ="diagonal",trace =0)  )rbind(  myPLN_wind$criteria,  myPLN_diagonal$criteria,  myPLN_final$criteria)%>% knitr::kable()