gnls with weights (form)
http://stackoverflow.com/questions/10508474/step-halving-issue-in-gnlsnlme
dat:
Site_Code SF
5 3
5 0
5 2
5 0
5 0
5 0
library(nlme)
g0 <- gnls(SF ~ a * Site_Code^b, data = dat,
weights = varPower(form = ~Site_Code),
start=list(a=30,b=-0.5))
This example indicates that averaged data should NOT be used for regression if weights
are used. Instead, original data should be used, because the noises in the un-averaged
data can be used as weights.
The following example shows varPower() can bring the fitting model closer to 'truth'.
> require(nlme)
> x1 = rnorm(20)
> y1 = x1 + rnorm(20)/10
>
> #weird ones
> x2 = rnorm(5)+5
> y2 = rnorm(5)
>
> y=c(y1,y2)
> x=c(x1,x2)
> mydata = data.frame(cbind(y,x))
>
> foo = function(a,b) { y = x*a + b }
> model1 = gnls( y ~foo(a,b), start=list(a=1,b=0))
> summary(model1)
Generalized nonlinear least squares fit
Model: y ~ foo(a, b)
Data: NULL
AIC BIC logLik
72.21925 75.87587 -33.10962
Coefficients:
Value Std.Error t-value p-value
a -0.0159375 0.08311006 -0.1917637 0.8496
b -0.3550085 0.20216324 -1.7560486 0.0924
Correlation:
a
b -0.346
Standardized residuals:
Min Q1 Med Q3 Max
-1.8712584 -0.7939610 0.2684611 0.7910997 1.2959748
Residual standard error: 0.9485101
Degrees of freedom: 25 total; 23 residual
>
> model2 = gnls( y ~foo(a,b), data=mydata, start=list(a=1,b=0), weights=varPower(form = ~x))
> summary(model2)
Generalized nonlinear least squares fit
Model: y ~ foo(a, b)
Data: mydata
AIC BIC logLik
20.85874 25.73424 -6.429369
Variance function:
Structure: Power of variance covariate
Formula: ~x
Parameter estimates:
power
1.468287
Coefficients:
Value Std.Error t-value p-value
a 0.9281320 0.06104714 15.203530 0.0000
b -0.0255823 0.00891448 -2.869751 0.0087
Correlation:
a
b 0.17
Standardized residuals:
Min Q1 Med Q3 Max
-2.0273323 -0.8480791 -0.1211550 0.3073090 2.2141947
Residual standard error: 0.4204751
Degrees of freedom: 25 total; 23 residual
> model1
Generalized nonlinear least squares fit
Model: y ~ foo(a, b)
Data: NULL
Log-likelihood: -33.10962
Coefficients:
a b
-0.0159375 -0.3550085
Degrees of freedom: 25 total; 23 residual
Residual standard error: 0.9485101
> model2
Generalized nonlinear least squares fit
Model: y ~ foo(a, b)
Data: mydata
Log-likelihood: -6.429369
Coefficients:
a b
0.92813202 -0.02558234
Variance function:
Structure: Power of variance covariate
Formula: ~x
Parameter estimates:
power
1.468287
Degrees of freedom: 25 total; 23 residual
Residual standard error: 0.4204751
No comments:
Post a Comment