Problems interpreting GAM output I have been advised to run General Additive Models to be able to describe trends in my data, my data being animal harvest numbers by year. I have done so, but have a problem with interpreting the model output. Hopefully it will be sufficient to show you parts of my script, and then three of my graphs with corresponding model output.
The synthax below, fitting both a linear and a nonlinear component:
library(mgcv)
model1 <- gam(Tot~Year+s(Year),family=poisson,data=ds)
model2 <- gam(Tot~Year+s(Year),family=poisson,data=ls)
model3 <- gam(Tot~Year+s(Year),family=poisson,data=nw)

Then I get the summary output for the models:
summary(model1)
summary(model2)
summary(model3)

Then I group the data in 5-year periods, before generating the data for the plots:
YearP=seq(1975,2015,by=5)
model1.pred=predict(model1,newdata=data.frame(Year=YearP),type="response",se.fit=T)
model2.pred=predict(model2,newdata=data.frame(Year=YearP),type="response",se.fit=T)
model3.pred=predict(model3,newdata=data.frame(Year=YearP),type="response",se.fit=T)

I graph the three models in the similar manner, showing both the data and the model output:
library(gplots)
plotCI(x=YearP, y=model1.pred$fit,uiw=2*model1.pred$se.fit,
    col="red",lwd=3,cex=1.2,las=1, 
    xlab="", ylab="Observed and fitted numbers")
points(ds$Year,ds$Tot,pch=19,cex=0.9)
text(1975,60,label="GRAPH 1",cex=1.4,adj=0)

Here's my three graphs xlab is years and ylab is "observed and fitted numbers", the data are the black dots and the model predictions are the red bars:

And here are the outputs for the three models:

I am trying to use the best method possible to explain my data, and I believe this is it. I have tried to google this, and I have read the relevant answers at CrossValidated, but have yet to find explanations that make sense to me. I need someone to tell me what the output means. 
I specifically need help to understand the estimate of the parametric coefficient and the p-value, what do they tell me? E.g. in graph 1 there is obviously an increase, and in graph 3 a decrease, but coefficient estimates for both are positive. This is difficult to understand.
Then I also need to know what the "Approximate significance of the smooth terms" tell me, especially the edf and the p-value here?
Extremely grateful for help to understand.
 A: Your model is performing a regression through the origin -- the (Intercept) term is 0.  It is unclear why mgcv is doing this, since as is shown below in an example, in the Gaussian case there is no difficulty to include both parametric and smoothed versions of a covariate.
The fact that the model, as fit, is regression through the origin explains why the Year coefficients as are all positive: at year 0, the poisson log-rate is 0, or 1 events occur per year.  Given that constraint, by year 1980, where the observed rate is approximately 30 events per year, the parametric "year" coefficients are all positive, and approximately $10^{-3}$, since $\exp(2000 \times 0.0013) \approx 35$.
It's important to understand how the intercept and parametric year coefficient are identified. The family of smooth functions generated by s certainly includes linear functions.  Identification of the parametric terms is achieved by constraining the smooth (s(Year)) term to sum to zero over the observed values of the covariates.  This can lead to some counter-intuitive behavior.  Consider two smooth functions with a linear component:
$$
\begin{align*}
Y_1 & = 4 + 2x + \sin(4\pi x), \\
Y_2 & = 4 + 2x + \cos(4\pi x).
\end{align*}
$$
In both cases, the non-linear part of the function integrates to zero for $x \in [0,1]$.  But when we fit the gams (dashed lines), and compare the parametric portions to the estimate you'd get from just estimating the linear trend in y (with lm, solid lines), for the sin function, gam overestimates the linear trend, both in terms of the coefficients of the true function, and compared to the (mis-specified) least-squares linear trend.    For the cos function, the parametric estimates are within a standard error of the true forms, and the linear part of the gam matches the linear trend estimated from the least-squares linear trend.  I don't have a good explanation for why.


Code
ex = data_frame(x = seq(0, 1, length.out = 100)) %>% mutate(
    epsilon = rnorm(100, sd = .05),
    y1 = 4 + 2 * x + sin(x * 4 * pi) + epsilon,
    y2 = 4 + 2 * x + cos(x * 4 * pi) + epsilon
)

plot(ex$x, ex$y1, col = 'black')
points(ex$x, ex$y2, col = 'red')

gfit1 = gam(y1 ~ 1 + x + s(x), data = ex)
gfit2 = gam(y2 ~ 1 + x + s(x), data = ex)
lin1 = lm(y1 ~ x, data = ex)
lin2 = lm(y2 ~ x, data = ex)
gam_linear = function(fit) {
    mm = model.matrix(fit)[, c('(Intercept)', 'x')]
    beta = coef(fit)[c('(Intercept)', 'x')]
    mm %*% beta
}

ex = ex %>% mutate(
    p1 = predict(gfit1),
    p2 = predict(gfit2),
    glin1 = gam_linear(gfit1),
    glin2 = gam_linear(gfit2),
    lin1 = predict(lin1),
    lin2 = predict(lin2)
)

lines(ex$x, ex$p1, col = 'black', lty = 2)
lines(ex$x, ex$p2, col = 'red', lty = 2)
lines(ex$x, ex$glin1, col = 'black', lty = 2)
lines(ex$x, ex$glin2, col = 'red', lty = 2)

lines(ex$x, ex$lin1, col = 'black', lty = 1)
lines(ex$x, ex$lin2, col = 'red', lty = 1)
legend(
    'topleft',
    legend = c(
        '4 + 2 * x + sin(x * 4 * pi) (gam)',
        '4 + 2 * x + cos(x * 4 * pi) (gam)',
        '4 + 2 * x + sin(x * 4 * pi) (lm)',
        '4 + 2 * x + cos(x * 4 * pi) (lm)'
    ),
    col = rep(c('black', 'red'), 2),
    lty = c(2, 2, 1, 1)
)


A: I'm not entirely clear on why your intercept is 0 for each of your models, as Andrew M pointed out.  The default gam() function provides an intercept, which is probably why your intercept did not change when you added + 1 to your model formula.
Why do you feel the need to specify both a parametric term and a smooth term for Year?  Including both terms in your model seems to complicate the interpretation unnecessarily.  A smooth term for Year by itself would allow you to visualize both the general trend in your data as well as any anomalies.  I may be missing something important regarding your field, but based off of what I know so far, I would specify my model as follows:
model <- gam(Tot ~ s(Year, k = 10), family = poisson, data = data)

Since you are trying to understand trends over the course of 40 years, I would also pay attention to the value of k for your smooth term.  k specifies the maximum degrees of freedom your smooth can take (which is k-1), and the gam() function automatically constrains k to be reasonably flexible given the number of observations you have.  The main thing is that you don't want k to be too small.  This can be easily checked with k.check(model).
The easiest way to display the resulting smooth (and hopefully draw conclusions from it) is to use the plot.gam() function:
plot(model, shaded = T)

This will plot the smooth terms in your model with shaded standard errors.
In my opinion, smooth terms are generally difficult to make hard conclusions about.  It's not like a parametric term where you can make a  statement like: "animal harvests increased by 15% a year."  However, if the association with a response is highly nonlinear, smooth terms can provide highly intuitive visualizations.  And if you desire to make bold, concrete statements about the predictor, you can always use smooths as a preliminary analysis to a parametric model, where you could specify higher order terms as needed.
