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: mac...@northnet.com.au -----Original Message----- From: R-help [mailto:r-help-boun...@r-project.org] On Behalf Of Wolfgang Raffelsberger Sent: Friday, 17 July 2015 03:41 To: r-help@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@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. ______________________________________________ R-help@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.
