[R-br] Método monte carlo

Fernando Antonio de souza nandodesouza em gmail.com
Terça Junho 2 15:27:39 BRT 2015


Walmes fiz algumas pequenas alterações no código que me enviou para
conseguir obter uma saída para o modelo lme. Desculpe ficar perguntando,
mas não conheco bem essa técnica.

A questão é: Baseada na saída abaixo, como eu a torno  interpretável? Faço
a média desses valores? Estou pensado em como avaliar se o modelo ajustou.
É possivel com essa metodologia obter estimativas de MSE, MAPE, CCC, etc.
Estou confuso sobre isso.

abraços


#Dados para o CMR abaixo:

modelo3<-lme(fixed=
CMS~PM+GPD+FDN+I(FDN^2),data=da,random=~1|Estudo,na.action=na.omit,method="REML")

L <- vector(mode="list", length=nrow(da))
for(i in 1:nrow(da)){
    L[[i]] <-lme(fixed=
CMS~PM+GPD+FDN+I(FDN^2),random=~1|Estudo,na.action=na.omit,method="REML",
data=da[-i,])
}

## Medidas de deviance leave-one-out.

sapply(L, FUN=function(x){-2*logLik(x)}) #Deviance

## log-veross leave-one-out.
sapply(L, logLik)

da<-structure(list(Estudo = structure(c(1L, 1L, 1L, 1L, 1L, 2L, 2L,
2L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 6L, 6L, 6L,
6L, 6L, 6L, 6L, 6L, 7L, 7L, 7L, 7L, 8L, 8L, 8L, 8L, 9L, 9L, 9L,
10L, 10L, 10L, 10L, 11L, 11L, 11L, 11L, 12L, 12L, 12L, 12L, 13L,
13L, 13L, 13L, 13L, 13L, 14L, 14L, 14L), .Label = c("1", "2",
"3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", "14",
"15", "16", "17", "18", "19", "20", "21", "22", "23", "24", "25",
"26", "27", "28", "29", "30", "31", "32", "33", "34", "35", "36",
"37", "38", "39", "40", "41", "42", "43", "44"), class = "factor"),
    NANIMAL = c(20L, 20L, 20L, 20L, 20L, 18L, 18L, 18L, 20L,
    20L, 20L, 20L, 20L, 20L, 20L, 20L, 20L, 18L, 18L, 18L, 16L,
    16L, 16L, 16L, 16L, 16L, 16L, 16L, 24L, 24L, 24L, 24L, 20L,
    20L, 20L, 20L, 18L, 18L, 18L, 24L, 24L, 24L, 24L, 24L, 24L,
    24L, 24L, 20L, 20L, 20L, 20L, 18L, 18L, 18L, 18L, 18L, 18L,
    10L, 10L, 10L), GENOTIPO = structure(c(2L, 2L, 2L, 2L, 2L,
    1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L,
    2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), .Label = c("1",
    "2"), class = "factor"), SEXO = structure(c(1L, 1L, 1L, 1L,
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
    1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L,
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), .Label = c("1",
    "2"), class = "factor"), VOL = structure(c(1L, 1L, 1L, 1L,
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L,
    2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L,
    2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L,
    2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L), .Label = c("1",
    "2"), class = "factor"), PVI = c(25, 25, 25, 25, 25, 20.17,
    20.1, 20.47, 21.3, 22, 22.9, 22.9, 18.25, 18.25, 22.47, 19.05,
    16.63, 22.62, 22.8, 22.95, 26.5, 26.5, 26.5, 26.5, 26.5,
    26.5, 26.5, 26.5, 19.1, 19.16, 19.96, 19.96, 20, 20, 20,
    20, 17, 16.25, 15.58, 16.4, 16.1, 16.6, 16.5, 18.74, 18.05,
    18.82, 18.15, 35.47, 35.7, 35.53, 35.66, 22.6, 22.6, 22.6,
    22.6, 22.6, 22.6, 15, 25, 35), PVF = c(36.8, 35.6, 33.88,
    35.8, 32.48, 31.87, 31.1, 31.9, 27.1, 27.6, 28.7, 31, 30.46,
    30.46, 32.15, 24.28, 19.62, 37.67, 39.48, 39.43, 30.15, 31.04,
    31.88, 31.92, 30.11, 29.44, 29.48, 30.87, 33.5, 32.66, 30.41,
    32.16, 25.23, 26.9, 28.65, 29.23, 24.66, 25.08, 27.75, 28.7,
    29.6, 31.1, 30.9, 31.7, 32.42, 34.97, 32.65, 44.21, 46.18,
    44.62, 45.55, 30, 27.73, 29.06, 32.7, 30.1, 28.8, 25, 35,
    45), PV = c(30.9, 30.3, 29.44, 30.4, 28.74, 26.02, 25.6,
    26.19, 24.2, 24.8, 25.8, 26.95, 24.36, 24.36, 27.31, 21.67,
    18.13, 30.15, 31.14, 31.19, 28.33, 28.77, 29.19, 29.21, 28.31,
    27.97, 27.99, 28.69, 26.3, 25.91, 25.19, 26.06, 22.62, 23.45,
    24.33, 24.62, 20.83, 20.67, 21.67, 22.55, 22.85, 23.85, 23.7,
    25.22, 25.24, 26.9, 25.4, 39.84, 40.94, 40.08, 40.61, 26.3,
    25.17, 25.83, 27.65, 26.35, 25.7, 20, 30, 40), PM = c(13.11,
    12.91, 12.64, 12.95, 12.41, 11.52, 11.38, 11.58, 10.91, 11.11,
    11.45, 11.83, 10.96, 10.96, 11.95, 10.04, 8.78, 12.87, 13.18,
    13.2, 12.28, 12.42, 12.56, 12.56, 12.27, 12.16, 12.17, 12.39,
    11.61, 11.48, 11.24, 11.53, 10.37, 10.66, 10.95, 11.05, 9.75,
    9.69, 10.04, 10.35, 10.45, 10.79, 10.74, 11.25, 11.26, 11.81,
    11.31, 15.86, 16.18, 15.93, 16.09, 11.61, 11.24, 11.46, 12.06,
    11.63, 11.41, 9.46, 12.82, 15.91), GPD = c(295, 265, 222,
    270, 187, 123, 137, 191, 86.5, 100.4, 92.1, 140.5, 195.83,
    195.83, 161.25, 87.08, 49.75, 268.75, 297.85, 294.28, 87,
    108, 128, 129, 86, 70, 71, 104, 206, 186, 149, 174, 93.4,
    123.2, 154.4, 164.8, 88.12, 101.53, 139.84, 195, 214.5, 229.9,
    228, 245, 271, 305, 274, 202, 243, 211, 212, 117.5, 81.5,
    102.6, 160.3, 119, 98.4, 205, 261, 183), D.EXP = c(40, 40,
    40, 40, 40, 98, 97, 77, 77, 77, 77, 77, 60, 60, 60, 60, 60,
    56, 56, 56, 42, 42, 42, 42, 42, 42, 42, 42, 70, 70, 70, 70,
    56, 56, 56, 56, 87, 87, 87, 63, 63, 63, 63, 53, 53, 53, 53,
    40, 40, 40, 40, 63, 63, 63, 63, 63, 63, 52, 55, 65), CMS = c(1490,
    1540, 1390, 1580, 1440, 900, 870, 870, 422, 501, 525, 533,
    371, 1377, 975.7, 596.7, 448.9, 1260, 1220, 1170, 889.56,
    842.96, 945.76, 1051.56, 888.93, 819.52, 906.88, 1044.11,
    946.8, 1236, 1126, 1192, 648.6, 764.4, 960, 941, 845.28,
    845.61, 780.59, 820.8, 830.6, 870, 887.8, 1098, 1159, 1289,
    1173, 1200, 1340, 1300, 1230, 1001, 987, 861, 896, 1037,
    785, 707.6, 977.3, 1032.9), FDN = c(41.52, 42.99, 53.49,
    51.2, 56.45, 58.3, 44.23, 29.96, 59.18, 55.36, 51.55, 47.73,
    43.91, 52.3, 52.8, 47.7, 44.9, 34.5, 34.5, 34.5, 64.21, 60.77,
    61.62, 59.32, 64.21, 60.77, 61.62, 59.32, 36.4, 41.6, 49.6,
    49.6, 54.7, 55.9, 54, 56.6, 35.46, 34.87, 34.26, 44.2, 43.44,
    42.4, 41.7, 43.36, 39.62, 35.88, 32.15, 31.46, 35.55, 39.54,
    44.07, 12.37, 25.47, 39.16, 12.37, 25.47, 39.16, 21.41, 21.41,
    21.41), CMSPV = c(48.22, 50.83, 47.21, 51.97, 50.1, 34.59,
    33.98, 33.23, 31.52, 35.31, 37.15, 37.12, 15.23, 56.54, 35.73,
    27.54, 24.77, 41.8, 39.18, 37.51, 31.4, 29.3, 32.4, 36, 31.4,
    29.3, 32.4, 36.39, 36, 47.7, 44.7, 45.74, 28.67, 32.6, 39.46,
    38.22, 40.58, 40.91, 36.02, 36.4, 36.35, 36.48, 37.46, 43.54,
    45.92, 47.92, 46.18, 30.12, 32.73, 32.44, 30.29, 38.06, 39.21,
    33.33, 32.41, 39.35, 30.54, 35.38, 32.58, 25.82), CMSPM = c(113.69,
    119.24, 109.98, 122.04, 116.01, 78.12, 76.44, 75.16, 64.76,
    72.67, 78.02, 78.42, 60.33, 125.6, 81.67, 59.42, 51.1, 97.94,
    92.55, 88.65, 72.44, 67.87, 75.3, 83.72, 72.45, 67.39, 74.52,
    84.27, 105.57, 107.63, 100.16, 103.35, 62.54, 71.73, 87.65,
    85.15, 87.68, 88.32, 79.27, 79.32, 79.47, 80.61, 82.65, 97.56,
    102.94, 109.14, 103.67, 80.17, 88.22, 86.12, 81.77, 78.46,
    83.78, 75.14, 74.31, 89.17, 68.77, 74.82, 76.24, 64.94)), .Names =
c("Estudo",
"NANIMAL", "GENOTIPO", "SEXO", "VOL", "PVI", "PVF", "PV", "PM",
"GPD", "D.EXP", "CMS", "FDN", "CMSPV", "CMSPM"), row.names = c(NA,
60L), class = "data.frame")


Em 1 de junho de 2015 23:18, Fernando Antonio de souza <
nandodesouza em gmail.com> escreveu:

> Obrigado, Walmes,
>
> Irei testar aqui e retornarei a lista. Irei estudar mais sobre o tema.
> Hoje li alguns tutoriais sobre o assunto e realmente tenho a impressão de
> ser muito simples. Mas não tinha ideía o quão simples era para
> implementá-la.
>
> Abçs
>
> Em 1 de junho de 2015 21:23, walmes . <walmeszeviani em gmail.com> escreveu:
>
>> Eu posso estar enganado então me corrijam em caso afirmativo. O
>> leave-ONE-out é mais simples do que imagina. Grosseiramente falando, crie
>> um laço for para deixar um caso de fora e ajuste o modelo. Depois estude
>> uma medida de ajuste. No caso de lm, não requer de fato fazer a exaustiva
>> tarefa de ajustar n modelos (deixando um caso de fora), pois se chega as
>> medidas leave-one-out por projeções matriciais, etc. Mas num caso mais
>> geral é algo como:
>>
>> da ## seu data.frame
>> L <- vector(mode="list", length=nrow(da))
>> for(i in 1:nrow(da)){
>>     L[[i]] <- sua_funcao_R(..., data=da[-i,])
>> }
>>
>> ## Medidas de deviance leave-one-out.
>> sapply(L, deviance)
>>
>> ## log-veross leave-one-out.
>> sapply(L, logLik)
>>
>> Código não testado.
>>
>> À disposição.
>> Walmes.
>>
>>
>> _______________________________________________
>> R-br mailing list
>> R-br em listas.c3sl.ufpr.br
>> https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br
>> Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forneça
>> código mínimo reproduzível.
>>
>
>
>
> --
> =======================================================================
> Fernando Souza
> Zootecnista, DSc. Produção Animal
> e-mail:nandodesouza em gmail.com
> https://producaoanimalcomr.wordpress.com/
> ========================================================================
>



-- 
=======================================================================
Fernando Souza
Zootecnista, DSc. Produção Animal
e-mail:nandodesouza em gmail.com
https://producaoanimalcomr.wordpress.com/
========================================================================
-------------- Próxima Parte ----------
Um anexo em HTML foi limpo...
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20150602/20f7a762/attachment.html>


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