Re: [R-br] Função estimate para o R
Caro Walmes, Os comando utilizados até agora foram os seguintes
RegFixo<-glm(Y~-1+as.factor(estudo)/X,family=gaussian(),data=Pierre3)#considera os efietos fixos porém considere o efeito do estudo library("car") Anova(RegFixo,type="III",test.statistic="F") #produz a soma de quadrados tipo III anova(RegFixo,test="F")#produz a soma de quadrados tipo I summary(RegFixo)
Acontece que necessito estimar o modelo baseado nos erros padrões estimados pela função glm e assim obter o intercepto total e a inclinação total do modelo. E não sei qual função ou como fazer isto.
Eu anexei o artigo para q você olhe caso seja necessário. trata-se de um passo a passo apresentado pelo autor para a realização de meta-analises.Ele fornece o banco de dados utilizado bem como apresenta os comando utilizados no SAS por ele. Eu estou repetindo os passo a passos para compreender o pensamento por trás deste tipo de análise, porém quero realizá-los no R!
Abaixo segue o banco de dados (também coloquei em anexo em formato csv)
Base de dados
obs estudo X Y
1 1 1.16391 -2.92296 2 1 1.94273 -1.03365 3 1 2.50229 -1.00383 4 1 3.98627 1.33325 5 1 4.36177 1.13752 6 1 4.98954 1.59278 7 2 2.1091 -1.33751 8 2 2.28976 -1.14376 9 2 2.45645 -1.09304 10 2 4.95084 1.4274 11 2 5.92572 2.99572 12 3 2.09605 -0.70073 13 3 2.30059 -1.11007 14 3 3.90388 0.17334 15 3 5.15574 2.38882 16 3 5.57553 1.78167 17 3 5.60523 3.09353 18 4 3.1955 0.18097 19 4 4.32691 0.7431 20 4 5.04092 1.947 21 4 5.17521 1.07029 22 4 5.57254 2.71423 23 4 5.59532 2.43982 24 5 2.31878 0.26487 25 5 4.06027 1.99979 26 6 2.48122 1.31151 27 6 2.73224 1.32218 28 6 3.66628 2.52286 29 6 4.15568 2.83423 30 6 4.57694 2.99086 31 6 4.86356 3.90049 32 7 2.767 1.18358 33 7 2.79859 1.70206 34 7 2.97074 1.84639 35 7 3.85289 3.13155 36 7 6.24981 5.27271 37 7 6.80389 5.85844 38 8 2.74107 -0.26691 39 8 3.72445 1.32567 40 8 4.51721 1.04086 41 8 5.18536 2.40105 42 9 3.24046 3.12334 43 9 3.2864 3.52558 44 9 4.52824 4.29043 45 9 5.10468 4.04917 46 9 5.40578 5.55135 47 9 7.25461 7.2311 48 10 3.09531 3.31464 49 10 3.13951 2.5785 50 10 3.86928 4.0819 51 10 6.40856 6.6731 52 10 7.35976 6.8001 53 10 7.95716 8.2098 54 11 3.24974 2.8825 55 11 3.25369 4.8511 56 11 4.4005 4.3538 57 11 4.74897 6.076 58 11 4.96576 5.9062 59 12 3.00403 2.5733 60 12 3.1542 3.9324 61 12 4.95652 5.8456 62 12 5.35746 5.2066 63 12 5.62651 6.3312 64 12 6.53297 6.6019 65 12 6.6958 6.9498 66 12 7.83844 8.1031 67 13 5.44713 6.4838 68 13 5.82682 6.2341 69 13 6.76986 7.8937 70 13 6.87949 6.9477 71 13 8.47819 8.1336 72 14 4.31493 5.228 73 14 5.70501 6.0941 74 14 6.70113 7.6858 75 15 4.11223 5.2129 76 15 4.55329 5.2169 77 15 4.75641 6.1755 78 15 5.0742 5.7492 79 15 5.90125 6.5407 80 15 6.84886 6.7004 81 15 8.87813 10.2425 82 16 4.25114 5.3359 83 16 6.50762 9.1897 84 16 6.52857 8.4259 85 16 6.67694 8.2782 86 16 6.79539 8.3206 87 16 8.62833 10.4385 88 17 6.47237 7.6132 89 17 6.73826 7.7381 90 17 7.57947 8.5328 91 17 8.11267 8.9045 92 17 9.18169 9.7732 93 18 5.49319 6.6492 94 18 5.62711 6.4464 95 18 6.67182 8.5909 96 18 7.4494 9.2077 97 18 8.13165 9.365 98 18 8.31896 9.5025 99 18 8.74606 10.3598 100 18 8.92038 9.4361 101 18 9.43227 11.0051 102 19 6.06067 7.9306 103 19 7.98874 10.1887 104 19 8.35784 10.9856 105 19 9.52789 11.507 106 19 9.68213 12.6051 107 20 7.15668 9.9928 108 20 9.61256 12.2149
grato pela atenção
Fernando, Para isso você precisa saber qual a definição dele de intercepto-inclinação geral. Penso que seja aquela obtida usando a restrição soma zero nas estimativas, mas pode ser outra. Para conseguir essas estimativas você nem precisa de pacote extra, é sou definir corretamente o tipo de contraste (ou restrição) que você quer. Evite enviar dados em anexo ou colados na mensagem. Isso não deixa o seu código reproduzível com um copia e cola. Eu teria que fazer download, blá blá, isso enrola e eu prefiro gerar dados do que ter que todo esse trabalhão. Para facilitar a vida de quem tá disposto a te ajudar, passe o dput() do seus dados, que assim com um copia e cola reproduzímos o seu código. Mostro isso no final. da <- expand.grid(estudo=gl(5,1), x=1:10) da$y <- round(rnorm(nrow(da),0,1),2) m0 <- lm(y~estudo*x, data=da) coef(m0) m1 <- lm(y~-1+estudo/x, data=da) # mesmo espaço coluna do anterior coef(m1) # são os interceptos e inclinações para cada estudo m2 <- lm(y~estudo*x, data=da, contrast=list(estudo=contr.sum)) coef(m2) coef(m2)[c(1,6)] # intercepto e inclinação geral coef(m2)[2:5] # interceptos desvios do geral para estudos de 1 a 4 -sum(coef(m2)[2:5]) # estudo 5 coef(m2)[7:10] # inclinações desvios da geral para estudos de 1 a 4 -sum(coef(m2)[7:10]) # estudo 5 summary(m2) # o resultado dessa função desse ver enviado para um # código reproduzível com copia e cola dput(da) À disposição. Walmes. ========================================================================== Walmes Marques Zeviani LEG (Laboratório de Estatística e Geoinformação, 25.450418 S, 49.231759 W) Departamento de Estatística - Universidade Federal do Paraná fone: (+55) 41 3361 3573 VoIP: (3361 3600) 1053 1173 e-mail: walmes@ufpr.br twitter: @walmeszeviani homepage: http://www.leg.ufpr.br/~walmes linux user number: 531218 ==========================================================================
Caro walmes, Mais uma vez meu muito obrigado realmente era o que eu desejava fazer. até a próxima Em 23 de setembro de 2011 11:03, Walmes Zeviani <walmeszeviani@gmail.com>escreveu:
Fernando,
Para isso você precisa saber qual a definição dele de intercepto-inclinação geral. Penso que seja aquela obtida usando a restrição soma zero nas estimativas, mas pode ser outra. Para conseguir essas estimativas você nem precisa de pacote extra, é sou definir corretamente o tipo de contraste (ou restrição) que você quer. Evite enviar dados em anexo ou colados na mensagem. Isso não deixa o seu código reproduzível com um copia e cola. Eu teria que fazer download, blá blá, isso enrola e eu prefiro gerar dados do que ter que todo esse trabalhão. Para facilitar a vida de quem tá disposto a te ajudar, passe o dput() do seus dados, que assim com um copia e cola reproduzímos o seu código. Mostro isso no final.
da <- expand.grid(estudo=gl(5,1), x=1:10) da$y <- round(rnorm(nrow(da),0,1),2)
m0 <- lm(y~estudo*x, data=da) coef(m0) m1 <- lm(y~-1+estudo/x, data=da) # mesmo espaço coluna do anterior coef(m1) # são os interceptos e inclinações para cada estudo
m2 <- lm(y~estudo*x, data=da, contrast=list(estudo=contr.sum)) coef(m2) coef(m2)[c(1,6)] # intercepto e inclinação geral coef(m2)[2:5] # interceptos desvios do geral para estudos de 1 a 4 -sum(coef(m2)[2:5]) # estudo 5 coef(m2)[7:10] # inclinações desvios da geral para estudos de 1 a 4 -sum(coef(m2)[7:10]) # estudo 5
summary(m2)
# o resultado dessa função desse ver enviado para um # código reproduzível com copia e cola dput(da)
À disposição. Walmes.
========================================================================== Walmes Marques Zeviani LEG (Laboratório de Estatística e Geoinformação, 25.450418 S, 49.231759 W) Departamento de Estatística - Universidade Federal do Paraná fone: (+55) 41 3361 3573 VoIP: (3361 3600) 1053 1173 e-mail: walmes@ufpr.br twitter: @walmeszeviani homepage: http://www.leg.ufpr.br/~walmes linux user number: 531218 ==========================================================================
_______________________________________________ R-br mailing list R-br@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.
participantes (2)
-
Fernando Antonio de souza -
Walmes Zeviani