Poisson lognormal regression model in R

 library(PLNmodels)

library(ggplot2)

library(corrplot)

Y=read.csv(choose.files())

X=read.csv(choose.files())

c=list(Y,X)

names(c) <- c("output", "input")

d <- prepare_data(c$output, c$input)

model <- PLN(Abundance ~ 

             LANE_WIDTH               +

             LENGTH                    +

             LOG_AVG_AADT                 +

             LOG_PAVEMENT_CONDITION         +

             HIGH_FREQ_TRANSIT           +

             eightyfivepercent_SPEED_ADJ   +

             P_SPEEDING_ADJ               +

             SL                         +

             NL, data = d)

print(model)      #variational lower bound of the ICL

coef(model)       #mu -vectors of means of the latent variable

sigma(model)      #covariance matrix of the latent variable

vcov(model)       #Variance-Covariance Matrix

fitted(model)

standard_error(model)

plot(fitted(model))

barplot(fitted(model))

#2 tailed z test

z<- coef(model)/standard_error(model)

p<- (1-pnorm(abs(z),0,1))*2

p        #pvalue