## 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 ``````
`` ``