What @Rolan is getting at is to use a random spline basis for the time by id random part. So your model would become:

m <- gam(y ~ s(V1) + V2 + s(time, id, bs = 'fs'),
         family=gaussian, data=dat, method = "REML")

This model says that the effect of time is smooth and varies by id, with a separate smooth being estimated for each id but each smooth is assumed to have the same wiggliness (a single smoothness penalty is estimated for all the time smoothers) but can differ in shape.

To estimate a separate "global" effect the model could be

m <- gam(y ~ s(V1) + V2 + s(time) + s(time, id, bs = 'fs'),
         family=gaussian, data=dat, method = "REML")

If you want similar models but where each smooth can have different wiggliness as well as shape, the the by smoothers can be used:

## without a "global" effect
m <- gam(y ~ s(V1) + V2 + s(id, bs = 're') + s(time, by = id),
         family=gaussian, data=dat, method = "REML")

## with a "global" effect
m <- gam(y ~ s(V1) + V2 +
         s(id, bs = 're') + s(time) + s(time, by = id, m = 1),
         family=gaussian, data=dat, method = "REML")

The m=1 means that the smoother uses a penalty on the squared first derivative, which penalises departure from a flat function of no effect. As this is on the subject specific smooths, the model is penalising deviations from the "global" smooth.

Some colleagues and I have described these models in some detail in a paper submitted to PeerJ, which is available as a preprint. A new version in response to reviewers comments should be up in a few days (we've submitted it to the journal).

What @Roland is getting at is to use a random spline basis for the time by id random part. So your model would become:

m <- gam(y ~ s(V1) + V2 + s(time, id, bs = 'fs'),
         family=gaussian, data=dat, method = "REML")

This model says that the effect of time is smooth and varies by id, with a separate smooth being estimated for each id but each smooth is assumed to have the same wiggliness (a single smoothness penalty is estimated for all the time smoothers) but can differ in shape.

To estimate a separate "global" effect the model could be

m <- gam(y ~ s(V1) + V2 + s(time) + s(time, id, bs = 'fs'),
         family=gaussian, data=dat, method = "REML")

If you want similar models but where each smooth can have different wiggliness as well as shape, then the by smoothers can be used:

## without a "global" effect
m <- gam(y ~ s(V1) + V2 + s(id, bs = 're') + s(time, by = id),
         family=gaussian, data=dat, method = "REML")

## with a "global" effect
m <- gam(y ~ s(V1) + V2 +
         s(id, bs = 're') + s(time) + s(time, by = id, m = 1),
         family=gaussian, data=dat, method = "REML")

The m=1 means that the smoother uses a penalty on the squared first derivative, which penalises departure from a flat function of no effect. As this is on the subject specific smooths, the model is penalising deviations from the "global" smooth.

Some colleagues and I have described these models in some detail in a paper submitted to PeerJ, which is available as a preprint. A new version in response to reviewers comments should be up in a few days (we've submitted it to the journal).