raster predict with gnls model

Dear List,

I have trouble with the predict function (raster package) when using a gnls model (nonlinear model fit using generalized least squares, from the nlme package).


I'm using R version 3.3.1 (2016-06-21), Platform: x86_64-w64-mingw32/x64 (64-bit), and raster package version 2.5-8.
I want to make predictions based on three continuous variables and one categorical variable. The predictors are layers in a RasterStack object (see description below: s20, s20_excludefactor). Using other models than gnls works. I tested three models (my_nls, my_lm, my_gnls, see description below):

1) Using nonlinear least squares (nls) fit excluding the factor as predictor works fine:
v20 <- predict(object=s20_excludefactor, model=my_nls) 
2) Using linear model (lm) including the factor as predictor works as well, except the warning:
?not sure if the correct factor levels are used here?. 
v20 <- predict(object=s20, model=my_lm) #  or
v20 <- predict(object=s20, model=my_lm, na.rm=T, factors=list(Gruppe=c("1","2","3")))

3) Using the glns model produces error:
?Error in p %*% beta[pmap[[nm]]] : non-conformable arguments?
v20 <- predict(object=s20, model=my_gnls) #  or
v20 <- predict(object=s20, model=my_gnls, na.rm=T, factors=list(Gruppe=c("1","2","3")))


I don't know what to do about the error in 3). Any help is highly appreciated! 
Best regards,
Johannes

s20_excludefactor (used in my_nls)

class       : RasterStack 
dimensions  : 50, 50, 2500, 3  (nrow, ncol, ncell, nlayers)
resolution  : 20, 20  (x, y)
extent      : 3390000, 3391000, 5296000, 5297000  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=3500000 +y_0=0 +datum=potsdam +units=m +no_defs +ellps=bessel +towgs84=598.1,73.7,418.2,0.202,0.045,-2.455,6.7 
names       : x1, x2, x3 
min values  :      0.0000,      0.0000,    207.3915 
max values  : 243534.0084,      1.0000,    220.1788

s20 (used in my_lm and my_gnls)

class       : RasterStack 
dimensions  : 50, 50, 2500, 4  (nrow, ncol, ncell, nlayers)
resolution  : 20, 20  (x, y)
extent      : 3390000, 3391000, 5296000, 5297000  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=3500000 +y_0=0 +datum=potsdam +units=m +no_defs +ellps=bessel +towgs84=598.1,73.7,418.2,0.202,0.045,-2.455,6.7 
names       : x1, x2, x3,      Gruppe 
min values  :      0.0000,      0.0000,    207.3915,      1.0000 
max values  : 243534.0084,      1.0000,    220.1788,      3.0000

my_nls

Nonlinear regression model
model: y ~ a/(1 + exp(b0 + b1 * x1 + b2 * x2 + b3 * x3))
data: datapply
a         b0         b1         b2         b3 
1.500e+03  4.544e+00 -9.134e-06 -1.084e+00 -6.121e-04 
weighted residual sum-of-squares: 56.55

my_lm

Call:
lm(formula = y ~ x1 + x2 + x3 + Gruppe, data = datapply)
Coefficients:
(Intercept) x1 x2 x3      Gruppe2      Gruppe3 
8.208e+01    2.641e-03   -3.463e+02    8.357e-02    6.749e+01    8.759e+01

my_gnls

Generalized nonlinear least squares fit
Model: y ~ a/(1 + exp(b0 + b1 * x1 + b2 * x2 + b3 * x3)) 
Data: datapply 
Log-likelihood: -4305.82
Coefficients:
a             b0             b1             b2 b3.(Intercept)     b3.Gruppe2     b3.Gruppe3 
1.312412e+03   4.410668e+00  -1.001456e-05  -1.082276e+00  -7.200440e-05  -4.635569e-04  -5.423707e-04 
Variance function:
Structure: Power of variance covariate
Formula: ~x1 
Parameter estimates:
power 
0.9695105