Grupo R : RMyP

Para aquellos que han leido Zuur (2009) “Mixed Effects Models and Extensions in Ecology with R” y en el capitulo 4 nos hace un plot de los predict del modelo sin intervalos de confianza al 95% y nos comenta que todavía no se ha desarrollado un código para calcularlos (predict.gls), pero que con un poco de código engorroso similar al de las ecuaciones de modelo lineales se podría obtener fácilmente. Ahí es donde por lo menos yo me quedé mirando la numeración de las páginas por si le faltaba alguna hoja, pero no, no se molesto en calcularlos ;)

Aquí esta la librería “rms”, que permite con la función Gls calcular los predict y graficarlos. El código cambia un poco pero esta basado en la librería “nlme”. Hay que fijarse en como se formulan los weights.

library(rms)
 
b3$BroodC <- NULL ##Eliminar una columna (datadist da problemas si en alguna columna hay NAs
dd <- datadist(b3); options(datadist='dd')  #opción de longitud necesaria para el Gls(  Gls {rms})
Form <- formula(log(Media) ~ logLc+fBrood2)
Ident<-varIdent( ~1 | fBrood2) 
M2<- Gls(Form,data=b3,weights = Ident)
 
###Ploteamos los residuales estandarizados
 
E <- resid(M, type = "normalized")
Fit <- fitted(M)
op <- par(mfrow = c(1, 3))
plot(x = Fit, y = E,
xlab = "Fitted values", ylab = "Residuals",
main = "Residuals versus fitted values")
text(Fit,E,b2$Cod,cex=0.6)
hist(E, nclass = 15)
plot(E ~fBrood2, data = b3, ylab = "Normalised residuals")
text(as.numeric(b3$fBrood2),E,b3$Cod,cex=0.6) ##Para saber identificar las observaciones
par(op) 
 
####Plot de los predict values and 95% confidence intervals in Gls
 
plot(Predict(M2))

This entry was posted on Friday, November 18th, 2011 at 10:57 am and is filed under Generalized Least Squares (NLS), Gráficos, Nonlinear Models. You can follow any responses to this entry through the RSS 2.0 feed. You can leave a response, or trackback from your own site.