The results and specifics from the 'qs' function in R I'm not entirely familiar with the results of the qs() function in R's package seasonal. But, by what I understand roughly, the results below show me that the data does not show any seasonality in the first place -- if I am wrong, please correct me. But could someone give me an idea of how the variable 'qsorievadj' is constructed and how these tests are performed in R within the seasonal package? 
> qs(var)
                 qs   p-val
qsori       0.27402 0.87196
qsorievadj  0.00000 1.00000
qsrsd       0.00000 1.00000
qssadj      0.27402 0.87196
qssadjevadj 0.00000 1.00000
qsirr       0.00000 1.00000
qsirrevadj  0.00000 1.00000
qssrsd      0.00000 1.00000

 A: These tests are not actually computed in R. The seasonal package is 
an interface to X-13-ARIMA-SEATS, which is where the computations are performed.
The QS statistics are based on a draft research report by Maravall (2012, Research Report, Banco de España). Details are also available in the documentation of X-13-ARIMA-SEATS.

These statistics are relatively straightforward to implement in R. Below I give you the general sketch to compute the statistic. (I give an example of the QS statistic computed for the original series).
# run the whole procedure for the "AirPassengers" series
require(seasonal)
m <- seas(AirPassengers)

The steps shown below follow the description given in the documentation of X-13-ARIMA-SEATS. First, the series for which the QS statistic is computed is differenced according to the chosen ARIMA model and the following rule:
$$
ndif = \max(1, \min(d + D, 2))
$$
where $ndif$ is the number of regular differences to be taken; $d$ and $D$ are respectively the number of regular and seasonal differences in the chosen ARIMA model. (If the QS statistic is computed for the series of residuals no differences are applied, $ndif=0$.)
require(polynom)
x <- AirPassengers
S <- frequency(x)
ndif <- max(1, min(sum(arimamodel(m)[c(2,5)]), 2))
dx <- filter(x, polynomial(c(1,-1))^ndif, sides=1)
dx <- window(dx, start=time(x)[ndif+1])
# alternatively, we can do it without package "polynom" by simply doing,
# if ndiff=1:
#dx <- diff(x)
# if ndiff=2:
#dx <- diff(diff(x)) 

Next, the first two autocorrelations of seasonal order (e.g. 12 and 24 in monthly data) are obtained. If these autocorrelations are lower or equal to zero, then they are set to zero.
R <- acf(dx, lag.max=S*2, plot=FALSE)$acf[-1,,1][c(S, 2*S)]
if (R[1] <= 0)
  R[1] <- 0
if (R[2] <= 0)
  R[2] <- 0

The statistic is defined as follows ($R_s$ and $R_{2s}$ denote the autocorrelations obtained in the previous step):
$$
QS = n(n+2)\left( \frac{R_s^2}{n-s} + \frac{R_{2s}^2}{n-2s} \right) \,,
$$
where $n$ is the number of observations in the differenced series and $s$ is 
the periodicity of the data (12 in this case with monthly data). According to simulations exercises, this statistic follows approximately the $\chi^2$ distribuion with 2 degrees of freedom.
Thus, the statistic and the corresponding p-value can be obtained as follows:
n <- length(dx)
QS <- n*(n+2)*(R[1]^2/(n-S) + R[2]^2/(n-2*S))
pvalue <- pchisq(q=QS, df=2, lower.tail=FALSE)
round(c(QS=QS, p.value=pvalue), 4)
#       QS  p.value 
# 167.6486   0.0000 

which agrees with the output returned by seasonal::qs:
qs(m)["qsori",]
#             qs p-val
# qsori 167.6486     0

In order to apply this test to other series, e.g. qsorievadj, you just need 
to take this series from the output files and apply the operations shown above to it.
Given the large value of the QS statistic (and the implied low p-value), we can 
conclude that there is seasonality in the series.
