Hi
your example is not reproducible. With ordinal regression the type of the y
values is important
sometimes an ordered factor is required.
Ordinal regression depends on your hypothesis see Ananth and Kleinbaum 1997
functions/packages to look at apart from ordinal
VGAM
polr::MASS
bayespolr::arm
lrm::rms
if you want to do GEE that is another matter.
Regards
Duncan
Duncan Mackay
Department of Agronomy and Soil Science
University of New England
Armidale NSW 2351
Email: home: mackay at northnet.com.au
-----Original Message-----
From: R-help [mailto:r-help-bounces at r-project.org] On Behalf Of Wolfgang
Raffelsberger
Sent: Friday, 17 July 2015 03:41
To: r-help at r-project.org
Subject: [R] (ordinal) logistic regression
Dear list,
I've been looking at previous posts on the list, but I haven't found any
close enough to my question/problem.
My data can be seen as a matrix of mutiple individuals (columns) with
(rather independent) measures (lines). Now based on supplemental
information, the individuals are organized in (multiple) ordered classes.
Now I would like to test for which type of measure (ie which line form my
matrix of data) the groups are distinct (eg different by group-mean). In
other words, I would like to see in which line of my input matrix the
measures for the groups of individuals associate to distinct group-values.
So I tried glm (or glm2) on each line of the matrix, but in my toy-example
(shown below) I'm surprised to get warnings about not finding convergence
in the "nice" toy-cases (ie distinct groups as I am looking for),e
even
with glm2 ! I see in such "nice" cases with glm() the
"Pr(>|z|)" is close
to 1, which in first sight is OK (since: H0 : coefficient =0), but I
suppose the test is not really set up right this way. When trying lrm (rms
package) I even get an error message (Unable to fit model using ?lrm.fit?)
with the "nice" cases.
In my real data with >4000 lines of data (ie >4000 glm tests) multiple
testing correction would transform everything from 1-p to end up at p=1, so
that?s another problem with this approach.
I suppose somehow I should transform my data (I don't see where I would
change the H0 ?) to obtain low and NOT high p-values (example below) in the
case I'm looking for, ie when group-means are distinct.
Any suggestions ?
Thank?s in advance,
Wolfgang
Here my toy-example :
datB1 <- c(12,14:16,18:21,20:22,20,22:24,19.5) # fit
partially/overlapping to 3grp model
datB2 <- c(11:15,21:25,31:36) # too beautiful to be real
...
datB3 <- c(15:12,11:15,12:14,15:12) # no fit to 3grp model
datB4 <- c(11:15,15:11,21:26) # grpA == grpB but grpA
!grpC
datB <- rbind(datB1,datB2,datB3,datB4)
set.seed(2015)
datB <- datB + round(runif(length(datB),-0.3,0.3),1) # add some noise
datB <- datB - rowMeans(datB) # centering
## here the definition of the groups
grpB <- gl(3,6,labels=LETTERS[1:3])[c(-6,-7)]
table(grpB)
## display
layout(1:4)
for(i in 1:4) plot(datB[i,],as.numeric(grpB))
## now the 'test'
glmLi <- function(dat,grp) {
## run glm : predict grp based on dat
dat <- data.frame(dat=dat,grp=grp)
glm(grp ~ dat, data=dat, family="binomial")}
logitB <- apply(datB,1,glmLi,grpB)
lapply(logitB,summary)
sapply(logitB,function(x) summary(x)$coefficients[,4]) # lines 1 & 2 are
designed to be 'positive' but give high p-values with convergence
problem
## for completness
sessionInfo()
R version 3.2.0 (2015-04-16)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 7 x64 (build 7601) Service Pack 1
locale:
[1] LC_COLLATE=French_France.1252 LC_CTYPE=French_France.1252
LC_MONETARY=French_France.1252
[4] LC_NUMERIC=C LC_TIME=French_France.1252
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] glm2_1.1.2 MASS_7.3-40 TinnRcom_1.0.18 formatR_1.2
svSocket_0.9-57
loaded via a namespace (and not attached):
[1] tools_3.2.0 svMisc_0.9-70 tcltk_3.2.0
[[alternative HTML version deleted]]
______________________________________________
R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.