[R-br] OR e IC usando comandos clm e clm2
Luciane Maria Pilotto
lutipilotto em yahoo.com.br
Terça Julho 14 22:00:51 BRT 2015
Olá,
usei os comandos abaixo!
##########################################
load(file.choose())#dataframe:"id3.rda"
attach(id3)
dput(head(id3,10))
structure(list(regiao = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L), .Label = c("Norte", "Nordeste", "Sudeste", "Sul",
"Centro-Oeste"), class = "factor"), estado = structure(c(2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("Rondonia", "Acre",
"Amazonas", "Roraima", "Para", "Amapa", "Tocantins", "Maranhao",
"Piaui", "Ceara", "Rio Grande Do Norte", "Paraiba", "Pernambuco",
"Alagoas", "Sergipe", "Bahia", "Minas Gerais", "Espirito Santo",
"Rio De Janeiro", "Sao Paulo", "Parana", "Santa Catarina", "Rio Grande Do Sul",
"Mato Grosso Do Sul", "Mato Grosso", "Goias", "Distrito Federal"
), class = "factor"), cod_mun = c(1200401L, 1200401L, 1200401L,
1200401L, 1200401L, 1200401L, 1200401L, 1200401L, 1200401L, 1200401L
), setor = c("05000056", "05000056", "05000056", "05000056",
"05000056", "05000056", "05000056", "05000056", "05000062", "05000062"
), cap_int = structure(c(2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L), .Label = c("Interior", "Capital"), class = "factor"), idade = c(15L,
15L, 17L, 17L, 18L, 18L, 18L, 19L, 19L, 15L), getario = structure(c(3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L), .Label = c("5 anos", "12 anos",
"15 a 19 anos", "35 a 44 anos", "65 a 74 anos"), class = "factor"),
sexo = structure(c(2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L
), .Label = c("Male", "Female"), class = "factor"), grp_etni = structure(c(1L,
4L, 4L, 4L, 2L, 4L, 4L, 4L, 4L, 4L), .Label = c("Branca",
"Preta", "Amarela", "Parda", "Indigena", "Sem Informacao"
), class = "factor"), quest_01 = c(8L, 2L, 4L, 4L, 3L, 4L,
3L, 3L, 3L, 4L), quest_02 = c(4L, 4L, 2L, 2L, 4L, 1L, 3L,
1L, 3L, 2L), density = c(2, 0.5, 2, 2, 0.75, 4, 1, 3, 1,
2), quest_03 = c(3L, 4L, 6L, 5L, 5L, 3L, 3L, 3L, 6L, 5L),
quest_04 = structure(c(2L, 3L, 3L, 4L, NA, 2L, 1L, NA, 3L,
3L), .Label = c("Ate 250", "251 a 500", "501 a 1.500", "1.501 a 2.500",
"2.501 a 4.500", "4.501 a 9.500", "Mais de 9.500", "Nao sabe/Nao respondeu"
), class = "factor"), inc_percapita1 = c(46.875, 500, 250,
500, NA, 93.75, 41.6666679382324, NA, 333.333343505859, 250
), inc_percapita2 = c(46.875, 500, 250, 500, NA, 93.75, 41.6666679382324,
NA, 333.333343505859, 250), inc_sqrt1 = c(132.58251953125,
707.106811523438, 500, 1000, NA, 187.5, 72.1687850952148,
NA, 577.350280761719, 500), inc_sqrt2 = c(132.58251953125,
707.106811523438, 500, 1000, NA, 187.5, 72.1687850952148,
NA, 577.350280761719, 500), q05 = c(11L, 4L, 12L, 6L, 11L,
10L, 7L, 12L, 8L, 7L), quest_06 = structure(c(2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("Nao", "Sim", "Nao se aplica",
"Nao sabe"), class = "factor"), quest_07 = structure(c(1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L), .Label = c("Nao", "Sim",
"Nao se aplica", "Nao sabe"), class = "factor"), quest_08 = c(9L,
9L, 9L, 9L, 9L, 9L, 9L, 9L, 3L, 9L), quest_09 = structure(c(1L,
1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("Nao", "Sim",
"Nao se aplica", "Nao sabe"), class = "factor"), quest_10 = structure(c(5L,
5L, 2L, 3L, 3L, 2L, 2L, 3L, 3L, 3L), .Label = c("Menos de 1 ano",
"1 a 2 anos", ">3 anos", "Outros", "Nïa se aplica", "Nao sabe"
), class = "factor"), quest_11 = structure(c(5L, 5L, 1L,
1L, 2L, 1L, 1L, 1L, 1L, 1L), .Label = c("Publico", "Particular",
"Plano de Saude/Convenios", "Outros", "Nao se aplica", "Nao sabe"
), class = "factor"), quest_12 = structure(c(6L, 6L, 2L,
4L, 1L, 3L, 1L, 4L, 2L, 1L), .Label = c("Revisao/Prevencao",
"Dor", "Extracao", "Tratamento", "Outros", "Nao se aplica",
"Nao sabe"), class = "factor"), quest_13 = structure(c(NA,
NA, 2L, 3L, 2L, 3L, 2L, 2L, 2L, 2L), .Label = c("Muito bom",
"Bom", "Regular", "Ruim", "Muito Ruim", "Nao se aplica",
"Nao sabe"), class = "factor"), quest_14 = structure(c(3L,
2L, 3L, 4L, 2L, 2L, 4L, 3L, 4L, 2L), .Label = c("Muito satisfeito",
"Satisfeito", "Nem satisfeito/insatisfeito", "Insatisfeito",
"Muito Insatisfeito", "Nao sabe"), class = "factor"), quest_15 = structure(c(1L,
1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), .Label = c("Nao", "Sim",
"Nao se aplica", "Nao sabe"), class = "factor"), exame = c("1",
"1", "1", "1", "1", "1", "1", "1", "1", "1"), cpod = c(3L,
0L, 8L, 3L, 3L, 12L, 1L, 6L, 6L, 3L), p_sang = c(0L, 1L,
1L, 0L, 0L, 1L, 0L, 1L, 0L, 1L), p_calc = c(0L, 1L, 1L, 0L,
0L, 1L, 0L, 1L, 0L, 0L), cpi_max = structure(c(1L, 3L, 3L,
1L, 1L, 3L, 1L, 3L, 1L, 2L), .Label = c("Higido", "Sangramento",
"Calculo", "Bolsa 4-5 mm", "Bolsa 6 mm ou +", "4", "A", "X"
), class = "factor"), dai = c(21L, 32L, 23L, 21L, 19L, 18L,
17L, 23L, 34L, 25L), trauma = c(NA_integer_, NA_integer_,
NA_integer_, NA_integer_, NA_integer_, NA_integer_, NA_integer_,
NA_integer_, NA_integer_, NA_integer_), n_higido = c(26L,
28L, 24L, 24L, 25L, 17L, 27L, 25L, 24L, 25L), n_cariado = c(3L,
0L, 6L, 2L, 0L, 2L, 1L, 2L, 4L, 3L), n_restcar = c(0L, 0L,
0L, 0L, 0L, 0L, 0L, 2L, 1L, 0L), n_restaur = c(0L, 0L, 2L,
0L, 3L, 6L, 0L, 2L, 0L, 0L), n_perdcar = c(0L, 0L, 0L, 1L,
0L, 4L, 0L, 0L, 1L, 0L), n_perdout = c(0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L), perdidos = c(0, 0, 0, 1, 0, 4, 0, 0,
1, 0), usaprots = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0), usaproti = structure(c(1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), .Label = c("Não Usa",
"Uma Ponte Fixa", "Mais de 1 PF", "Prótese Parcial Removível",
"Prótese Fixa + Removível", "Prótese Total", "Sem Informação"
), class = "factor"), necprots = structure(c(1L, 1L, 3L,
1L, 1L, 1L, 1L, 1L, 1L, 1L), .Label = c("Não necessita",
"Prótese 1 elemento", "Mais de 1 elemento", "Combinação de próteses",
"Prótese Total", "Sem Informação"), class = "factor"), necproti = structure(c(1L,
1L, 4L, 2L, 1L, 3L, 1L, 1L, 2L, 1L), .Label = c("Não necessita",
"Prótese 1 elemento", "Mais de 1 elemento", "Combinação de próteses",
"Prótese Total", "Sem Informação"), class = "factor"), necprot = structure(c(1L,
1L, 3L, 2L, 1L, 2L, 1L, 1L, 2L, 1L), .Label = c("Não necessita",
"Parcial 1 maxilar", "Parcial 2 maxilar", "Total 1 maxilar",
"Parcial + Total", "Total 2 maxilar", "Sem Informação"), class = "factor"),
oidp = c(0L, 0L, 0L, 0L, 0L, 2L, 1L, 0L, 0L, 1L), f1 = c(1,
1, 1, 1, 1, 1, 1, 1, 1, 1), f2 = c(0.0518900007009506, 0.0518900007009506,
0.0518900007009506, 0.0518900007009506, 0.0518900007009506,
0.0518900007009506, 0.0518900007009506, 0.0518900007009506,
0.106059998273849, 0.106059998273849), f3 = c(0.0662299990653992,
0.0662299990653992, 0.0662299990653992, 0.0662299990653992,
0.0662299990653992, 0.0662299990653992, 0.0662299990653992,
0.0662299990653992, 0.138889998197556, 0.138889998197556),
f = c(0.00343999988399446, 0.00343999988399446, 0.00343999988399446,
0.00343999988399446, 0.00343999988399446, 0.00343999988399446,
0.00343999988399446, 0.00343999988399446, 0.0147299999371171,
0.0147299999371171), bwgr_et = c(291.019989013672, 291.019989013672,
291.019989013672, 291.019989013672, 291.019989013672, 291.019989013672,
291.019989013672, 291.019989013672, 67.879997253418, 67.879997253418
), cluster = c(120040105000056, 120040105000056, 120040105000056,
120040105000056, 120040105000056, 120040105000056, 120040105000056,
120040105000056, 120040105000062, 120040105000062), cluster2 = c("01120308",
"01120308", "01120308", "01120308", "01120308", "01120308",
"01120308", "01120308", "01120309", "01120309"), getni = c(1L,
4L, 4L, 4L, 2L, 4L, 4L, 4L, 4L, 4L), q04 = c(2L, 3L, 3L,
4L, NA, 2L, 1L, NA, 3L, 3L), q06 = c(2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L), q07 = c(1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 2L, 1L), q08 = c(0, 0, 0, 0, 0, 0, 0, 0, 3, 0), q10 = c(NA,
NA, 2L, 3L, 3L, 2L, 2L, 3L, 3L, 3L), q11 = c(NA, NA, 1L,
1L, 2L, 1L, 1L, 1L, 1L, 1L), q12 = c(NA, NA, 2L, 4L, 1L,
3L, 1L, 4L, 2L, 1L), q13 = structure(c(NA, NA, 2L, 3L, 2L,
3L, 2L, 2L, 2L, 2L), .Label = c("1", "2", "3", "4", "5"), class = "factor"),
q14 = c(3L, 2L, 3L, 4L, 2L, 2L, 4L, 3L, 4L, 2L), q15 = c(1L,
1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), edcat = c(3, 1, 4, 2,
3, 3, 2, 4, 2, 2), peso = c(291, 291, 291, 291, 291, 291,
291, 291, 68, 68), cariado = c(3L, 0L, 6L, 2L, 0L, 2L, 1L,
4L, 5L, 3L), getar = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1)), .Names = c("regiao",
"estado", "cod_mun", "setor", "cap_int", "idade", "getario",
"sexo", "grp_etni", "quest_01", "quest_02", "density", "quest_03",
"quest_04", "inc_percapita1", "inc_percapita2", "inc_sqrt1",
"inc_sqrt2", "q05", "quest_06", "quest_07", "quest_08", "quest_09",
"quest_10", "quest_11", "quest_12", "quest_13", "quest_14", "quest_15",
"exame", "cpod", "p_sang", "p_calc", "cpi_max", "dai", "trauma",
"n_higido", "n_cariado", "n_restcar", "n_restaur", "n_perdcar",
"n_perdout", "perdidos", "usaprots", "usaproti", "necprots",
"necproti", "necprot", "oidp", "f1", "f2", "f3", "f", "bwgr_et",
"cluster", "cluster2", "getni", "q04", "q06", "q07", "q08", "q10",
"q11", "q12", "q13", "q14", "q15", "edcat", "peso", "cariado",
"getar"), row.names = c("1241", "1242", "1243", "1244", "1245",
"1246", "1247", "1248", "1256", "1268"), class = "data.frame")
################################################################################
id3$q13<-as.factor(id3$q13)
m1 <- polr(q13 ~ q11 + q10 + q12 + edcat + q08 + q06 + q14, data=id3, Hess=TRUE)
summary(m1)
ordinal.or.display(m1)
####################################################################################
#Verificando os ajustes do modelo:
fm1 <- clm(q13 ~ q11 + q10 + q12 + edcat + q08 + q06 + q14, data=id3)
summary(fm1)
exp(coef(fm1))
(ci <- confint(fm1))
exp(cbind(OR = coef(fm1), ci))
Obs: o IC não está correto, está usando os coeficientes tb!!!
############################################################################################
nominal_test(fm1)# test partial proportional odds assumption
#A variável edcat violou o pressuposto de ordinalidade
##clm2 - partial proprotional odds
fm.nom <- clm2(q13 ~ q11 + q10 + q12 + q08 + q06 + q14, data=id3, nominal=~ edcat)
summary(fm.nom)
exp(coef(fm.nom))
exp(cbind(OR = coef(fm.nom), ci))
Obs: o IC não está correto, está usando os coeficientes tb!!!
Obrigada,
Luciane
__________________________________________________
Luciane Maria Pilotto
Mestre e Doutoranda em Saúde Bucal Coletiva - FO/UFRGS
NDE Odontologia - UNIVATES
Telefone: (51) 84512344
--------------------------------------------
Em qua, 8/7/15, Daniel Marcelino <dmsilva.br em gmail.com> escreveu:
Assunto: Re: [R-br] OR e IC usando comandos clm e clm2
Para: "R-br Lista" <r-br em listas.c3sl.ufpr.br>
Data: Quarta-feira, 8 de Julho de 2015, 21:05
Sem exemplo não dá para
falar muito. Mas deve haver uma função
`confint` para cada função de ajuste,
não?
Eu não conheço esse pacote, quando
preciso ajustar um modelo ordinal,
eu uso o
a função `` do pacote MASS.
modelo <- polr(y ~ x1 +
x2, data = dados, Hess=TRUE)
summary(modelo);
(ci <- confint(modelo));
confint.default(modelo); # ci assumindo
normalidade dos dados
#
Odds Ratios
exp(coef(modelo))
# or e ci
exp(cbind(or = coef(modelo), ci))
# Tabela
(ctable <- coef(summary(modelo)));
# calcular e salvar os p
values
p <- pnorm(abs(ctable[, "t
value"]), lower.tail = FALSE) * 2
(ctable <- cbind(ctable, "p value"
= p))
Daniel
2015-07-08 20:27 GMT-03:00 Luciane Maria
Pilotto <lutipilotto em yahoo.com.br>:
> Olá pessoal,
>
> Estou trabalhando com regressão
logística ordinal e realizando os ajustes com a função
"clm" do pacote ordinal e gostaria de obter o
intervalo de confiaça. De acordo com o Tutorial on fitting
Cumulative Link Models with the ordinal Package, de Rune
Haubo B Christensen (January 21, 2015) é possível rodar o
OR, mas não o intervalo de confiaça. O mesmo ocorre com a
função "clm2" para partial proportional odds.
>
> Agradeço qualquer
ajuda!!
>
>
> Luciane
>
>
_______________________________________________
> 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.
_______________________________________________
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.
Mais detalhes sobre a lista de discussão R-br