Partição de estimativas em GLM
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.
Basicamente não existe segredo para reproduzir a mesma análise. O modo mais direto é declarar o modelo para estimar os parâmetros separados para cada nível do fator, depois montar as matrizes correspondentes às funções lineares e submetê-las à glht() ou qualquer outra função que permita inferências para funções lineares de parâmetros. O pacote car, contrast e gmodels têm funções equivalentes à glht(), minha preferida. ##----------------------------------------------------------------------------- m0 <- lm(Y~0+estudo/X, DATA) summary(m0) m0$assign tapply(coef(m0), m0$assign, mean) ## Matriz de pesos. m <- rbind(rep(1, nlevels(DATA$estudo))/nlevels(DATA$estudo)) require(multcomp) ## Intercepto médio. X <- cbind(m, 0*m) summary(glht(m0, linfct=X), test=adjusted(type="none")) ## Inclinação média. X <- cbind(0*m, m) summary(glht(m0, linfct=X), test=adjusted(type="none")) ## As inclinações para cada nível saem no próprio summary() do modelo ## quando se declara o modelo na forma ~0+fator/numérica. Para ## conhecimento, pode ser feito assim. M <- diag(nlevels(DATA$estudo)) X <- cbind(M, 0*M) summary(glht(m0, linfct=X), test=adjusted(type="fdr")) X <- cbind(0*M, M) summary(glht(m0, linfct=X), test=adjusted(type="fdr")) ##----------------------------------------------------------------------------- À disposição. Walmes.
Obrigado Walmes, era exatamente o que eu queria. Eu tinha em mente que deveria utilizar a função glht mas não sabia como. A questão é criar as matrizes de comparação. Preciso aprender mais sobre isso. Você teria algum material sobre a construção dessas matrizes ou recomenda algum texto? Obrigado pela atenção. Em Sáb 10 Mai 2014 22:29:10 BRT, walmes . escreveu:
Basicamente não existe segredo para reproduzir a mesma análise. O modo mais direto é declarar o modelo para estimar os parâmetros separados para cada nível do fator, depois montar as matrizes correspondentes às funções lineares e submetê-las à glht() ou qualquer outra função que permita inferências para funções lineares de parâmetros. O pacote car, contrast e gmodels têm funções equivalentes à glht(), minha preferida.
##-----------------------------------------------------------------------------
m0 <- lm(Y~0+estudo/X, DATA) summary(m0)
m0$assign tapply(coef(m0), m0$assign, mean)
## Matriz de pesos. m <- rbind(rep(1, nlevels(DATA$estudo))/nlevels(DATA$estudo))
require(multcomp)
## Intercepto médio. X <- cbind(m, 0*m) summary(glht(m0, linfct=X), test=adjusted(type="none"))
## Inclinação média. X <- cbind(0*m, m) summary(glht(m0, linfct=X), test=adjusted(type="none"))
## As inclinações para cada nível saem no próprio summary() do modelo ## quando se declara o modelo na forma ~0+fator/numérica. Para ## conhecimento, pode ser feito assim.
M <- diag(nlevels(DATA$estudo)) X <- cbind(M, 0*M) summary(glht(m0, linfct=X), test=adjusted(type="fdr"))
X <- cbind(0*M, M) summary(glht(m0, linfct=X), test=adjusted(type="fdr"))
##-----------------------------------------------------------------------------
À disposição. Walmes.
_______________________________________________ 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.
Oi Walmes aproveitando a questão, mais uma dúvida. É possível fazer um teste dos parâmetros do modelo de efeito aleatório utilizando a função lme? Em 10-05-2014 22:29, walmes . escreveu:
Basicamente não existe segredo para reproduzir a mesma análise. O modo mais direto é declarar o modelo para estimar os parâmetros separados para cada nível do fator, depois montar as matrizes correspondentes às funções lineares e submetê-las à glht() ou qualquer outra função que permita inferências para funções lineares de parâmetros. O pacote car, contrast e gmodels têm funções equivalentes à glht(), minha preferida.
##-----------------------------------------------------------------------------
m0 <- lm(Y~0+estudo/X, DATA) summary(m0)
m0$assign tapply(coef(m0), m0$assign, mean)
## Matriz de pesos. m <- rbind(rep(1, nlevels(DATA$estudo))/nlevels(DATA$estudo))
require(multcomp)
## Intercepto médio. X <- cbind(m, 0*m) summary(glht(m0, linfct=X), test=adjusted(type="none"))
## Inclinação média. X <- cbind(0*m, m) summary(glht(m0, linfct=X), test=adjusted(type="none"))
## As inclinações para cada nível saem no próprio summary() do modelo ## quando se declara o modelo na forma ~0+fator/numérica. Para ## conhecimento, pode ser feito assim.
M <- diag(nlevels(DATA$estudo)) X <- cbind(M, 0*M) summary(glht(m0, linfct=X), test=adjusted(type="fdr"))
X <- cbind(0*M, M) summary(glht(m0, linfct=X), test=adjusted(type="fdr"))
##-----------------------------------------------------------------------------
À disposição. Walmes.
_______________________________________________ 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 Souza -
walmes .