#====================================================== # Descricao do Estudo #======================================================= # # Coorte aberta de mes HIV positiva - Desenho desbalanceado # # Resposta medida: comprimento do Recem-nascido nos primeiros # 18 meses de vida # # Covariaveis: Sexo (0 menina, 1- menino) # Status (0 sororevertor,1- infectado) # #====================================================== # Leitura Geral #======================================================= # rm(list=ls(all=TRUE)) # limpar a console require(nlme) require(growth) source("longit.R") #scripts do R dados<-read.table("alturalong.txt",header=T) # ler dados # attach(dados) names(dados) dim(dados) dados[1:15,] # ja esta no formato longo # #======================================================== # Manipulando Banco de Dados Longitudinais # Banco Longo (uma medida por linha) #========================================================= # length(unique(dados$ident)) # qtos individuos no estudo as.data.frame(table(dados$ident)) #mostra qtas vezes cada individuo aparece max(table(dados$ident)) #numero maximo de medidas por individuo min(table(dados$ident)) #numero minimo de medidas por individuo # crianca386<-subset(dados,ident==386) # seleciona crianca com minimo crianca386 # tam<-function(x) {length(unique(x))} tam(dados$status) tapply(dados$ident, dados$status, tam) #encontra o numero de individuos por status tapply(dados$ident, dados$sexo, tam) #encontra o n?mero de indiv?duos por sexo # median(tapply(dados$Idade,list(dados$ident), max)) #tempo mediano de acompanhamento min(tapply(dados$Idade,list(dados$ident), max)) # tempo m?nimo de acompanhamento median(tapply(dados$Idade[dados$sexo==0],list(dados$ident[dados$sexo==0]),max)) # mediano pra sexo=0 # datalong <- as.ldframe(dados, ident, Idade) class(datalong) # attach(datalong) freq <- as.data.frame(table(datalong$Subject)) freq #mostra qtas vezes cada individuo aparece hist((table(dados$ident)),xlab="number of measurements", main = paste("")) summary(freq[, 2]) hist(freq[, 2]) # #================================================================== # Graficos de Perfis #================================================================== plot(datalong, Altura, n.full = 0) # #attach(dados) # function(choice,xl=range(time),yl=range(y)) # {# Plots profiles for subjects in vector choice. m=length(unique(dados$ident)) n=length(dados$Altura)/m xl=range(dados$Idade) yl=range(dados$Altura) choice=dados$ident j=dados$ident==choice[1] plot(dados$Idade[j],dados$Altura[j],n.full=0, type="l",ylim=yl,xlim=xl,ylab="Altura",xlab=" Idade", main="Perfis das Criancas") if(length(choice)>1){ for(it in 2:length(choice)){ j=dados$ident==choice[it] lines(dados$Idade[j],dados$Altura[j],n.full=0)} } cat("\n") # } # Fazendo um alisamento no grafico de perfis lines(lowess(dados$Idade,dados$Altura), col = 2, type="l",ylab="Altura",xlab="Idade") #======================================================================== # Graficos Alisados para Meninos e Meninas #======================================================================== par(mfrow=c(1,2)) meninos<-subset(dados,sexo==0) # attach(meninos) plot(lowess(meninos$Idade,meninos$Altura),type="l",ylab="height (cms)", xlab="age (months)", main="Pefis por Sexo" ) meninas<-subset(dados,sexo==1) attach(meninas) lines(lowess(meninas$Idade,meninas$Altura),lty=2,type="l") legend(0, 78, legend = c("girls","boys"), lty = c(1,2), bty = "n", cex = 0.8) # # grupo1<-subset(dados,status==0) # attach(grupo1) plot(lowess(grupo1$Idade,grupo1$Altura),type="l",ylab="height (cms)", xlab="age (months)", # main="Pefis por Grupo" ) grupo2<-subset(dados,status==1) # attach(grupo2) lines(lowess(grupo2$Idade,grupo2$Altura),lty=2,type="l") legend(0, 78, legend = c("uninfected","infected"), lty = c(1,2), bty = "n", cex = 0.8) # #==================================================================== # Estrutura da Media #======================================================================= # Estrutura1: Termo Quadratico # # Idade centrada na media dados$Idade1<-dados$Idade-7 # # Termo quadratico dados$Idade2<-dados$Idade1*dados$Idade1 # # Estrutura2: Regressao com duas inclinacoes # em tempo = 5 meses # S<-length(dados$Idade) dados$Age<-rep(0,S) for(i in 1:S) { if(dados$Idade[i]>5) {dados$Age[i]<-dados$Idade[i]-5} } # #========================================================================= # Estrutura de covariancia/correlacao # usando o variograma #========================================================================= # Modelo polinomial cheio modelo11<-gls(Altura ~ Idade+I(Idade^2)+sexo+status+Idade:status+ I(Idade^2):status+Idade:sexo+I(Idade^2):sexo,data=dados) summary(modelo11) plot(Variogram(modelo11,form=~Idade,resType="r", robust=T)) # modelo12<-gls(Altura ~ Idade+I(Idade^2)+sexo+status+Idade:status+ I(Idade^2):status+Idade:sexo+I(Idade^2):sexo,correlation= corCompSymm(form=~1|ident),data=dados) plot(Variogram(modelo12,form=~Idade,resType="r",robust=T)) # # Modelo com Duas Inclinacoes modelo21<-gls(Altura ~ Idade+Age+sexo+status+Idade:status+Age:status+ Idade:sexo+Age:sexo,data=dados) plot(Variogram(modelo21,form=~Idade,resType="r", robust=T)) # modelo22<-gls(Altura ~ Idade+Age+sexo+status+Idade:status+Age:status+ Idade:sexo+Age:sexo,correlation= corCompSymm(form=~1|ident),data=dados) plot(Variogram(modelo22,form=~Idade,resType="r", robust=T)) # # #========================================================================== # Modelos de Efeitos Fixos com Estrutura de Covariancia Exp(|t_j - t_l|) # e outras estruturas espaciais: AR1, Gauss, Esf?rica. # EMVR #=========================================================================== # ## MODELO 1 ## Modelo linear de efeito fixo com duas inclinacoes e especificao da ## estrutura de correlao # # Exponencial modelo1<-gls(Altura ~ Idade+Age+sexo+status+Idade:status+Age:status+ Idade:sexo+Age:sexo,correlation=corExp(form=~Idade|ident),data=dados) summary(modelo1) plot(modelo1) plot(modelo1,resid(.,type="p")~fitted(.)|status,id=0.05) plot(modelo1,Altura~fitted(.),adj=-0.3) # predito vs obs. plot(Variogram(modelo1,form=~Idade,resType="n", robust=T)) # AR1 continuo modelo11<-gls(Altura ~ Idade+Age+sexo+status+Idade:status+Age:status+ Idade:sexo+Age:sexo,correlation=corAR1(form=~Idade|ident),data=dados) summary(modelo11) # Gaussiano modelo12<-gls(Altura ~ Idade+Age+sexo+status+Idade:status+Age:status+ Idade:sexo+Age:sexo,correlation=corGaus(form=~Idade|ident),data=dados) summary(modelo12) # Spherica modelo13<-gls(Altura ~ Idade+Age+sexo+status+Idade:status+Age:status+ Idade:sexo+Age:sexo,correlation=corSpher(form=~Idade|ident),data=dados) summary(modelo13) # # # Retirando os termos sem signific?ncia # modelo2<-gls(Altura ~ Idade+Age+sexo+status+Idade:status+Age:status+ Idade:sexo,correlation=corExp(form=~Idade|ident),data=dados) summary(modelo2) plot(modelo2) plot(Variogram(modelo2,form=~Idade+Age,resType="n", robust=T)) # # Retirando os termos sem signific?ncia # modelo3<-gls(Altura ~ Idade+Age+sexo+status+Idade:status+Idade:sexo, correlation=corExp(form=~Idade|ident),data=dados) summary(modelo3) plot(modelo3) plot(modelo3,resid(.,type="p")~fitted(.)|status,id=0.05,adj=-0.3) plot(modelo3,Altura~fitted(.),adj=-0.3) # predito vs obs. plot(Variogram(modelo3,form=~Idade+Age,resType="n", robust=T)) # ## MODELO 2 ## Modelo linear de efeito fixo com termo quadratico e especificacao da ## estrutura de correlacao # # Modelo OLS # modelo0<-gls(Altura ~ Idade1+I(Idade1^2)+sexo+status+Idade1:status+ I(Idade1^2):status+Idade1:sexo+I(Idade1^2):sexo,data=dados) summary(modelo0) plot(modelo0) plot(Variogram(modelo0,form=~Idade1,resType="n", robust=T)) # Exponencial modelo4<-gls(Altura ~ Idade1+I(Idade1^2)+sexo+status+Idade1:status+ I(Idade1^2):status+Idade1:sexo+I(Idade1^2):sexo,correlation= corExp(form=~Idade1|ident),data=dados) summary(modelo4) plot(modelo4) plot(modelo4,resid(.,type="p")~fitted(.)|status) plot(modelo4,Altura~fitted(.),adj=-0.3) # predito vs obs. plot(Variogram(modelo4,form=~Idade1,resType="n", robust=T)) modelo41<-update(modelo4,weights=varPower(),data=dados) summary(modelo41) anova(modelo4,modelo41) plot(modelo41) plot(modelo41,resid(.,type="p")~fitted(.)|status,id=0.05,adj=-0.3) plot(modelo41,Altura~fitted(.),adj=-0.3) # predito vs obs. plot(Variogram(modelo41,form=~Idade1,resType="n", robust=T)) # modelo5<-gls(Altura ~ Idade+I(Idade^2)+sexo+status+Idade:status+ I(Idade^2):status+Idade:sexo+I(Idade^2):sexo,correlation= corCompSymm(form=~1|ident),data=dados) summary(modelo5) # #========================================================= # Analise de residuos #========================================================== res = modelo4$residuals pred = modelo4$fitted # Histogramas com curva normal x = seq(min(res)-1,max(res)+1,by=1) hist(res,prob=T,xlab="Res?duos",ylab="Densidade",main="") lines(x,dnorm(x,mean=mean(res),sd=sd(res))) # Gr?fico de Probabilidade Normal dos Res?duos qqnorm(res,xlab="Quantis da Normal Padr?o",ylab="Quantis dos Res?duos",main="") qqline(res) # # Gr?fico de dispers?o dos Res?duos versus Valores ajustados plot(pred,res,xlab="Valores preditos",ylab="Res?duos") abline(h=0) lines(lowess(pred,res),type="l",lwd=2,col="red") # # Gr?fico de dispers?o dos Res?duos versus Vari?vel Tempo plot(dados$Idade-7,res,xlab="Idad",ylab="Res?duos") abline(h=0) lines(lowess(Idade,res),type="l",lwd=2,col="red") # Idade_centrada<-dados$Idade-7 plot(Variogram(modelo4,form=~Idade_centrada,resType="n", robust=T)) # #===================================================================== # Gr?ficos dos Preditos #======================================================================= fitted<-modelo4$fitted Idade<-dados$Idade-7 #out1<-lme(Altura~Idade+I(Idade^2)+sexo+status+Idade*sexo+Idade*status, random= ~Idade+I(Idade^2)|ident,method="ML") # equivalente #fitted1<-out1$fitted[,1] #dados1<-data.frame(dados,fitted,fitted1) dados1<-data.frame(dados,fitted) meninos<-subset(dados1,sexo==0) attach(meninos) par(mfrow=c(1,2)) plot(lowess(meninos$Idade,meninos$Altura),type="l",col=2,ylab="Altura (cms)", xlab="Idade (m?s)",main="G?nero") lines(lowess(meninos$Idade,meninos$fitted),lty=2,col=3,type="l") #lines(lowess(meninos$Idade,meninos$fitted1),lty=3,col=4,type="l") meninas<-subset(dados1,sexo==1) attach(meninas) lines(lowess(meninas$Idade,meninas$Altura),type="l",col=2) lines(lowess(meninas$Idade,meninas$fitted),lty=2,col=3,type="l") #lines(lowess(meninas$Idade,meninas$fitted1),lty=3,col=4,type="l") casos<-subset(dados1,status==0) attach(casos) plot(lowess(casos$Idade,casos$Altura),type="l",col=2,ylab="Altura (cms)", xlab="Idade (m?s)",main="Grupo") lines(lowess(casos$Idade,casos$fitted),lty=2,col=3,type="l") control<-subset(dados1,status==1) attach(control) lines(lowess(control$Idade,control$Altura),type="l",col=2) lines(lowess(control$Idade,control$fitted),lty=2,col=3,type="l") # #==================================================================== # Fazendo os Ajustes de Modelos de Efeito Aleatorio #==================================================================== # # Modelo Inadequado (linear para a m?dia) out<-lme(Altura~Idade+sexo+status, random= ~1|ident,data=dados) summary(out) plot(out) # mostra a inadequacao do modelo linear # # Modelo Inadequado (linear pra m?dia) out<-lme(Altura~Idade+sexo+status+Idade:sexo+Idade:status, random= ~Idade|ident,data=dados) plot(out) # # Modelo com Duas Inclina??es # modelo1<-lme(Altura ~ Idade+Age+sexo+status+Idade:status+Age:status+ Idade:sexo+Age:sexo,random=~Idade+Age|ident,data=dados) summary(modelo1) plot(modelo1) plot(Variogram(modelo1,form=~Idade+Age,resType="n", robust=T)) modelo11<-update(modelo1,weights = varPower()) summary(modelo11) plot(modelo11) plot(modelo11,resid(.,type="p")~fitted(.)|status,id=0.05,adj=-0.3) plot(modelo11,Altura~fitted(.),adj=-0.3) # predito vs obs. plot(Variogram(modelo11,form=~Idade,resType="n", robust=T)) anova(modelo1,modelo11) # # Modelo Quadr?tico # out<-lme(Altura~Idade+I(Idade^2)+sexo+status+Idade:sexo+Idade:status+ I(Idade^2):sexo+I(Idade^2):status, random= ~Idade+I(Idade^2)| ident,data=dados) plot(out) summary(out) plot(Variogram(out,form=~Idade,resType="n", robust=T)) out1<-lme(Altura~Idade+I(Idade^2)+sexo+status+Idade*sexo+Idade*status+ I(Idade^2)*sexo+I(Idade^2)*status, random= ~Idade|ident,data=dados) plot(out1) plot(Variogram(out1,form=~Idade,resType="n", robust=T)) anova(out,out1) plot(out,Altura~fitted(.),adj=-0.3) # predito vs obs. plot(Variogram(out, form = ~ Idade | ident, nint = 20, robust = TRUE)) out2<-lme(Altura~Idade+I(Idade^2)+sexo+status+Idade:sexo+Idade:status+ I(Idade^2):sexo+I(Idade^2):status, random= ~Idade+I(Idade^2)| ident,weights=varPower(),data=dados) plot(out2) summary(out2) plot(out2,Altura~fitted(.),adj=-0.3) # predito vs obs. plot(Variogram(out2, form = ~ Idade | ident, nint = 20, robust = TRUE)) anova(out,out2) # #========================================================= # An?lise de res?duos #========================================================== # # Res?duos n?o-transformados res = out$residuals[,1] pred = out$fitted[,1] length(res) # # Res?duos transformados g11 = 2.13539904^2 g22 = 0.56059959^2 g33 = 0.03098822^2 g12 = -0.129*2.13539904*0.56059959 g13 = 0.003 g23 = -0.918*0.56059959*0.03098822 # G = matrix(c(g11,g12,g13,g12,g22,g23,g13,g23,g33),3,3,byrow=T) sigma2 = (out$sigma)^2 # Z1 = rep(1,1318) Z2 = dados$Idade-7 Z3 = I((dados$Idade-7)^2) ID<-dados$ident Z = cbind(ID,Z1,Z2,Z3) resID = cbind(res,ID) predID = cbind(pred,ID) timeID = cbind(dados$Idade-7,ID) rest = NULL predt = NULL timet = NULL # j=0 for(i in 1:139){ j=unique(dados$ident)[i] a<-dim(subset(dados,dados$ident==j)) i<-ifelse(a[1]>1,i,i+1) j=unique(dados$ident)[i] covYi = Z[ID==j,2:4]%*%G%*%t(Z[ID==j,2:4])+sigma2*diag(sum(ID==j)) rest = c(rest,as.vector(solve(chol(covYi))%*%resID[ID==j,1])) predt = c(predt,as.vector(solve(chol(covYi))%*%predID[ID==j,1])) timet = c(timet,as.vector(solve(chol(covYi))%*%timeID[ID==j,1])) } # # Histogramas com curva normal par(mfrow=c(1,2)) x = seq(min(res)-1,max(res)+1,by=1) hist(res,prob=T,xlab="Res?duos",ylab="Densidade",main="") lines(x,dnorm(x,mean=mean(res),sd=sd(res))) x = seq(min(rest)-1,max(rest)+1,by=0.05) hist(rest,prob=T,xlab="Res?duos transformados",ylab="Densidade",main="") lines(x,dnorm(x,mean=mean(rest),sd=sd(rest))) # # Gr?fico de Probabilidade Normal dos Res?duos par(mfrow=c(1,2)) qqnorm(res,xlab="Quantis da Normal Padr?o",ylab="Quantis dos Res?duos",main="") qqline(res) qqnorm(rest,xlab="Quantis da Normal Padr?o",ylab="Quantis dos Res?duos Transformados",main="") qqline(rest) # # Gr?fico de dispers?o dos Res?duos versus Valores ajustados par(mfrow=c(1,2)) plot(pred,res,xlab="Valores preditos",ylab="Res?duos") abline(h=0) lines(lowess(pred,res),type="l",lwd=2,col="red") plot(predt,rest,xlab="Valores preditos transformados",ylab="Res?duos transformados") abline(h=0) lines(lowess(predt,rest),type="l",lwd=2,col="red") # # Gr?fico de dispers?o dos Res?duos versus Vari?vel Tempo par(mfrow=c(1,2)) plot(Idade,res,xlab="Idad",ylab="Res?duos") abline(h=0) lines(lowess(Idade,res),type="l",lwd=2,col="red") plot(timet,rest,xlab="Tempo transformado",ylab="Res?duos transformados") abline(h=0) lines(lowess(timet,rest),type="l",lwd=2,col="red") # Idade_centrada<-dados$Idade-7 plot(Variogram(out,form=~Idade_centrada,resType="n", robust=T)) # #===================================================================== # Gr?ficos dos Preditos #======================================================================= fitted<-out$fitted[,1] Idade<-Idade0-7 #out1<-lme(Altura~Idade+I(Idade^2)+sexo+status+Idade*sexo+Idade*status, random= ~Idade+I(Idade^2)|ident,method="ML") # equivalente #fitted1<-out1$fitted[,1] #dados1<-data.frame(dados,fitted,fitted1) dados1<-data.frame(dados,fitted) meninos<-subset(dados1,sexo==0) attach(meninos) par(mfrow=c(1,2)) plot(lowess(meninos$Idade,meninos$Altura),type="l",col=2,ylab="Altura (cms)", xlab="Idade (m?s)",main="G?nero") lines(lowess(meninos$Idade,meninos$fitted),lty=2,col=3,type="l") #lines(lowess(meninos$Idade,meninos$fitted1),lty=3,col=4,type="l") meninas<-subset(dados1,sexo==1) attach(meninas) lines(lowess(meninas$Idade,meninas$Altura),type="l",col=2) lines(lowess(meninas$Idade,meninas$fitted),lty=2,col=3,type="l") #lines(lowess(meninas$Idade,meninas$fitted1),lty=3,col=4,type="l") casos<-subset(dados1,status==0) attach(casos) plot(lowess(casos$Idade,casos$Altura),type="l",col=2,ylab="Altura (cms)", xlab="Idade (m?s)",main="Grupo") lines(lowess(casos$Idade,casos$fitted),lty=2,col=3,type="l") control<-subset(dados1,status==1) attach(control) lines(lowess(control$Idade,control$Altura),type="l",col=2) lines(lowess(control$Idade,control$fitted),lty=2,col=3,type="l") # #===================================================================== # GAM #======================================================================= require(mgcv) attach(dados) b<-gam(Altura~s(Idade)+status) # sem efeito aleat?rio + spline idade plot(b) summary(b) plot(b,pages=1,residuals=TRUE) plot(b,pages=1,seWithMean=TRUE) anova(b) b<-gam(Altura~s(Idade,bs="cr")+status) # cubic spline b<-gam(Altura~s(Idade,k=20)+status) # n?mero de basis, default 10 b<-gam(Altura~s(Idade)+status, gamma=1.4) # controla a suaviaza??o, default=1 plot(b,all,terms=T) anova(b) dados$fem <- -(dados$sexo-1) b.sexo<-gam(Altura~s(Idade,by=sexo)+s(Idade,by=dados$fem)) # sem efeito aleat?rio + spline momento par(mfrow=c(1,2)) plot(b.sexo) # # # Estratificado por G?nero e Grupo # dados$fem <- -(dados$sexo-1) b.sexo<-gam(Altura~s(Idade,by=sexo,k=5)+s(Idade,by=dados$fem,k=5)) # sem efeito aleat?rio + spline momento par(mfrow=c(1,2)) #plot(b.sexo,shade=T) plot(b.sexo,ylab="Altura") dados$grupo <- -(dados$status - 1) b.grupo<-gam(Altura~s(Idade,by=status,k=6)+s(Idade,by=dados$grupo,k=6)) # sem efeito aleat?rio + spline momento par(mfrow=c(1,2)) plot(b.grupo,ylab="Altura") # # NAO ESTA FUNCIONANDO A FUNCAO GAMM---> Problemas na estimacao... # a<-gamm(Altura~s(Idade)+sexo+status+sexo:Idade+status:Idade+I(Idade^2):sexo+I(Idade^2):status,random=list(ident=~Idade+I(Idade^2)),data=dados) a<-gamm(Altura~s(Idade,k=5)+sexo+status,knots=list(c(3,6,9)), random=list(ident=~1),data=dados) summary(a) summary(a$lme) plot(a$lme) summary(a$gam) anova(a$gam) plot(a$gam) # # #==================================================================== # Fazendo os Ajustes de Modelos de Efeito Aleatorio # Resposta Transformada para SQRT(idade + 1) #==================================================================== # require(nlme) attach(dados) # #Avaliando a transforma??o raiz + 1 out<-gamm(Altura~s(sqrt(Idade+1))+status+sexo+status:sqrt(Idade+1)+sexo:sqrt(Idade+1),random=list(ident=~1)) plot(out$gam) out<-gamm(Altura~s(sqrt(Idade+1))+status+sexo+status:sqrt(Idade+1)+sexo:sqrt(Idade+1),random=list(ident=~sqrt(Idade+1))) plot(out$gam,ylab="Altura") # #Fazendo um alisamento no grafico de perfis plot(lowess(dados$Idade,Altura), col = 2, type="l",ylab="Altura",xlab="Idade") lines(lowess(sqrt(dados$Idade+1),Altura), col = 4, type="l",ylab="Altura",xlab="Idade") # # Utilizando a transforma??o raiz + 1 # out<-lme(Altura~sqrt(Idade+1)+status+sexo+ status:sqrt(Idade+1) +sexo:sqrt(Idade+1),random=list(ident=~sqrt(Idade+1))) plot(out) summary(out) # #================================================================== # Fazendo o Ajuste com Duas Inclina??es #================================================================== S<-length(dados$Idade) Age<-rep(0,S) for(i in 1:S) { if(dados$Idade[i]>5) {Age[i]<-dados$Idade[i]-7} } # dados1<-data.frame(Idade,Age,sexo,status) attach(dados1) Idade<-Idade - 7 # out<-lme(Altura~Idade+Age+status+sexo+status:Idade +sexo:Idade+status:Age+sexo:Age,random=list(ident=~Idade+Age)) plot(out) #========================================================================== # GEEPACK #============================================================================ # # Modelos Marginais sao ajustados utilizando a biblioteca geepack # library(geepack) attach(dados) Idade<-Idade-mean(Idade) # # Regressao Linear.. # out0<-lm(Altura ~ Idade+status+sexo+I(Idade^2)+Idade:status+I(Idade^2):status+ Idade:sexo+I(Idade^2):sexo) summary(out0) # Autoregressiva (dois componentes) # NAO PODE POIS OS TEMPOS NAO SAO IGUALMENTE ESPACADOS # out<-geeglm(Altura ~ Idade*status+I(Idade^2)*status+Idade*sexo+ # I(Idade^2)*sexo, id = ident, corstr = "AR1") # summary(out) # Esferica (dois componentes) # ESTRANHA POIS TEMOS 470 TEMPOS DISTINTOS out1<-geeglm(Altura ~ Idade*status+I(Idade^2)*status+Idade*sexo+ I(Idade^2)*sexo, id = ident, corstr = "exchangeable", data=dados) summary(out1) # Independencia (um componente) out2<-geeglm(Altura ~ Idade*status+I(Idade^2)*status+Idade*sexo+ I(Idade^2)*sexo, id = ident, corstr = "independence", data=dados) summary(out2) # #========================================================================== # Modelando a Matriz de Covari?ncia - GEE #============================================================================ # # Modelos Marginais sao ajustados utilizando a biblioteca gee # library(gee) attach(dados) Idade<-Idade-mean(Idade) # Autoregressiva (dois parametros) # Nao pode pois nao ? igualmente espacado # out<-gee(Altura ~ Idade*status+I(Idade^2)*status+Idade*sexo+ # I(Idade^2)*sexo, id = ident, corstr = "AR-M",Mv=1) # summary(out) # Esferica (dois Parametros) out1<-gee(Altura ~ Idade*status+I(Idade^2)*status+Idade*sexo+ I(Idade^2)*sexo, id = ident, data=dados, corstr = "exchangeable") summary(out1) # Independencia (um parametro) out<-gee(Altura ~ Idade*status+I(Idade^2)*status+Idade*sexo+ I(Idade^2)*sexo, id = ident, corstr = "independence",data=dados) summary(out) # # #============================================================================ # Modelando a Matriz de Covari?ncia - GLS (library -nlme) #========================================================================== # # Modelos Marginais ajustados pela biblioteca nlme # library(nlme) attach(dados) #Idade<-Idade-mean(Idade) # Simetria Composta - Mesma Vari?ncia (dois componentes de vari?ncia) out<-gls(Altura ~ Idade*status+I(Idade^2)*status+Idade*sexo+I(Idade^2)*sexo, correlation = corCompSymm(form = ~1 |ident)) # weights = varIdent(form = ~1| Idade)) # N?o ? possivel vari?ncias diferentes. summary(out) # AR(1) (5 componentes de Vari?ncia) out2<-gls(Altura ~ Idade*status+I(Idade^2)*status+Idade*sexo+I(Idade^2)*sexo,correlation=corCAR1(form=~Idade|ident)) out2 summary(out2) anova(out,out2) # #=================================================================================