[R] (ordinal) logistic regression

[Duncan Mackay]

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.