Dear all,
I have fitted a gls model where the response is a B-spline function of a
continuous variable (age), and variance is modelled with an exponential
correlation structure and a power function of the fitted values. I have
no problem to get the fitted values and I can compute the standard error
of the fitted mean. However, prediction fails: predict.gls does not seem
to work with B-spline transformations of the predictor. Is there a way
to solve this ? Furthermore, I would like to compute the standard error
of the predicted values (se of new observations in the range of the
observed values). How can I do this ?
Below are some code and results.
Thank you for your attention,
Renaud
> version
Professional Edition Version 6.0.3 Release 2 for Microsoft Windows :
2001
>
> summary(Data[, match(c("numani", "age", "poids"), table = names(Data))])
numani age poids
SNLOOV15809: 17 Min.: 0 Min.: 2.30
SNLOOV14696: 17 1st Qu.: 44 1st Qu.: 8.60
SNLOOV14053: 17 Median:115 Median:13.60
SNLOOV02497: 17 Mean:146 Mean:14.06
SNLOOV15464: 16 3rd Qu.:239 3rd Qu.:18.90
SNLOOV15177: 16 Max.:398 Max.:36.40
(Other):1545
> fm <- gls(poids ~ bs(age, df = 4),
+ correlation = corExp(form = ~ age | numani),
+ weights = varPower(form = ~ fitted(.)),
+ data = Data, method = "ML")
>
> summary(fm)
Generalized least squares fit by maximum likelihood
Model: poids ~ bs(age, df = 4)
Data: Data
AIC BIC logLik
7863.04 7906.284 -3923.52
Correlation Structure: Exponential spatial correlation
Formula: ~ age | numani
Parameter estimate(s):
range
0.9156445
Variance function:
Structure: Power of variance covariate
Formula: ~ fitted(.)
Parameter estimates:
power
1.152723
Coefficients:
Value Std.Error t-value p-value
(Intercept) 3.30796 0.0776681 42.59097 <.0001
bs(age, df = 4)1 6.06705 0.2474282 24.52045 <.0001
bs(age, df = 4)2 17.48736 0.6028137 29.00956 <.0001
bs(age, df = 4)3 14.85172 0.8443603 17.58931 <.0001
bs(age, df = 4)4 21.44970 0.7188427 29.83921 <.0001
Correlation:
(Intr) b(,d=4)1 b(,d=4)2 b(,d=4)3
bs(age, df = 4)1 -0.762
bs(age, df = 4)2 0.241 -0.615
bs(age, df = 4)3 -0.309 0.508 -0.779
bs(age, df = 4)4 0.010 -0.160 0.422 -0.640
Standardized residuals:
Min Q1 Med Q3 Max
-2.638735 -0.6878759 -0.09464405 0.5930081 4.360848
Residual standard error: 0.1418703
Degrees of freedom: 1645 total; 1640 residual
>
>## Fitted values
> Data$Fit <- fitted(fm)
> unique(Data[is.element(Data$age, c(30, 60, 90)),
+ match(c("age", "Fit"), table = names(Data))])
age Fit
180 90 13.231687
409 60 10.770752
670 30 7.519429
>
>## Predicted values
> data.frame(age = c(90, 60, 30),
Pred = predict(fm, newdata = data.frame(age = c(90, 60,
30))))
age Pred
1 90 24.75766
2 60 17.28133
3 30 3.30796
--
Dr Renaud Lancelot, vétérinaire
CIRAD, Département Elevage et Médecine Vétérinaire (CIRAD-Emvt)
Programme Productions Animales
http://www.cirad.fr/presentation/programmes/prod-ani.shtml (Français)
http://www.cirad.fr/presentation/en/program-eng/prod-ani.shtml (English)
ISRA-LNERV tel (221) 832 49 02
BP 2057 Dakar-Hann fax (221) 821 18 79 (CIRAD)
Senegal e-mail renaud.lancelot@cirad.fr
|