[R-br] Modelo efeito misto - decomposição R2 marginal

Fernando Souza nandodesouza em gmail.com
Segunda Maio 26 23:34:32 BRT 2014 

Caros Amigos,
Fiz uma análise de regressão utilizando o modelo de efeito misto. Eu 
desejo fazer o coeficiente de determinação parcial conforme descrito no 
livro do Kutner, et al (Applied linear statistical Models, ) 
<http://www.amazon.com/Applied-Linear-Statistical-Models-Michael/dp/007310874X> 
5° edição,   paginas 268-269 para obter uma estimativa de quanto cada 
variável incluída no modelo contribuiu para o R2 geral ( ou seja a 
contribuição marginal de cada variável para a redução da variação quando 
todas as outras variáveis estão no modelo- No livro é chamado de 
coeficiente parcial de determinação). O R2 do modelo deve ser igual a 
soma do R2 parcial de cada variável individualment (R2modelo= 
R2variável1+R2variável2...). Fiz uma função para automatizar os passos 
descritos no livro, como vocês poderão ver a soma do R2 parcial de cada 
variável difere do R2 geral do modelo. Acredito que eu possa ter 
entendido algo errado, por isso gostaria de ajuda. Alguém da lista ja 
fez este tipo de avaliação?


DATA<- structure(list(DATA = structure(c(25L, 27L, 1L, 2L, 3L, 4L, 5L,
6L, 7L, 8L, 9L, 10L, 11L, 12L, 13L, 14L, 15L, 16L, 17L, 18L,
19L, 20L, 21L, 22L, 23L, 24L, 26L, 28L, 25L, 27L), .Label = c("01/09/13",
"02/09/13", "03/09/13", "04/09/13", "05/09/13", "07/09/13", "08/09/13",
"10/09/13", "11/09/13", "14/09/13", "15/09/13", "16/09/13", "17/09/13",
"18/09/13", "19/09/13", "20/09/13", "21/09/13", "22/09/13", "23/09/13",
"24/09/13", "25/09/13", "26/09/13", "27/09/13", "28/09/13", "29/08/13",
"29/09/13", "30/08/13", "30/09/13"), class = "factor"), GEST = c(63L,
64L, 66L, 67L, 68L, 69L, 70L, 72L, 73L, 75L, 76L, 79L, 80L, 81L,
82L, 83L, 84L, 85L, 86L, 87L, 88L, 89L, 90L, 91L, 92L, 93L, 94L,
95L, 63L, 64L), MANEJO = c(2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L), DIAS = c(1L, 2L, 4L, 5L, 6L, 7L, 8L, 10L,
11L, 13L, 14L, 17L, 18L, 19L, 20L, 21L, 22L, 23L, 24L, 25L, 26L,
27L, 28L, 29L, 30L, 31L, 32L, 33L, 1L, 2L), ANIMAL = structure(c(2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L), .Label = c("5",
"4"), class = c("ordered", "factor")), CONSUMO = c(1.43, 1.24,
2, 1.29, 0.94, 1.99, 1.48, 1.25, 1.8, 2.16, 1.36, 1.33, 2.25,
1.1, 1.78, 1.47, 2, 2.07, 1.35, 1.84, 1.86, 1.5, 1.31, 1.45,
1.46, 2.51, 2.77, 1.89, 1.39, 1.68), PV = c(40, 40.06, 40.18,
40.24, 40.3, 40.36, 40.42, 40.54, 40.6, 40.72, 40.78, 40.99,
41.04, 41.08, 41.13, 41.17, 41.22, 41.26, 41.31, 41.35, 41.4,
41.45, 41.49, 41.54, 41.58, 41.63, 41.67, 41.9, 32.3, 32.47),
     TEMPMIN = c(11.85, 11.85, 17.9, 18.45, 18.6, 19.6, 19.2,
     17.8, 18.7, 15.5, 15.55, 18.7, 18.15, 17.55, 17.85, 18.75,
     20.95, 19.3, 19.35, 19.4, 18.8, 17.65, 17.75, 16.45, 16.95,
     18.4, 21.1, 18.2, 11.85, 11.85), TEMPMAX = c(28.05, 28.05,
     30, 30.9, 30.8, 30.35, 31.2, 27.9, 31.3, 28.6, 28.25, 31.95,
     32.75, 30.7, 31.65, 32.85, 32.65, 33.2, 33.05, 24.7, 24.35,
     24.95, 28, 28.15, 30.9, 33.1, 34.95, 26.65, 28.05, 28.05),
     URAMIN = c(38, 38, 40, 41.5, 43, 42.5, 41, 61.5, 44.5, 41.5,
     41, 18, 19.5, 26.5, 19.5, 19.5, 20, 25.5, 29, 70, 71.5, 42.5,
     50, 44, 31, 23, 15.65, 46.5, 38, 38), URAMAX = c(89, 89,
     94, 91, 91, 96.5, 94.5, 99, 95.5, 91, 93, 74.5, 74.5, 75.5,
     77, 76, 80, 88, 92, 94, 93.5, 94, 91.5, 85.5, 83.5, 79, 51.75,
     77.5, 89, 89), TEMPMEC = c(19.95, 19.95, 23.95, 24.68, 24.7,
     24.98, 25.2, 22.85, 25, 22.05, 21.9, 25.33, 25.45, 24.13,
     24.75, 25.8, 26.8, 26.25, 26.2, 22.05, 21.58, 21.3, 22.88,
     22.3, 23.93, 25.75, 28.03, 22.43, 19.95, 19.95), URAMED = c(63.5,
     63.5, 66.25, 69.5, 67.75, 80.25, 70, 67, 47.75, 47, 48.25,
     56.75, 60.5, 82.5, 68.25, 70.75, 64.75, 51, 63.5, 62, 56.75,
     84.25, 84.25, 64.25, 64.5, 57.75, 72.25, 74.75, 63.5, 63.5
     ), ITH = c(4, 4, 151, 220, 212, 304, 267, 82, 133, 20, 16,
     201, 231, 246, 218, 317, 343, 230, 302, 40, 21, 49, 147,
     51, 143, 238, 409, 82, 4, 4), IMS = c(0.9, 0.62, 0.69, 1.62,
     0.81, 1.16, 1.45, 0.87, 0.86, 0.86, 1.04, 1.68, 0.68, 0.56,
     1.47, 1.04, 1.14, 0.66, 0.97, 1.35, 1.14, 1.05, 1.52, 0.77,
     0.95, 0.84, 1.31, 1.52, 1, 1.37), cook = c(0, 0, 0, 0, 0,
     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
     0, 0, 0, 0, 0, 0)), .Names = c("DATA", "GEST", "MANEJO",
"DIAS", "ANIMAL", "CONSUMO", "PV", "TEMPMIN", "TEMPMAX", "URAMIN",
"URAMAX", "TEMPMEC", "URAMED", "ITH", "IMS", "cook"), row.names = c("178",
"179", "181", "182", "183", "184", "185", "187", "188", "190",
"191", "194", "195", "196", "197", "198", "199", "200", "201",
"202", "203", "204", "205", "206", "207", "208", "209", "210",
"211", "212"), class = c("nffGroupedData", "nfGroupedData", "groupedData",
"data.frame"), formula = CONSUMO ~ I(PV^0.75) + IMS + GEST +
     TEMPMAX + TEMPMIN + URAMIN + URAMAX + as.numeric(ITH) | ANIMAL, FUN 
= function (x)
max(x, na.rm = TRUE), order.groups = TRUE)

modelocompleto<-lme(fixed=CONSUMO~GEST+as.numeric(ITH),data=DATA,correlation=corAR1(form=~GEST|ANIMAL),na.action=na.omit,weights=varIdent(form=~1|ANIMAL),random=~1|ANIMAL,method="REML")
modeloreduzidoGEST<-update(modelocompleto,fixed=CONSUMO~as.numeric(ITH),random=~1|ANIMAL,method="REML",control=ctrl)
modeloreduzidoITH<-update(modelocompleto,fixed=CONSUMO~GEST,random=~1|ANIMAL,method="REML",control=ctrl)

R2partial<-function(modelocompleto,modeloreduzidoGEST,modeloreduzidoITH){
SSRGEST<-(sum(resid(modeloreduzidoGEST)^2) - sum(resid(modelocompleto)^2))
SSRITH<-(sum(resid(modeloreduzidoITH)^2) - sum(resid(modelocompleto)^2))
SSEGEST<-sum(resid(modeloreduzidoGEST)^2)
SSEITH<-sum(resid(modeloreduzidoITH)^2)
R2partialGEST<-SSRGEST/SSEGEST
R2partialITH<-SSRITH/SSEITH
resultado1<-paste("R2parcial-GEST",R2partialGEST,sep=":")
resultado2<-paste("R2parcial-ITH",R2partialITH,sep=":")
R2total<-paste("R2total",sum(R2partialGEST,R2partialITH),sep=":")
return(c(resultado1,resultado2,R2total))
}
R2partial(modelocompleto,modeloreduzidoGEST,modeloreduzidoITH) # R2 
parcial para cada uma das vaiáveis

install.packages("MuMIn")
library(MuMIn)
r.squaredGLMM(modeloreduzidoGEST) #R2 Geral

-------------- Próxima Parte ----------
Um anexo em HTML foi limpo...
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20140526/cad66115/attachment.html>

Mais detalhes sobre a lista de discussão R-br