[R-br] Partição de estimativas em GLM
Fernando Souza
nandodesouza em gmail.com
Sábado Maio 10 18:08:06 BRT 2014
Caros Amigos,
Estou refazendo os comandos apresentado em um artigo de meta-análise.
Acontece que o autor está utilizando o SAS e gostaria de refazer o
comando utilizando o R:
O código utilizado pelo autor no SAS foi:
PROC GLM DATA=Dataregs;
CLASS Study;
MODEL Y = X Study X*Study/SOLUTION;
RANDOM Study;
LSMEANS Study/at X=0 STDERR;
ESTIMATE 'Overall Intercept' INTERCEPT 20
Study 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1/DIVISOR=20;
ESTIMATE 'Overall Slope' X 20
X*Study 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1/DIVISOR=20;
ESTIMATE 'Slope Study 1' X 1
X*Study 1;
ESTIMATE 'Slope Study 2' X 1
X*Study 0 1;
ESTIMATE 'Slope Study 3' X 1
X*Study 0 0 1;
.
.
.
ESTIMATE 'Slope Study 18' X 1
X*Study 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1;
ESTIMATE 'Slope Study 19' X 1
X*Study 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1;
ESTIMATE 'Slope Study 20' X 1
X*Study 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1;
RUN;
O comando ESTIMATE, segundo o autor estima os coeficientes do intercepto
e da inclinação através de todos os estudos (linha 1 e 2) e os demais
calcula estimativas de inclinação para cada estudo.
Estou utilizando o seguinte código
RegFixo<-lm(Y~X*estudo,contrasts=list(estudo="contr.SAS"),data=DATA)
RegFixo<-lm(Y~X*estudo,contrasts=list(estudo="contr.sum"),data=DATA) #
Este comando é o mesmo que fazer o Estimate?
library(car)
Anova(RegFixo,type="III")
*Minha dúvida é como calcular os ESTIMATES no R conforme comando acima
do SAS?*
O dput() da base da dados que estou utilizando é fornecido abaixo.
Também segue o artigo para eventual dúvida sobre da saída obtida pelo
autor no SAS. A minha dúvida refere-se ao comando das páginas 745,746,747
DATA<-structure(list(OBS = 1:108, estudo = structure(c(1L, 1L, 1L,
1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L,
4L, 4L, 4L, 4L, 5L, 5L, 6L, 6L, 6L, 6L, 6L, 6L, 7L, 7L, 7L, 7L,
7L, 7L, 8L, 8L, 8L, 8L, 9L, 9L, 9L, 9L, 9L, 9L, 10L, 10L, 10L,
10L, 10L, 10L, 11L, 11L, 11L, 11L, 11L, 12L, 12L, 12L, 12L, 12L,
12L, 12L, 12L, 13L, 13L, 13L, 13L, 13L, 14L, 14L, 14L, 15L, 15L,
15L, 15L, 15L, 15L, 15L, 16L, 16L, 16L, 16L, 16L, 16L, 17L, 17L,
17L, 17L, 17L, 18L, 18L, 18L, 18L, 18L, 18L, 18L, 18L, 18L, 19L,
19L, 19L, 19L, 19L, 20L, 20L), .Label = c("1", "2", "3", "4",
"5", "6", "7", "8", "9", "10", "11", "12", "13", "14", "15",
"16", "17", "18", "19", "20"), class = "factor"), X = c(1.16,
1.94, 2.5, 3.99, 4.36, 4.99, 2.11, 2.29, 2.46, 4.95, 5.93, 2.1,
2.3, 3.9, 5.16, 5.58, 5.61, 3.2, 4.33, 5.04, 5.18, 5.57, 5.6,
2.32, 4.06, 2.48, 2.73, 3.67, 4.16, 4.58, 4.86, 2.77, 2.8, 2.97,
3.85, 6.25, 6.8, 2.74, 3.72, 4.52, 5.19, 3.24, 3.29, 4.53, 5.1,
5.41, 7.25, 3.1, 3.14, 3.87, 6.41, 7.36, 7.96, 3.25, 3.25, 4.4,
4.75, 4.97, 3, 3.15, 4.96, 5.36, 5.63, 6.53, 6.7, 7.84, 5.45,
5.83, 6.77, 6.88, 8.48, 4.31, 5.71, 6.7, 4.11, 4.55, 4.76, 5.07,
5.9, 6.85, 8.88, 4.25, 6.51, 6.53, 6.68, 6.8, 8.63, 6.47, 6.74,
7.58, 8.11, 9.18, 5.49, 5.63, 6.67, 7.45, 8.13, 8.32, 8.75, 8.92,
9.43, 6.06, 7.99, 8.36, 9.53, 9.68, 7.16, 9.61), Y = c(-2.92,
-1.03, -1, 1.33, 1.14, 1.59, -1.34, -1.14, -1.09, 1.43, 3, -0.7,
-1.11, 0.17, 2.39, 1.78, 3.09, 0.18, 0.74, 1.95, 1.07, 2.71,
2.44, 0.26, 2, 1.31, 1.32, 2.52, 2.83, 2.99, 3.9, 1.18, 1.7,
1.85, 3.13, 5.27, 5.86, -0.27, 1.33, 1.04, 2.4, 3.12, 3.53, 4.29,
4.05, 5.55, 7.23, 3.31, 2.58, 4.08, 6.67, 6.8, 8.21, 2.88, 4.85,
4.35, 6.08, 5.91, 2.57, 3.93, 5.85, 5.21, 6.33, 6.6, 6.95, 8.1,
6.48, 6.23, 7.89, 6.95, 8.13, 5.23, 6.09, 7.69, 5.21, 5.22, 6.18,
5.75, 6.54, 6.7, 10.24, 5.34, 9.19, 8.43, 8.28, 8.32, 10.44,
7.61, 7.74, 8.53, 8.9, 9.77, 6.65, 6.45, 8.59, 9.21, 9.37, 9.5,
10.36, 9.44, 11.01, 7.93, 10.19, 10.99, 11.51, 12.61, 9.99, 12.21
), YADJ = c(0.89, 0.64, 3, 4.36, 3.81, 4.36, 1.16, 1.58, 2.62,
4.43, 6.46, 1.47, 1.62, 2.92, 5.38, 5.41, 6.39, 3.77, 3.64, 4.88,
4.81, 5.12, 5.79, 1.65, 4.02, 2.03, 2.21, 3.65, 3.91, 4.51, 4.62,
2.13, 2.52, 2.72, 4.64, 5.31, 6.99, 2.35, 3.13, 4.78, 5.01, 3.31,
2.79, 4.08, 5.19, 5.6, 7.01, 2.15, 3.3, 4, 6.4, 7.9, 7.82, 2.97,
3.7, 3.95, 4.53, 4.47, 2.05, 3.67, 4.21, 6.01, 5.75, 5.88, 7.1,
7.87, 4.84, 5.41, 7.26, 7.11, 8.93, 4.04, 5.45, 7.07, 3.39, 5.02,
4.42, 5.19, 5.31, 7.49, 8.98, 4.3, 7.1, 6.06, 6.79, 6.9, 8.49,
6.39, 6.98, 8.21, 8.44, 8.84, 5.16, 5.32, 6.93, 8.7, 8.22, 8.63,
8.76, 9.27, 9.14, 5.89, 7.56, 9.06, 10.09, 10.1, 7.61, 9.78)), .Names =
c("OBS",
"estudo", "X", "Y", "YADJ"), row.names = c(NA, -108L), class = "data.frame")
Eu inseri um link para 1 arquivo nesta mensagem:
stpierre.pdf <https://www.box.com/shared/ym97esmv6biljsb725j6>(387
KB)Box
<https://www.box.com/thunderbird>https://www.box.com/shared/ym97esmv6biljsb725j6
O Mozilla Thunderbird <http://www.getthunderbird.com> torna fácil
compartilhar grandes arquivos pelo e-mail.
-------------- Próxima Parte ----------
Um anexo em HTML foi limpo...
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20140510/128cbea2/attachment.html>
-------------- Próxima Parte ----------
Um anexo não-texto foi limpo...
Nome: não disponível
Tipo: image/png
Tamanho: 641 bytes
Descrição: não disponível
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20140510/128cbea2/attachment.png>
-------------- Próxima Parte ----------
Um anexo não-texto foi limpo...
Nome: não disponível
Tipo: image/png
Tamanho: 766 bytes
Descrição: não disponível
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20140510/128cbea2/attachment-0001.png>
Mais detalhes sobre a lista de discussão R-br