Wednesday, December 4, 2013

Fitting with weights, gnls, weights and form, in R



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