Hi all
I want to construct a generalised nonlinear model (binomial family) using gnm,
of the form:
Response = a + b variable1 + c variable2 + d variable3 - d b variable4 - d c
variable5,
with the parameters b, c, and d appearing more than once. Hence, I think I
need to use a custom nonlin function with gnm.
One of my predictor variables is categorical, so I have created a dummy
variable for each factor level (i.e. zeros and ones), and incorporated each of
these variables into the model. I have previously done something similar with
nonlinear least squares models (nls) and it worked very well, however gnm
doesn't seem to like these dummy variables, and doesn't estimate the parameters.
It works if I put the dummy variables into a linear model using gnm, like so:
(where LithologyB is a series of ones and zeros, representing the occurrence of
Lithology B, etc.)
> summary(gnm(response~LithologyB+LithologyC+LithologyD+LithologyE+LithologyF,
> family=binomial))
Linear predictor - using glm.fit
Call:
gnm(formula = response ~ LithologyB + LithologyC + LithologyD + LithologyE
+ LithologyF, family = binomial)
Deviance Residuals:
Min 1Q Median 3Q Max
-17.676 -8.433 -1.181 7.825 23.415
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.026747 0.005900 -174.03 <2e-16 ***
LithologyB -0.715508 0.020844 -34.33 <2e-16 ***
LithologyC 0.272983 0.008588 31.79 <2e-16 ***
LithologyD -1.082133 0.020952 -51.65 <2e-16 ***
LithologyE 1.699259 0.007617 223.10 <2e-16 ***
LithologyF 1.505159 0.009932 151.54 <2e-16 ***
This gives the same result as putting in the original categorical Lithology
variable. But something goes wrong when I create the same model using a custom
nonlin function:
> gnm_function<- function(B1, B2, B3, B4, B5){
+ list(
+ predictors = list(B1=1, B2=1, B3=1, B4=1, B5=1),
+ variables = list(substitute(LithologyB), substitute(LithologyC),
substitute(LithologyD), substitute(LithologyE), substitute(LithologyF)),
+ term = function(predLabels, varLabels) {
+ paste(
+ predLabels[1], "*", varLabels[1], "+", predLabels[2], "*", varLabels[2], "+",
predLabels[3], "*", varLabels[3], "+", predLabels[4], "*", varLabels[4], "+",
predLabels[5], "*", varLabels[5]
+ )})}
> class(gnm_function) <- "nonlin"
> summary(gnm(response~gnm_function(B1, B2, B3, B4, B5), family=binomial))
Initialising
Running main iterations
Done
Call:
gnm(formula = response ~ gnm_function(B1, B2, B3, B4, B5), family =
binomial)
Deviance Residuals:
Min 1Q Median 3Q Max
-17.676 -9.787 -1.347 8.547 17.946
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.722412 0.003469 -208.2 <2e-16 ***
B1 NA NA NA NA
B2 NA NA NA NA
B3 NA NA NA NA
B4 1.394925 0.005936 235.0 <2e-16 ***
B5 NA NA NA NA
Can anyone suggest what is going wrong and a possible way around the problem?
Thanks
Brett Murphy
ARC Postdoctoral Research Fellow
School of Plant Science
University of Tasmania
c/o CSIRO Tropical Ecosystems Research Centre
PMB 44, Winnellie, NT 0822
Ph: (08) 8944 8447
Fax: (08) 8944 8444
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
