Re: [R] Simple question about function with glm
Thank you akk. I know it is not statistically sounded to check the distribution of response before glm. I will check the distribution of xmodel$residuals later on. About the program problem. It can print summary(xmodel) but not confint(xmodel) by amending my code as suggested by Bill Venables. Regards, CH On 5/7/07, Ben Bolker [EMAIL PROTECTED] wrote: Bill.Venables at csiro.au writes: Finally, I'm a bit puzzled why you use glm() when the simpler lm() would have done the job. You are fitting a linear model and do not need the extra paraphernaila that generalized linear models require. Bill Venables. Perhaps the original poster is confused about the difference between general (a la PROC GLM) and generalized (glm) linear models? The code is also a little puzzling because the same tests seem to be run whether p0.05 or not. Perhaps the code will eventually be written to log-transform the data if it fails the normality test? [ hint: ?boxcox in the MASS package might be a better way to go ] Ben Bolker __ R-help@stat.math.ethz.ch mailing list 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. -- The scientists of today think deeply instead of clearly. One must be sane to think clearly, but one can think deeply and be quite insane. Nikola Tesla http://www.macgrass.com __ R-help@stat.math.ethz.ch mailing list 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] different hights centering in one device region
Hello, I have a question. I creat a PDF file with four rows and two cols. Is it possible to: -create a plot regions with different hights (example: rows 1 2-4) -ro center an image in the whole width (example: rows 4) Thank's a lot. Felix example: Title |Text |Text | |Text |Text | - | | | | | | |image | image | | | | | | | - | | | | | | |image | image | | | | | | | - | | | | | | | | | |image | | | | | | | | | | - R-Code: --- pdf(pPDF, height=13.7, paper=special) par(oma=c(0,0,1,0), mfrow=c(4,2) ) #Text field plot.new() text(0, 0.6, pos=4, cex=1.2, paste(Text) ) text(0, 0.4, pos=4, cex=1.2, paste(Text) ) plot.new() text(0, 0.6, pos=4, cex=1.2, paste(Text) ) text(0, 0.4, pos=4, cex=1.2, paste(Text) ) image(zMEDIAN_1) image(zMEDIAN_2) image(zSD_1) #Title par(mfrow=c(1,1), oma=c(0,0,1,0)) mtext(Darstellung der Differenz zweier Exportfiles V0.03, 3, outer = TRUE, cex = par(cex.main)) dev.off() __ R-help@stat.math.ethz.ch mailing list 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.
Re: [R] different hights centering in one device region
Felix Wave wrote: Hello, I have a question. I creat a PDF file with four rows and two cols. Is it possible to: -create a plot regions with different hights (example: rows 1 2-4) -ro center an image in the whole width (example: rows 4) Check the help page for layout(). Read it carefully, the mechanism is somewhat elaborate. You probably need to do it by parcelling out a 4x4 matrix with an allocation matrix like this: 1122 3344 5566 0770 Thank's a lot. Felix example: Title |Text |Text | |Text |Text | - | | | | | | |image| image | | | | | | | - | | | | | | |image| image | | | | | | | - | | | | | | | | | |image | | | | | | | | | | - R-Code: --- pdf(pPDF, height=13.7, paper=special) par(oma=c(0,0,1,0), mfrow=c(4,2) ) #Text field plot.new() text(0, 0.6, pos=4, cex=1.2, paste(Text) ) text(0, 0.4, pos=4, cex=1.2, paste(Text) ) plot.new() text(0, 0.6, pos=4, cex=1.2, paste(Text) ) text(0, 0.4, pos=4, cex=1.2, paste(Text) ) image(zMEDIAN_1) image(zMEDIAN_2) image(zSD_1) #Title par(mfrow=c(1,1), oma=c(0,0,1,0)) mtext(Darstellung der Differenz zweier Exportfiles V0.03, 3, outer = TRUE, cex = par(cex.main)) dev.off() __ R-help@stat.math.ethz.ch mailing list 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. -- O__ Peter Dalgaard Øster Farimagsgade 5, Entr.B c/ /'_ --- Dept. of Biostatistics PO Box 2099, 1014 Cph. K (*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918 ~~ - ([EMAIL PROTECTED]) FAX: (+45) 35327907 __ R-help@stat.math.ethz.ch mailing list 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.
Re: [R] to draw a smooth arc
There is now an xspline() function in R-devel, with an example showing how to add arrows. I thought a bit more about a 'circular arc' function, but there really is a problem with that. Few R plot regions have a 1:1 aspect ratio including some that are intended to do so (see the rw-FAQ). symbols() is designed to draw circles in device coordinates, but attempting to specify circular arcs by endpoints in user coordinates is fraught. On Wed, 2 May 2007, Paul Murrell wrote: Hi Paulo Barata wrote: Dr. Murrell and all, One final suggestion: a future function arc() in package graphics, with centre-radius-angle parameterisation, could also include an option to draw arrows at either end of the arc, as one can find in function arrows(). ... and in grid.xspline() and grid.curve(). Paul Thank you. Paulo Barata --- Paul Murrell wrote: Hi Paulo Barata wrote: Dr. Snow and Prof. Ripley, Dr. Snow's suggestion, using clipplot (package TeachingDemos), is maybe a partial solution to the problem of drawing an arc of a circle (as long as the line width of the arc is not that large, as pointed out by Prof. Ripley). If the arc is symmetrical around a vertical line, then it is not so difficult to draw it that way. But an arc that does not have this kind of symmetry would possibly require some geometrical computations to find the proper rectangle to be used for clipping. I would like to suggest that in a future version of R some function be included in the graphics package to draw smooth arcs with given center, radius, initial and final angles. I suppose that the basic ingredients are available in function symbols (graphics). Just to back up a few previous posts ... There is something like this facility already available via the grid.xspline() function in the grid package. This provides very flexible curve drawing (including curves very close to Bezier curves) based on the X-Splines implemented in xfig. The grid.curve() function provides a convenience layer that allows for at least certain parameterisations of arcs (you specify the arc end points and the angle). These functions are built on functionality within the core graphics engine, so exposing a similar interface (e.g., an xspline() function) within traditional graphics would be relatively straightforward. The core functionality draws the curves as line segments (but automatically figures out how many segments to use so that the curve looks smooth); it does NOT call curve-drawing primitives in the graphics device (like PostScript's curveto). In summary: there is some support for smooth curves, but we could still benefit from a specific arc() function with the standard centre-radius-angle parameterisation and we could also benefit from exposing the native strengths of different graphics devices (rather than the current lowest-common-denominator approach). Paul Thank you very much. Paulo Barata (Rio de Janeiro - Brazil) --- Prof Brian Ripley wrote: On Tue, 1 May 2007, Greg Snow wrote: Here is an approach that clips the circle you like from symbols down to an arc (this will work as long as the arc is less than half a circle, for arcs greater than half a circle, you could draw the whole circle then use this to draw an arc of the bacground color over the section you don't want): library(TeachingDemos) plot(-5:5, -5:5, type='n') clipplot( symbols(0,0,circles=2, add=TRUE), c(0,5), c(0,5) ) I had considered this approach: clipping a circle to a rectangle isn't strictly an arc, as will be clear if the line width is large. Consider clipplot(symbols(0, 0 ,circles=2, add=TRUE, lwd=5), c(-1,5), c(-1,5)) Note too that what happens with clipping is device-dependent. If R's internal clipping is used, the part-circle is converted to a polygon. __ R-help@stat.math.ethz.ch mailing list 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. -- Brian D. Ripley, [EMAIL PROTECTED] Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UKFax: +44 1865 272595 __ R-help@stat.math.ethz.ch mailing list 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.
Re: [R] Neural Nets (nnet) - evaluating success rate of predictions
On 5/6/07, nathaniel Grey [EMAIL PROTECTED] wrote: Hello R-Users, I have been using (nnet) by Ripley to train a neural net on a test dataset, I have obtained predictions for a validtion dataset using: PP-predict(nnetobject,validationdata) Using PP I can find the -2 log likelihood for the validation datset. However what I really want to know is how well my nueral net is doing at classifying my binary output variable. I am new to R and I can't figure out how you can assess the success rates of predictions. table(PP, binaryvariable) should get you started. Also if you're using nnet with random starts, I strongly suggest taking the best out of several hundred (or maybe thousand) trials - it makes a big difference! Hadley __ R-help@stat.math.ethz.ch mailing list 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] contour and density
Hello, I would like to draw the contour of a posterior distribution. However I only have simulated points from this posterior distribution. So I estimate the density of my posterior distribution : density(points) I get the mean, the max, the median I know I can plot this density estimate. However if I try : contour(z=density(points)) It does not work. How can I do that ? Besides how I can get the value of my density estimate at one point ? when I try density(points)(vector), R tells me it is not a function. Thank you very much. _ météo et bien plus encore ! [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list 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] [R-pkgs] New Package Reliability
Dear R useRs, A new package 'Reliability' is now available on CRAN. It is mainly a set of functions functions for estimating parameters in software reliability models. Only infinite failure models are implemented so far. This is the first version of the package. The canonical reference is: J.D. Musa, A. Iannino, and K. Okumoto. Software Reliability: Measurement, Prediction, Application. McGraw-Hill, 1987. Michael R. Lyu. Handbook of Software Realibility Engineering. IEEE Computer Society Press, 1996. http://www.cse.cuhk.edu.hk/~lyu/book/reliability/ Suggestions, bug reports and other comments are very welcome. enjoy and best regards Andreas ___ R-packages mailing list [EMAIL PROTECTED] https://stat.ethz.ch/mailman/listinfo/r-packages __ R-help@stat.math.ethz.ch mailing list 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.
Re: [R] Neural Nets (nnet) - evaluating success rate of predictions
well, how to do you know which ones are the best out of several hundreds? I will average all results out of several hundreds. On 5/7/07, hadley wickham [EMAIL PROTECTED] wrote: On 5/6/07, nathaniel Grey [EMAIL PROTECTED] wrote: Hello R-Users, I have been using (nnet) by Ripley to train a neural net on a test dataset, I have obtained predictions for a validtion dataset using: PP-predict(nnetobject,validationdata) Using PP I can find the -2 log likelihood for the validation datset. However what I really want to know is how well my nueral net is doing at classifying my binary output variable. I am new to R and I can't figure out how you can assess the success rates of predictions. table(PP, binaryvariable) should get you started. Also if you're using nnet with random starts, I strongly suggest taking the best out of several hundred (or maybe thousand) trials - it makes a big difference! Hadley __ R-help@stat.math.ethz.ch mailing list 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. -- WenSui Liu A lousy statistician who happens to know a little programming (http://spaces.msn.com/statcompute/blog) __ R-help@stat.math.ethz.ch mailing list 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.
Re: [R] Neural Nets (nnet) - evaluating success rate of predictions
Pick the one with the lowest error rate on your training data? Hadley On 5/7/07, Wensui Liu [EMAIL PROTECTED] wrote: well, how to do you know which ones are the best out of several hundreds? I will average all results out of several hundreds. On 5/7/07, hadley wickham [EMAIL PROTECTED] wrote: On 5/6/07, nathaniel Grey [EMAIL PROTECTED] wrote: Hello R-Users, I have been using (nnet) by Ripley to train a neural net on a test dataset, I have obtained predictions for a validtion dataset using: PP-predict(nnetobject,validationdata) Using PP I can find the -2 log likelihood for the validation datset. However what I really want to know is how well my nueral net is doing at classifying my binary output variable. I am new to R and I can't figure out how you can assess the success rates of predictions. table(PP, binaryvariable) should get you started. Also if you're using nnet with random starts, I strongly suggest taking the best out of several hundred (or maybe thousand) trials - it makes a big difference! Hadley __ R-help@stat.math.ethz.ch mailing list 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. -- WenSui Liu A lousy statistician who happens to know a little programming (http://spaces.msn.com/statcompute/blog) __ R-help@stat.math.ethz.ch mailing list 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.
Re: [R] decimal values
options(digits = 2) ?options Atenciosamente, Ana Patricia Martins --- Serviço Métodos Estatísticos Departamento de Metodologia Estatística INE - Portugal Telef: 218 426 100 - Ext: 3210 E-mail: [EMAIL PROTECTED] -Original Message- From: elyakhlifi mustapha [mailto:[EMAIL PROTECTED] Sent: sexta-feira, 4 de Maio de 2007 13:59 To: R-help@stat.math.ethz.ch Subject: [R] decimal values hello, how can I do to drop decimal after the comma please for example for tthis line print(P) [1] 62.00 1.00 7.661290 5.20 17.10 2.318801 how canI do to keep only 62 1 7.66 5.2 17.12.32 thanks __ ble contre les messages non sollicités [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list 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.
Re: [R] Predicted Cox survival curves - factor coding problems..
The combination of survfit, coxph, and factors is getting confused. It is not smart enough to match a new data frame that contains a numeric for sitenew to a fit that contained that variable as a factor. (Perhaps it should be smart enough to at least die gracefully -- but it's not). The simple solution is to not use factors. site1 - 1*(coxsnps$sitenew==1) site2 - 1*(coxsnps$sitenew==2) test1 - coxph(Surv(time, censor) ~ snp1 + sex + site1 + site2 + gene + eth.self + strata(edu), data= coxsnps) output profile1 - data.frame(snp1=c(0,1), site2=c(0,0), sex=c(0,0), site1=c(0,0), site2=c(0,0), geno=c(0,0) eth.self=c(0,0)) plot(survfit(test1, newdata=profile1)) Note that you do not have to explicitly make edu a factor. Since it is included in a strata statement, the coxph routine must treat it as discrete groups. Terry Therneau __ R-help@stat.math.ethz.ch mailing list 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.
Re: [R] A function for raising a matrix to a power?
At 15:48 06/05/2007, Ron E. VanNimwegen wrote: Hi, Is there a function for raising a matrix to a power? For example if you like to compute A%*%A%*%A, is there an abbreviation similar to A3? Atte Tenkanen I may be revealing my ignorance here, but is MatrixExp in the msm package (available from CRAN) not relevant here? [snip other suggestion] Michael Dewey http://www.aghmed.fsnet.co.uk __ R-help@stat.math.ethz.ch mailing list 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] finding trajectories
hello R-users, i like to know if there exists a tool in R to find different trajectories of a variable in a dataset with three resp. four points of measurement. i search something like PROC TRAJ for SAS, but i don't have SAS nor do i know to use it. other alternatives would be M-plus or latent gold. (but actually i don't want to learn the use of M-plus or latent gold.) what i think as an important element of PROC TRAJ is the calculation of BIC, a criterion to decide what number of trajectories fitting best the data. thanks for every help. -- Egon Werlen Fachpsychologe fuer Gesundheitspsychologie FSP Forschungszentrum fuer Rehabilitations- und Gesundheitspsychologie Universitaet Freiburg Englisberg 7 CH-1763 Granges-Paccot Tel: +41 (0)26 300 76 60 Fax: +41 (0)26 300 96 34 __ R-help@stat.math.ethz.ch mailing list 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.
Re: [R] Predicted Cox survival curves - factor coding problems..
On Mon, 7 May 2007, Terry Therneau wrote: The combination of survfit, coxph, and factors is getting confused. It is not smart enough to match a new data frame that contains a numeric for sitenew to a fit that contained that variable as a factor. (Perhaps it should be smart enough to at least die gracefully -- but it's not). The 'standard' model-fitting functions in R do make an attempt to match the new data to that used for fitting, or die gracefully. Perhaps Thomas could look into adding this to survift and coxph (see http://developer.r-project.org/model-fitting-functions.txt). The simple solution is to not use factors. site1 - 1*(coxsnps$sitenew==1) site2 - 1*(coxsnps$sitenew==2) test1 - coxph(Surv(time, censor) ~ snp1 + sex + site1 + site2 + gene + eth.self + strata(edu), data= coxsnps) output profile1 - data.frame(snp1=c(0,1), site2=c(0,0), sex=c(0,0), site1=c(0,0), site2=c(0,0), geno=c(0,0) eth.self=c(0,0)) plot(survfit(test1, newdata=profile1)) Note that you do not have to explicitly make edu a factor. Since it is included in a strata statement, the coxph routine must treat it as discrete groups. Terry Therneau -- Brian D. Ripley, [EMAIL PROTECTED] Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UKFax: +44 1865 272595 __ R-help@stat.math.ethz.ch mailing list 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.
Re: [R] summing values according to a factor
Howdy!
I guess what you want to do is compare Q1/T1 among the sections? If you
want to compute the sum of Q1/T1 by Section, you can do something like:
sum.by.section - with(mydata, tapply(Q1/T1, section, sum))
Substitute sum with anything you want to compute.
Cheers,
Andy
From: Salvatore Enrico Indiogine
Greetings!
I have exam scores of students in several sections. The data
looks like:
StuNum Section Q1 T1
111 502 45 123
112 502 23123
113 503 58123
114 504 63 123
115 504 83 123
..
where Q1 is the score for question 1 and T1 is the maximum possible
score for question 1
I need to check whether the section has an effect on the scores. I
thought about using chisq.test and calculate the sums of scores per
section.
I think that I have to use apply() but I am lost here.
Thanks in advance,
Enrico
--
Enrico Indiogine
Mathematics Education
Texas AM University
[EMAIL PROTECTED]
__
R-help@stat.math.ethz.ch mailing list
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.
--
Notice: This e-mail message, together with any attachments,...{{dropped}}
__
R-help@stat.math.ethz.ch mailing list
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] Representing a statistic as a colour on a 2d plot
Hello. I have a 2d plot which looks like this: http://www.nabble.com/file/8242/1.JPG This plot is derived from a file that holds statistics about each point on the plot and looks like this: abc d e a00.4980.4730.524 0.528 b 0.498 0 0 0 0 c 0.473 0 0 0 0 d 0.524 0 0 0 0 e 0.528 0 0 0 0 However, I have another file called 2.txt, with the following contents: a b c d e 0.5 0.7 0.32 0.34 0.01 What I would like to know is how do I convert these values in 2.txt to colours or colour intensities so that the x's in the diagram above can be colour coded as such. Many thanks. -- View this message in context: http://www.nabble.com/Representing-a-statistic-as-a-colour-on-a-2d-plot-tf3703885.html#a10357828 Sent from the R help mailing list archive at Nabble.com. __ R-help@stat.math.ethz.ch mailing list 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] creating a new column
hie l would like to create a 6th column actual surv time from the following
data
the condition being
if censoringTimesurvivaltime then actual survtime =survival time
else actual survtime =censoring time
the code l used to create the data is
s=2
while(s!=0){ n=20
m-matrix(nrow=n,ncol=4)
colnames(m)=c(treatmentgrp,strata,censoringTime,survivalTime)
for(i in 1:20)
m[i,]-c(sample(c(1,2),1,replace=TRUE),sample(c(1,2),1,replace=TRUE),rexp(1,.007),rexp(1,.002))
m-cbind(m,0)
m[m[,3]m[,4],5]-1
colnames(m)[5]-censoring
print(m)
s=s-1
treatmentgrp strata censoringTime survivalTime censoring
[1,] 1 1 1.0121591137.80922 0
[2,] 2 2 32.971439 247.21786 0
[3,] 2 1 85.758253 797.04949 0
[4,] 1 1 16.999171 78.92309 0
[5,] 2 1 272.909896 298.21483 0
[6,] 1 2 138.230629 935.96765 0
[7,] 2 2 91.529859 141.08405 0
l keep getting an error message when i try to create the 6th column
-
[[alternative HTML version deleted]]
__
R-help@stat.math.ethz.ch mailing list
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.
Re: [R] Representing a statistic as a colour on a 2d plot
check ?rainbow to generate the colours (which also shows you other related functions like 'heat.colors' that you may also find useful. To plot the coloured points, check ?points. hope this helps a bit Jose Quoting mister_bluesman [EMAIL PROTECTED]: Hello. I have a 2d plot which looks like this: http://www.nabble.com/file/8242/1.JPG This plot is derived from a file that holds statistics about each point on the plot and looks like this: abc d e a00.4980.4730.524 0.528 b 0.498 0 0 0 0 c 0.473 0 0 0 0 d 0.524 0 0 0 0 e 0.528 0 0 0 0 However, I have another file called 2.txt, with the following contents: a b c d e 0.5 0.7 0.32 0.34 0.01 What I would like to know is how do I convert these values in 2.txt to colours or colour intensities so that the x's in the diagram above can be colour coded as such. Many thanks. -- View this message in context: http://www.nabble.com/Representing-a-statistic-as-a-colour-on-a-2d-plot-tf3703885.html#a10357828 Sent from the R help mailing list archive at Nabble.com. __ R-help@stat.math.ethz.ch mailing list 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. -- Dr. Jose I. de las Heras Email: [EMAIL PROTECTED] The Wellcome Trust Centre for Cell BiologyPhone: +44 (0)131 6513374 Institute for Cell Molecular BiologyFax: +44 (0)131 6507360 Swann Building, Mayfield Road University of Edinburgh Edinburgh EH9 3JR UK __ R-help@stat.math.ethz.ch mailing list 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.
Re: [R] R² in a non-linear regresion analisys
Douglas Bates: It may seem obvious how to define the multiple correlation coefficient R^2 for a non-linear regression model but it's not. A couple of articles explaining why, that may be of interest: Anderson-Sprecher R. (1994). ‘Model comparisons and R²’. The American Statistician, volume 48, no. 2, pages 113–117. http://dx.doi.org/10.2307/2684259 Kvålseth T.O. (1985). ‘Cautionary note about R²’. The American Statistician, volume 39, no. 4, pages 279–285. http://dx.doi.org/10.2307/2683704 -- Karl Ove Hufthammer __ R-help@stat.math.ethz.ch mailing list 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.
Re: [R] creating a new column
To create 6th column in the matrix m, you should use the cbind function.
To calculate the vector of pairwise min or max values, you should use the
pmin and pmax functions:
act.surv.time-pmin(m[,censoringTime],m[,survivalTime])
m-cbind(m,act.surv.time)
raymond chiruka wrote:
hie l would like to create a 6th column actual surv time from the
following data
the condition being
if censoringTimesurvivaltime then actual survtime =survival time
else actual survtime =censoring time
the code l used to create the data is
s=2
while(s!=0){ n=20
m-matrix(nrow=n,ncol=4)
colnames(m)=c(treatmentgrp,strata,censoringTime,survivalTime)
for(i in 1:20)
m[i,]-c(sample(c(1,2),1,replace=TRUE),sample(c(1,2),1,replace=TRUE),rexp(1,.007),rexp(1,.002))
m-cbind(m,0)
m[m[,3]m[,4],5]-1
colnames(m)[5]-censoring
print(m)
s=s-1
treatmentgrp strata censoringTime survivalTime censoring
[1,] 1 1 1.0121591137.80922 0
[2,] 2 2 32.971439 247.21786 0
[3,] 2 1 85.758253 797.04949 0
[4,] 1 1 16.999171 78.92309 0
[5,] 2 1 272.909896 298.21483 0
[6,] 1 2 138.230629 935.96765 0
[7,] 2 2 91.529859 141.08405 0
l keep getting an error message when i try to create the 6th column
-
__
R-help@stat.math.ethz.ch mailing list
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.
--
View this message in context:
http://www.nabble.com/creating-a-new-column-tf3704043.html#a10358728
Sent from the R help mailing list archive at Nabble.com.
__
R-help@stat.math.ethz.ch mailing list
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] looking for equivalent of matlab's medfilt1 function
Dear all, I have several files with Matlab code, which I am translating to R. For the zero-level approach, I took the very old shell script from R-help archives, which has made some obvious leg-work such as replacement of = with -. Now I am translating indexing, matrix operations and function call using this table http://37mm.no/mpy/octave-r.html The problem is, I cannot find the R equivalent of the matlab's function medfilt1, performing 1-dimensional median filtering of a vector. Its summary is here http://www-ccs.ucsd.edu/matlab/toolbox/signal/medfilt1.html -- View this message in context: http://www.nabble.com/looking-for-equivalent-of-matlab%27s-medfilt1-function-tf3704211.html#a10358868 Sent from the R help mailing list archive at Nabble.com. __ R-help@stat.math.ethz.ch mailing list 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.
Re: [R] looking for equivalent of matlab's medfilt1 function
Vladimir == Vladimir Eremeev [EMAIL PROTECTED] on Mon, 7 May 2007 07:58:48 -0700 (PDT) writes: Vladimir Dear all, Vladimir I have several files with Matlab code, which I am translating to R. Vladimir For the zero-level approach, I took the very old Vladimir shell script from R-help archives, which has made Vladimir some obvious leg-work such as replacement of = Vladimir with -. Vladimir Now I am translating indexing, matrix operations and function call using Vladimir this table Vladimir http://37mm.no/mpy/octave-r.html You should also look at the 'matlab' package which defines quite a few R functions such as eyes(), zeros(), repmat(), to behave as the Matlab functions do. Vladimir The problem is, I cannot find the R equivalent of the matlab's function Vladimir medfilt1, performing 1-dimensional median filtering of a vector. Its summary Vladimir is here http://www-ccs.ucsd.edu/matlab/toolbox/signal/medfilt1.html To statisticians, this has been known as running medians, thanks to John Tukey. The smooth() function contains smart variations of running median of 3 which seems to be the matlab default. For 'k 3', I'd recommend the fast runmed(x, k) function which also has a bit more sophisticated end-point handling than Matlab's medfilt1() seems to provide. Martin __ R-help@stat.math.ethz.ch mailing list 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.
Re: [R] factanal AIC?
The log likelihood is
const + n/2* [ log(det(Sigma^-1)) - trace(Sigma^-1*S) ] where Sigma is the
population covariance matrix
-Original Message-
From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Spencer Graves
Sent: Friday, May 04, 2007 9:20 PM
To: R-devel mailing list
Cc: Jens Oehlschlägel; r-help@stat.math.ethz.ch
Subject: Re: [R] factanal AIC?
Dear R Developers:
What is the proper log(likelihood) for 'factanal'? I believe it should
be something like the following:
(-n/2)*(k*(log(2*pi)+1)+log(det(S)))
or without k*(log(2*pi)-1):
(-n/2)*log(det(S)),
where n = the number of (multivariate) observations.
The 'factanal' help file say the factanal component criteria:
The results of the optimization: the value of the negative log-likelihood and
information on the iterations used.
However, I'm not able to get this. If it's a log(likelihood), then
replacing a data set m1 by rbind(m1, m1) should double the log(likelihood).
However it has no impact on it. Consider the following modification of the
first example in the 'factanal' help page:
v1 - c(1,1,1,1,1,1,1,1,1,1,3,3,3,3,3,4,5,6)
v2 - c(1,2,1,1,1,1,2,1,2,1,3,4,3,3,3,4,6,5)
v3 - c(3,3,3,3,3,1,1,1,1,1,1,1,1,1,1,5,4,6)
v4 - c(3,3,4,3,3,1,1,2,1,1,1,1,2,1,1,5,6,4)
v5 - c(1,1,1,1,1,3,3,3,3,3,1,1,1,1,1,6,4,5)
v6 - c(1,1,1,2,1,3,3,3,4,3,1,1,1,2,1,6,5,4)
m1 - cbind(v1,v2,v3,v4,v5,v6)
tst - factanal(m1, factors=3)
fit1 - factanal(m1, factors=3)
fit2 - factanal(rbind(m1, m1), factors=3) fit1$criteria[objective]
objective
0.4755156
fit2$criteria[objective]
objective
0.4755156
From the following example, it seems that the k*(log(2*pi)-1) term is
omitted:
x2 - c(-1,1)
X2.4 - as.matrix(expand.grid(x2, x2, x2, x2)) factanal(X2.4, 1)$criteria
objective counts.function counts.gradient
0 7 7
However, I can't get the 'fit1$criteria' above, even ignoring the sample
size. Consider the following:
S3 - tcrossprod(fit1$loadings)+diag(fit1$uniquenesses)
log(det(S3))
[1] -6.725685
Shouldn't this be something closer to the 0.4755 for fit1$criteria?
Thanks,
Spencer Graves
JENS: For your purposes, I suggest you try to compute
(-n/2)*log(det(tcrossprod(fit$loadings)+diag(fit$uniquenesses)). If you do
this, you might get more sensible results with your examples.
Hope this helps.
Spencer Graves
Jens Oehlschlägel wrote:
Dear list members,
Could any expert on factor analysis be so kind to explain how to calculate
AIC on the output of factanal. Do I calculate AIC wrong or is
factanal$criteria[objective] not a negative log-likelihood?
Best regards
Jens Oehlschlägel
The AIC calculated using summary.factanal below don't appear correct to me:
n items factors total.df rest.df model.df LL AIC
AICc BIC
1 100020 1 210 170 40 -5.192975386 90.38595
93.80618 286.6962
2 100020 2 210 151 59 -3.474172303 124.94834
132.48026 414.5059
3 100020 3 210 133 77 -1.745821627 157.49164
170.51984 535.3888
4 100020 4 210 116 94 -0.120505369 188.24101
207.97582 649.5700
5 100020 5 210 100 110 -0.099209921 220.19842
247.66749 760.0515
6 100020 6 210 85 125 -0.072272695 250.14455
286.18574 863.6140
7 100020 7 210 71 139 -0.054668588 278.10934
323.36515 960.2873
8 100020 8 210 58 152 -0.041708051 304.08342
358.99723 1050.0622
9 100020 9 210 46 164 -0.028686298 328.05737
392.87174 1132.9292
10 100020 10 210 35 175 -0.015742783 350.03149
424.78877 1208.8887
11 100020 11 210 25 185 -0.007007901 370.01402
454.55947 1277.9487
12 100020 12 210 16 194 -0.005090725 388.01018
481.99776 1340.1147
summary.factanal - function(object, ...){
if (inherits(object, try-error)){
c(n=NA, items=NA, factors=NA, total.df=NA, rest.df=NA, model.df=NA,
LL=NA, AIC=NA, AICc=NA, BIC=NA)
}else{
n - object$n.obs
p - length(object$uniquenesses)
m - object$factors
model.df - (p*m) + (m*(m+1))/2 + p - m^2
total.df - p*(p+1)/2
rest.df - total.df - model.df # = object$dof
LL - -as.vector(object$criteria[objective])
k - model.df
aic - 2*k - 2*LL
aicc - aic + (2*k*(k+1))/(n-k-1)
bic - k*log(n) - 2*LL
c(n=n, items=p, factors=m, total.df=total.df, rest.df=rest.df,
model.df=model.df, LL=LL, AIC=aic, AICc=aicc, BIC=bic)
}
}
multifactanal - function(factors=1:3, ...){
names(factors) - factors
ret - lapply(factors, function(factors){
try(factanal(factors=factors, ...))
})
class(ret) - multifactanal
ret
}
Re: [R] Neural Nets (nnet) - evaluating success rate of predictions
Folks:
If I understand correctly, the following may be pertinent.
Note that the procedure:
min.nnet = nnet[k] such that error rate of nnet[k] = min[i] {error
rate(nnet(training data) from ith random start) }
does not guarantee a classifier with a lower error rate on **new** data than
any single one of the random starts. That is because you are using the same
training set to choose the model (= nnet parameters) as you are using to
determine the error rate. I know it's tempting to think that choosing the
best among many random starts always gets you a better classifier, but it
need not. The error rate on the training set for any classifier -- be it a
single one or one derived in some way from many -- is a biased estimate of
the true error rate, so that choosing a classifer on this basis does not
assure better performance for future data. In particular, I would guess that
choosing the best among many (hundreds/thousands) random starts is probably
almost guaranteed to produce a poor predictor (ergo the importance of
parsimony/penalization). I would appreciate comments from anyone, pro or
con, with knowledge and experience of these things, however, as I'm rather
limited on both.
The simple answer to the question of obtaining the error rate using
validation data is: Do whatever you like to choose/fit a classifier on the
training set. **Once you are done,** the estimate of your error rate is the
error rate you get on applying that classifier to the validation set. But
you can do this only once! If you don't like that error rate and go back to
finding a a better predictor in some way, then your validation data have now
been used to derive the classifier and thus has become part of the training
data, so any further assessment of the error rate of a new classifier on it
is now also a biased estimate. You need yet new validation data for that.
Of course, there are all sort of cross validation schemes one can use to
avoid -- or maybe mitigate -- these issues: most books on statistical
classification/machine learning discuss this in detail.
Bert Gunter
Genentech Nonclinical Statistics
-Original Message-
From: [EMAIL PROTECTED]
[mailto:[EMAIL PROTECTED] On Behalf Of hadley wickham
Sent: Monday, May 07, 2007 5:26 AM
To: Wensui Liu
Cc: r-help@stat.math.ethz.ch
Subject: Re: [R] Neural Nets (nnet) - evaluating success rate of predictions
Pick the one with the lowest error rate on your training data?
Hadley
On 5/7/07, Wensui Liu [EMAIL PROTECTED] wrote:
well, how to do you know which ones are the best out of several hundreds?
I will average all results out of several hundreds.
On 5/7/07, hadley wickham [EMAIL PROTECTED] wrote:
On 5/6/07, nathaniel Grey [EMAIL PROTECTED] wrote:
Hello R-Users,
I have been using (nnet) by Ripley to train a neural net on a test
dataset, I have obtained predictions for a validtion dataset using:
PP-predict(nnetobject,validationdata)
Using PP I can find the -2 log likelihood for the validation datset.
However what I really want to know is how well my nueral net is doing
at classifying my binary output variable. I am new to R and I can't figure
out how you can assess the success rates of predictions.
table(PP, binaryvariable)
should get you started.
Also if you're using nnet with random starts, I strongly suggest
taking the best out of several hundred (or maybe thousand) trials - it
makes a big difference!
Hadley
__
R-help@stat.math.ethz.ch mailing list
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.
--
WenSui Liu
A lousy statistician who happens to know a little programming
(http://spaces.msn.com/statcompute/blog)
__
R-help@stat.math.ethz.ch mailing list
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@stat.math.ethz.ch mailing list
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.
Re: [R] heatmap and phylogram / dendogram ploting problem, ape package
Emmanuel Paradis, the maintainer of the ape package was very helpful in solving this problem. It seems that it heatmap does not reorder the rows, so you must reorder them or change the heatmap code to do so. The heatmap maintainers will document this, but not change the behavior. The following fix works for me: On 4/30/07, Emmanuel wrote: . . . Thanks for your interesting case study. I found out that, apparently, your problem comes from a bug in heatmap: it assumes that the row names of the matrix are in the same order than the labels of the dendrogram. In your case, the following fix seems to work. You have to edit heatmap with fix(heatmap), find the line: x - x[rowInd, colInd] and change it for: x - x[labels(ddr)[rowInd], colInd] I think this will not solve the problem in all cases. . . On 4/24/07, jamesrh [EMAIL PROTECTED] wrote: I am having trouble displaying a dendrogram of evolutionary relationships (a phylogram imported from the ape package) as the vertical component of a heatmap, but keeping the hierarchical clustering of the horizontal component. The relationships of the vertical component in the generated heatmap are not that of the dendrogram, although the ordering is. In more detail, I am attempting to generate a heatmap from a table that contains the abundance of different bacteria at different locations, with a dendrogram that groups the environments by the pattern of bacterial abundance. This is easy, thanks to a nice code snippet at the R Graph Gallery (http://addictedtor.free.fr/graphiques/RGraphGallery.php?graph=66): env - read.table(env.dat) x - as.matrix(env) hv - heatmap(x, col = cm.colors(256), scale=none, xlab = Environments, ylab= Species, main = heatmap of species present in environments) However, instead of a dendrogram that groups the rows (vertical axis) by the abundance pattern (as above), I would like to force it to order and display a dendrogram indicating their evolutionary relatedness. I am using the excellent ape package (http://cran.r-project.org/src/contrib/Descriptions/ape.html) to import the evolutionary dendrograms. I have already manipulated the dendrogram to be ultrameric, with branches all the same length, to prevent an error, although I would prefer not to have to do so: library(ape) mytree - read.tree(file = ultra.newick, text = NULL, tree.names = NULL, skip = 0, comment.char = #) #I then change them into a hclust: tree - as.hclust(mytree) #and make this into a dendrogram dend - as.dendrogram(tree) However, when I use this dendrogram as part of the heatmap, the relationships in the dendrogram that I'm loading are not preserved, although the order of bacteria in the heatmap changes: hv - heatmap(x, col = cm.colors(256), scale=none, Rowv=dend, keep.dendro=TRUE, xlab = Environments, ylab= Species, main = heatmap of species present in environments) Is there something obvious I am missing? When I plot the hclust and dendrogram, they seem to preserve the relationships that I am attempting to show, but not when the dendrogram is used in the heatmap. I'm not sure I really understand the datatype of R dendrogram objects, and this may be the source of my problems. The heatmap documentation (http://addictedtor.free.fr/graphiques/help/index.html?pack=statsalias=h eatmapfun=heatmap.html) is somewhat opaque to someone with my programing skills. Would it be better to reorder the heatmap and then somehow add in the dendrogram with add.expr command? If anyone could suggest something, or point me in the right direction I would greatly appreciate it. jamesh Here are the contents of the two files used as examples above: env.dat: soil1 soil2 soil3 marine1 marine2 marine3 One 23 24 26 0 0 0 Two 43 43 43 3 6 8 Three 56 78 45 23 45 56 Four 6 6 3 34 56 34 Five 2 17 12 76 45 65 ultra.newick: (((One:1,Four:1):1,(Three:1,Two:1):1):1,Five:3):0.0; __ R-help@stat.math.ethz.ch mailing list 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] Centering a legend
Using layout I am plotting 5 boxplots on top of each other, all of them using colored boxes in different order. In a 6th box below I want to give the legend explaining the box colors, and I want this box to be centered in an otherwise empty plot. I am creating this plot by plot(0:1,0:1,type=n) Is there an easy way to find the parameters x and y for calling legend which will center the legend in this plot? -- Erich Neuwirth, University of Vienna Faculty of Computer Science Computer Supported Didactics Working Group Visit our SunSITE at http://sunsite.univie.ac.at Phone: +43-1-4277-39464 Fax: +43-1-4277-39459 __ R-help@stat.math.ethz.ch mailing list 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] graphing with barchart question
I am graphing data using barchart (barchart(DV ~ IV | subject). I have 2 groups of 9 subjects each. How can I easily identify which group each subject belongs to? I have been trying to color code them, but can't seem to get that to work. Any suggestions would be appreciated. Thanks, Matt Bridgman __ R-help@stat.math.ethz.ch mailing list 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.
Re: [R] Centering a legend
On Mon, 2007-05-07 at 19:44 +0200, Erich Neuwirth wrote: Using layout I am plotting 5 boxplots on top of each other, all of them using colored boxes in different order. In a 6th box below I want to give the legend explaining the box colors, and I want this box to be centered in an otherwise empty plot. I am creating this plot by plot(0:1,0:1,type=n) Is there an easy way to find the parameters x and y for calling legend which will center the legend in this plot? Eric, legend() has positional arguments to enable you to do what you require. Something along the lines of the following: # Do your plot plot(0:1, 0:1, type = n, ann = FALSE, axes = FALSE) # Do the legend, centered legend(center, This is a test legend) See the Details section of ?legend, third paragraph. Also, the last set of examples there. HTH, Marc Schwartz __ R-help@stat.math.ethz.ch mailing list 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.
Re: [R] graphing with barchart question
On 5/7/07, Matthew Bridgman [EMAIL PROTECTED] wrote: I am graphing data using barchart (barchart(DV ~ IV | subject). I have 2 groups of 9 subjects each. How can I easily identify which group each subject belongs to? I have been trying to color code them, but can't seem to get that to work. Any suggestions would be appreciated. Not sure what you are looking for (is your group one of DV and IV, or is it something entirely different?). A reproducible example would have helped, but try looking at ?panel.barchart and use the 'groups' and 'stack' arguments. -Deepayan __ R-help@stat.math.ethz.ch mailing list 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] like apply(x,1,sum), but using multiplication?
Hi, I need to multiply all columns in a matrix so something like apply(x,2,sum), but using multiplication should do. I have tried apply(x,2,*) I know this must be trivial, but I get: Error in FUN(newX[, i], ...) : invalid unary operator The help for apply states that unary operators must be quoted. I tried single quotes too, with the same results. Thanks, -Jose -- Jose Quesada, PhD. http://www.andrew.cmu.edu/~jquesada __ R-help@stat.math.ethz.ch mailing list 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.
Re: [R] like apply(x,1,sum), but using multiplication?
Try: apply(x, 2, prod) -- Henrique Dallazuanna Curitiba-Paraná-Brasil 25° 25' 40 S 49° 16' 22 Ohttp://maps.google.com/maps?f=qhl=enq=Curitiba,+Brazillayer=ie=UTF8z=18ll=-25.448315,-49.276916spn=0.002054,0.005407t=kom=1 On 07/05/07, Jose Quesada [EMAIL PROTECTED] wrote: Hi, I need to multiply all columns in a matrix so something like apply(x,2,sum), but using multiplication should do. I have tried apply(x,2,*) I know this must be trivial, but I get: Error in FUN(newX[, i], ...) : invalid unary operator The help for apply states that unary operators must be quoted. I tried single quotes too, with the same results. Thanks, -Jose -- Jose Quesada, PhD. http://www.andrew.cmu.edu/~jquesada __ R-help@stat.math.ethz.ch mailing list 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. [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list 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.
Re: [R] like apply(x,1,sum), but using multiplication?
see ?prod b On May 7, 2007, at 2:25 PM, Jose Quesada wrote: Hi, I need to multiply all columns in a matrix so something like apply(x,2,sum), but using multiplication should do. I have tried apply(x,2,*) I know this must be trivial, but I get: Error in FUN(newX[, i], ...) : invalid unary operator The help for apply states that unary operators must be quoted. I tried single quotes too, with the same results. Thanks, -Jose -- Jose Quesada, PhD. http://www.andrew.cmu.edu/~jquesada __ R-help@stat.math.ethz.ch mailing list 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@stat.math.ethz.ch mailing list 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.
Re: [R] like apply(x,1,sum), but using multiplication?
Jose Quesada said the following on 5/7/2007 11:25 AM: Hi, I need to multiply all columns in a matrix so something like apply(x,2,sum), but using multiplication should do. I have tried apply(x,2,*) I know this must be trivial, but I get: Error in FUN(newX[, i], ...) : invalid unary operator The help for apply states that unary operators must be quoted. I tried single quotes too, with the same results. Thanks, -Jose Try: apply(x,2,prod) HTH, --sundar __ R-help@stat.math.ethz.ch mailing list 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] TukeyHSD fails on my data
Howdo folks, So I have my data (attached). There are two columns I'm interested in; algname and dur. I'd like to know how dur changes with algname. algname is nominal and there are 7 possibilities. There are two more nominal independents, task and id, so my model is: dur ~ algname+task+id From the research I've done, a TukeyHSD seems to be what I need, so I do: TukeyHSD(aov(dur ~ algname+task+id, cstats), ordered=T) But I get this back: Error in rep.int(n, length(means)) : unimplemented type 'NULL' in 'rep' In addition: Warning messages: 1: non-factors ignored: task in: replications(paste(~, xx), data = mf) 2: non-factors ignored: id in: replications(paste(~, xx), data = mf) Can anyone shed any light on the situation? Gav __ R-help@stat.math.ethz.ch mailing list 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] H-scatter plot in geostatistics
Hello: I would like to make h-scatter plot. My data looks like as follow: x, y, z, 12.0, 11.2, 12, 10.21, 5.42, 8, 5.12, -8.25, 7, I want to make h-scatter plot for the z values by difference distance h from 1 to 20. There are 1023 observations. Do I need to change data class from data-frame to others? I am not sure we can make h-scatter plot in R. Have a nice day. HongSU. __ R-help@stat.math.ethz.ch mailing list 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.
Re: [R] summing values according to a factor
Dear Andy and all others who have replied: On 07/05/07, Liaw, Andy [EMAIL PROTECTED] wrote: I guess what you want to do is compare Q1/T1 among the sections? If you want to compute the sum of Q1/T1 by Section, you can do something like: sum.by.section - with(mydata, tapply(Q1/T1, section, sum)) Substitute sum with anything you want to compute. That worked perfectly. Thanks! Enrico -- Enrico Indiogine Mathematics Education Texas AM University [EMAIL PROTECTED] __ R-help@stat.math.ethz.ch mailing list 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.
Re: [R] graphing with barchart question
Sorry. I have attached my data frame: DV = dv; IV = bins; subject = id, Group = group. barchart(dv ~ bins | id + group, groups = group, data = matt.df) The two suggestions you offered give me error messages regarding invalid line type and do not plot all of the data. If I drop the 'groups' argument I get all of the data, but it plots all combinations of group and id (yielding 36 plots instead of 18). Is there a way to eliminate some of the combinations? Thanks again for your help. Matt It doesn't have to be, you can do barchart(DV ~ IV | subject + Group, groups = Group) or barchart(DV ~ IV | subject:Group, groups = Group, auto.key = TRUE) Of course, I can't check this because you still haven't given me a reproducible example (and please CC r-help so that follow ups remain archived). -Deepayan __ R-help@stat.math.ethz.ch mailing list 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.
Re: [R] TukeyHSD fails on my data
We don't have the data (nothing useful was attached - see the posting guide for what you can attach), but it looks like you should be using the 'which' argument to TukeyHSD. On Mon, 7 May 2007, Gav Wood wrote: Howdo folks, So I have my data (attached). There are two columns I'm interested in; algname and dur. I'd like to know how dur changes with algname. algname is nominal and there are 7 possibilities. There are two more nominal independents, task and id, so my model is: dur ~ algname+task+id From the research I've done, a TukeyHSD seems to be what I need, so I do: TukeyHSD(aov(dur ~ algname+task+id, cstats), ordered=T) But I get this back: Error in rep.int(n, length(means)) : unimplemented type 'NULL' in 'rep' In addition: Warning messages: 1: non-factors ignored: task in: replications(paste(~, xx), data = mf) 2: non-factors ignored: id in: replications(paste(~, xx), data = mf) Can anyone shed any light on the situation? Gav -- Brian D. Ripley, [EMAIL PROTECTED] Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UKFax: +44 1865 272595 __ R-help@stat.math.ethz.ch mailing list 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.
Re: [R] A function for raising a matrix to a power?
Michael Dewey wrote: I may be revealing my ignorance here, but is MatrixExp in the msm package (available from CRAN) not relevant here? Never heard about it. Searching with help.search(matrix) didn't show any function that might be used as Matrix^n. Alberto Monteiro __ R-help@stat.math.ethz.ch mailing list 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] logrank test
hie how do you compute the logrank test using R what commands do you use thanks - Don't pick lemons. [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list 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] computing logrank statistic/test
hie how do you compute the logrank test using R what commands do you use my data looks something like just an example treatmentgrp strata censoringTime survivalTime censoring act.surv.time [1,] 2 2 42.89005 1847.3358 1 42.89005 [2,] 1 1 74.40379 440.3467 1 74.40379 [3,] 2 2 35.17913 344.1113 1 35.17913 [4,] 2 2 55.79590 741.2375 1 55.79590 [5,] 2 1 384.73976 189.7997 0 189.79968 [6,] 1 2 112.76215 375.8227 1 112.76215 [7,] 2 2 20.82441 658.0410 1 20.82441 [8,] 2 2 30.58497 1537.5776 1 30.58497 [9,] 2 1 128.48140 306.9908 1 128.48140 l would like to compare the result with SAS - The fish are biting. [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list 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.
Re: [R] TukeyHSD fails on my data
Prof Brian Ripley ripley at stats.ox.ac.uk writes: We don't have the data (nothing useful was attached - see the posting guide for what you can attach), but it looks like you should be using the 'which' argument to TukeyHSD. But I get this back: Error in rep.int(n, length(means)) : unimplemented type 'NULL' in 'rep' In addition: Warning messages: 1: non-factors ignored: task in: replications(paste(~, xx), data = mf) 2: non-factors ignored: id in: replications(paste(~, xx), data = mf) Can anyone shed any light on the situation? The problem was, as Jeff Laake adroitly pointed out, I had failed to tell R that the latter two independents were factors. It was easliy fixed with as.factor(). Thanks folks, Gav __ R-help@stat.math.ethz.ch mailing list 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.
Re: [R] A function for raising a matrix to a power?
On Mon, 7 May 2007, Alberto Vieira Ferreira Monteiro wrote: Michael Dewey wrote: I may be revealing my ignorance here, but is MatrixExp in the msm package (available from CRAN) not relevant here? Never heard about it. Searching with help.search(matrix) didn't show any function that might be used as Matrix^n. That only searches packages you have installed: MatrixExp certainly shows up on my system. RSiteSearch is a better way to find out about packages that you may not have installed: it shows several possibilities. -- Brian D. Ripley, [EMAIL PROTECTED] Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UKFax: +44 1865 272595 __ R-help@stat.math.ethz.ch mailing list 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.
Re: [R] computing logrank statistic/test
On Mon, 2007-05-07 at 14:02 -0700, raymond chiruka wrote: hie how do you compute the logrank test using R what commands do you use my data looks something like just an example treatmentgrp strata censoringTime survivalTime censoring act.surv.time [1,] 2 2 42.89005 1847.3358 1 42.89005 [2,] 1 1 74.40379 440.3467 1 74.40379 [3,] 2 2 35.17913 344.1113 1 35.17913 [4,] 2 2 55.79590 741.2375 1 55.79590 [5,] 2 1 384.73976 189.7997 0 189.79968 [6,] 1 2 112.76215 375.8227 1 112.76215 [7,] 2 2 20.82441 658.0410 1 20.82441 [8,] 2 2 30.58497 1537.5776 1 30.58497 [9,] 2 1 128.48140 306.9908 1 128.48140 l would like to compare the result with SAS Raymond, You asked the same question back on May 1, with several answers by myself, Peter Dalgaard and Roland Rau. Did you not see the replies? https://stat.ethz.ch/pipermail/r-help/2007-May/130956.html HTH, Marc Schwartz __ R-help@stat.math.ethz.ch mailing list 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.
Re: [R] computing logrank statistic/test
2007/5/7, raymond chiruka [EMAIL PROTECTED]: hie how do you compute the logrank test using R what commands do you use my data looks something like just an example treatmentgrp strata censoringTime survivalTime censoring act.surv.time [1,] 2 2 42.89005 1847.3358 1 42.89005 [2,] 1 1 74.40379 440.3467 1 74.40379 [3,] 2 2 35.17913 344.1113 1 35.17913 [4,] 2 2 55.79590 741.2375 1 55.79590 [5,] 2 1 384.73976 189.7997 0 189.79968 [6,] 1 2 112.76215 375.8227 1 112.76215 [7,] 2 2 20.82441 658.0410 1 20.82441 [8,] 2 2 30.58497 1537.5776 1 30.58497 [9,] 2 1 128.48140 306.9908 1 128.48140 l would like to compare the result with SAS - The fish are biting. [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list 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. Please, check this page: http://www.ats.ucla.edu/stat/R/examples/asa/asa_ch2_r.htm Rod. __ R-help@stat.math.ethz.ch mailing list 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.
Re: [R] A function for raising a matrix to a power?
You might try this, from 9/22/2006 with subject line Exponentiate a matrix:
I am getting a bit rusty on some of these things, but I seem to recall
that there is a numerical advantage (speed and/or accuracy?) to
diagonalizing:
expM - function(X,e) { v - La.svd(X); v$u %*% diag(v$d^e) %*% v$vt }
P - matrix(c(.3,.7, .7, .3), ncol=2)
P %*% P %*% P
[,1] [,2]
[1,] 0.468 0.532
[2,] 0.532 0.468
expM(P,3)
[,1] [,2]
[1,] 0.468 0.532
[2,] 0.532 0.468
I think this also works for non-integer, negative, large, and complex
exponents:
expM(P, 1.5)
[,1] [,2]
[1,] 0.3735089 0.6264911
[2,] 0.6264911 0.3735089
expM(P, 1000)
[,1] [,2]
[1,] 0.5 0.5
[2,] 0.5 0.5
expM(P, -3)
[,1][,2]
[1,] -7.3125 8.3125
[2,] 8.3125 -7.3125
expM(P, 3+.5i)
[,1] [,2]
[1,] 0.4713+0.0141531i 0.5287-0.0141531i
[2,] 0.5287-0.0141531i 0.4713+0.0141531i
Paul Gilbert
Doran, Harold wrote:
Suppose I have a square matrix P
P - matrix(c(.3,.7, .7, .3), ncol=2)
I know that
P * P
Returns the element by element product, whereas
P%*%P
Returns the matrix product.
Now, P2 also returns the element by element product. But, is there a
slick way to write
P %*% P %*% P
Obviously, P3 does not return the result I expect.
Thanks,
Harold
Atte Tenkanen wrote:
Hi,
Is there a function for raising a matrix to a power? For example if you like
to compute A%*%A%*%A, is there an abbreviation similar to A^3?
Atte Tenkanen
A=rbind(c(1,1),c(-1,-2))
A
[,1] [,2]
[1,]11
[2,] -1 -2
A^3
[,1] [,2]
[1,]11
[2,] -1 -8
But:
A%*%A%*%A
[,1] [,2]
[1,]12
[2,] -2 -5
__
R-help@stat.math.ethz.ch mailing list
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.
La version française suit le texte anglais.
This email may contain privileged and/or confidential inform...{{dropped}}
__
R-help@stat.math.ethz.ch mailing list
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] getting informative error messages
Certain errors seem to generate messages that are less informative than
most -- they just tell you which function an error happened in, but
don't indicate which line or expression the error occurred in.
Here's a toy example:
f - function(x) {a - 1; y - x[list(1:3)]; b - 2; return(y)}
options(error=NULL)
f(1:3)
Error in f(1:3) : invalid subscript type
traceback()
1: f(1:3)
In this function, it's clear that the error is in subscripting 'x', but
it's not always so immediately obvious in lengthier functions.
Is there anything I can do to get a more informative error message in
this type of situation? I couldn't find any help in the section
Debugging R Code in R-exts (or anything at all relevant in R-intro).
(Different values for options(error=...) and different formatting of the
function made no difference.)
-- Tony Plate
sessionInfo()
R version 2.5.0 (2007-04-23)
i386-pc-mingw32
locale:
LC_COLLATE=English_United States.1252;LC_CTYPE=English_United
States.1252;LC_MONETARY=English_United
States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1252
attached base packages:
[1] stats graphics grDevices utils datasets methods
[7] base
other attached packages:
tap.misc
1.0
__
R-help@stat.math.ethz.ch mailing list
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] Bad optimization solution
Dear All
I am trying to perform the below optimization problem, but getting
(0.5,0.5) as optimal solution, which is wrong; the correct solution
should be (1,0) or (0,1).
Am I doing something wrong? I am using R 2.5.0 on Fedora Core 6 (Linux).
Thanks in advance,
Paul
--
myfunc - function(x) {
x1 - x[1]
x2 - x[2]
abs(x1-x2)
}
optim(c(0.5,0.5),myfunc,lower=c(0,0),upper=c(1,1),method=L-BFGS-B,control=list(fnscale=-1))
__
R-help@stat.math.ethz.ch mailing list
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.
Re: [R] Bad optimization solution
Paul,
I think the problem is the starting point. I do not remember the details
of the BFGS method, but I am almost sure the (.5, .5) starting point is
suspect, since the abs function is not differentiable at 0. If you perturb
the starting point even slightly you will have no problem.
Andy
__
Andy Jaworski
518-1-01
Process Laboratory
3M Corporate Research Laboratory
-
E-mail: [EMAIL PROTECTED]
Tel: (651) 733-6092
Fax: (651) 736-3122
Paul Smith
[EMAIL PROTECTED]
To
Sent by: R-help r-help@stat.math.ethz.ch
[EMAIL PROTECTED] cc
at.math.ethz.ch
Subject
[R] Bad optimization solution
05/07/2007 04:30
PM
Dear All
I am trying to perform the below optimization problem, but getting
(0.5,0.5) as optimal solution, which is wrong; the correct solution
should be (1,0) or (0,1).
Am I doing something wrong? I am using R 2.5.0 on Fedora Core 6 (Linux).
Thanks in advance,
Paul
--
myfunc - function(x) {
x1 - x[1]
x2 - x[2]
abs(x1-x2)
}
optim(c(0.5,0.5),myfunc,lower=c(0,0),upper=c(1,1),method=L-BFGS-B,control=list(fnscale=-1))
__
R-help@stat.math.ethz.ch mailing list
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@stat.math.ethz.ch mailing list
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.
Re: [R] A function for raising a matrix to a power?
Paul,
Your solution based on SVD does not work even for the matrix in your example
(the reason it worked for e=3, was that it is an odd power and since P is a
permutation matrix. It will also work for all odd powers, but not for even
powers).
However, a spectral decomposition will work for all symmetric matrices (SVD
based exponentiation doesn't work for any matrix).
Here is the function to do exponentiation based on spectral decomposition:
expM.sd - function(X,e){Xsd - eigen(X); Xsd$vec %*% diag(Xsd$val^e) %*%
t(Xsd$vec)}
P - matrix(c(.3,.7, .7, .3), ncol=2)
P
[,1] [,2]
[1,] 0.3 0.7
[2,] 0.7 0.3
P %*% P
[,1] [,2]
[1,] 0.58 0.42
[2,] 0.42 0.58
expM(P,2)
[,1] [,2]
[1,] 0.42 0.58
[2,] 0.58 0.42
expM.sd(P,2)
[,1] [,2]
[1,] 0.58 0.42
[2,] 0.42 0.58
Exponentiation based on spectral decomposition should be generally more
accurate (not sure about the speed). A caveat is that it will only work for
real, symmetric or complex, hermitian matrices. It will not work for
asymmetric matrices. It would work for asymmetric matrix if the
eigenvectors are made to be orthonormal.
Ravi.
---
Ravi Varadhan, Ph.D.
Assistant Professor, The Center on Aging and Health
Division of Geriatric Medicine and Gerontology
Johns Hopkins University
Ph: (410) 502-2619
Fax: (410) 614-9625
Email: [EMAIL PROTECTED]
Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html
-Original Message-
From: [EMAIL PROTECTED]
[mailto:[EMAIL PROTECTED] On Behalf Of Paul Gilbert
Sent: Monday, May 07, 2007 5:16 PM
To: Atte Tenkanen
Cc: r-help@stat.math.ethz.ch
Subject: Re: [R] A function for raising a matrix to a power?
You might try this, from 9/22/2006 with subject line Exponentiate a matrix:
I am getting a bit rusty on some of these things, but I seem to recall
that there is a numerical advantage (speed and/or accuracy?) to
diagonalizing:
expM - function(X,e) { v - La.svd(X); v$u %*% diag(v$d^e) %*% v$vt }
P - matrix(c(.3,.7, .7, .3), ncol=2)
P %*% P %*% P
[,1] [,2]
[1,] 0.468 0.532
[2,] 0.532 0.468
expM(P,3)
[,1] [,2]
[1,] 0.468 0.532
[2,] 0.532 0.468
I think this also works for non-integer, negative, large, and complex
exponents:
expM(P, 1.5)
[,1] [,2]
[1,] 0.3735089 0.6264911
[2,] 0.6264911 0.3735089
expM(P, 1000)
[,1] [,2]
[1,] 0.5 0.5
[2,] 0.5 0.5
expM(P, -3)
[,1][,2]
[1,] -7.3125 8.3125
[2,] 8.3125 -7.3125
expM(P, 3+.5i)
[,1] [,2]
[1,] 0.4713+0.0141531i 0.5287-0.0141531i
[2,] 0.5287-0.0141531i 0.4713+0.0141531i
Paul Gilbert
Doran, Harold wrote:
Suppose I have a square matrix P
P - matrix(c(.3,.7, .7, .3), ncol=2)
I know that
P * P
Returns the element by element product, whereas
P%*%P
Returns the matrix product.
Now, P2 also returns the element by element product. But, is there a
slick way to write
P %*% P %*% P
Obviously, P3 does not return the result I expect.
Thanks,
Harold
Atte Tenkanen wrote:
Hi,
Is there a function for raising a matrix to a power? For example if you
like to compute A%*%A%*%A, is there an abbreviation similar to A^3?
Atte Tenkanen
A=rbind(c(1,1),c(-1,-2))
A
[,1] [,2]
[1,]11
[2,] -1 -2
A^3
[,1] [,2]
[1,]11
[2,] -1 -8
But:
A%*%A%*%A
[,1] [,2]
[1,]12
[2,] -2 -5
__
R-help@stat.math.ethz.ch mailing list
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.
La version française suit le texte anglais.
This email may contain privileged and/or confidential inform...{{dropped}}
__
R-help@stat.math.ethz.ch mailing list
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@stat.math.ethz.ch mailing list
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.
Re: [R] graphing with barchart question
On 5/7/07, Matthew Bridgman [EMAIL PROTECTED] wrote: Sorry. I have attached my data frame: DV = dv; IV = bins; subject = id, Group = group. barchart(dv ~ bins | id + group, groups = group, data = matt.df) The two suggestions you offered give me error messages regarding invalid line type and do not plot all of the data. If I drop the 'groups' argument I get all of the data, but it plots all combinations of group and id (yielding 36 plots instead of 18). Is there a way to eliminate some of the combinations? This seems to work for me: barchart(dv ~ factor(bins) | interaction(id, group), data = matt.df, groups = group, origin = 0, auto.key = list(columns = 2)) If you get the invalid line type thing, try upgrading R and lattice. -Deepayan __ R-help@stat.math.ethz.ch mailing list 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.
Re: [R] Bad optimization solution
On 5/7/07, [EMAIL PROTECTED] [EMAIL PROTECTED] wrote:
I think the problem is the starting point. I do not remember the details
of the BFGS method, but I am almost sure the (.5, .5) starting point is
suspect, since the abs function is not differentiable at 0. If you perturb
the starting point even slightly you will have no problem.
Paul Smith
[EMAIL PROTECTED]
To
Sent by: R-help r-help@stat.math.ethz.ch
[EMAIL PROTECTED] cc
at.math.ethz.ch
Subject
[R] Bad optimization solution
05/07/2007 04:30
PM
Dear All
I am trying to perform the below optimization problem, but getting
(0.5,0.5) as optimal solution, which is wrong; the correct solution
should be (1,0) or (0,1).
Am I doing something wrong? I am using R 2.5.0 on Fedora Core 6 (Linux).
Thanks in advance,
Paul
--
myfunc - function(x) {
x1 - x[1]
x2 - x[2]
abs(x1-x2)
}
optim(c(0.5,0.5),myfunc,lower=c(0,0),upper=c(1,1),method=L-BFGS-B,control=list(fnscale=-1))
Yes, with (0.2,0.9), a correct solution comes out. However, how can
one be sure in general that the solution obtained by optim is correct?
In ?optim says:
Method 'L-BFGS-B' is that of Byrd _et. al._ (1995) which allows
_box constraints_, that is each variable can be given a lower
and/or upper bound. The initial value must satisfy the
constraints. This uses a limited-memory modification of the BFGS
quasi-Newton method. If non-trivial bounds are supplied, this
method will be selected, with a warning.
which only demands that the initial value must satisfy the constraints.
Paul
__
R-help@stat.math.ethz.ch mailing list
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.
Re: [R] Bad optimization solution
On 5/7/07, Paul Smith [EMAIL PROTECTED] wrote:
I think the problem is the starting point. I do not remember the details
of the BFGS method, but I am almost sure the (.5, .5) starting point is
suspect, since the abs function is not differentiable at 0. If you perturb
the starting point even slightly you will have no problem.
Paul Smith
[EMAIL PROTECTED]
To
Sent by: R-help r-help@stat.math.ethz.ch
[EMAIL PROTECTED] cc
at.math.ethz.ch
Subject
[R] Bad optimization solution
05/07/2007 04:30
PM
Dear All
I am trying to perform the below optimization problem, but getting
(0.5,0.5) as optimal solution, which is wrong; the correct solution
should be (1,0) or (0,1).
Am I doing something wrong? I am using R 2.5.0 on Fedora Core 6 (Linux).
Thanks in advance,
Paul
--
myfunc - function(x) {
x1 - x[1]
x2 - x[2]
abs(x1-x2)
}
optim(c(0.5,0.5),myfunc,lower=c(0,0),upper=c(1,1),method=L-BFGS-B,control=list(fnscale=-1))
Yes, with (0.2,0.9), a correct solution comes out. However, how can
one be sure in general that the solution obtained by optim is correct?
In ?optim says:
Method 'L-BFGS-B' is that of Byrd _et. al._ (1995) which allows
_box constraints_, that is each variable can be given a lower
and/or upper bound. The initial value must satisfy the
constraints. This uses a limited-memory modification of the BFGS
quasi-Newton method. If non-trivial bounds are supplied, this
method will be selected, with a warning.
which only demands that the initial value must satisfy the constraints.
Furthermore, X^2 is everywhere differentiable and notwithstanding the
reported problem occurs with
myfunc - function(x) {
x1 - x[1]
x2 - x[2]
(x1-x2)^2
}
optim(c(0.2,0.2),myfunc,lower=c(0,0),upper=c(1,1),method=L-BFGS-B,control=list(fnscale=-1))
Paul
__
R-help@stat.math.ethz.ch mailing list
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.
Re: [R] Bad optimization solution
Paul Smith said the following on 5/7/2007 3:25 PM:
On 5/7/07, Paul Smith [EMAIL PROTECTED] wrote:
I think the problem is the starting point. I do not remember the details
of the BFGS method, but I am almost sure the (.5, .5) starting point is
suspect, since the abs function is not differentiable at 0. If you perturb
the starting point even slightly you will have no problem.
Paul Smith
[EMAIL PROTECTED]
To
Sent by: R-help r-help@stat.math.ethz.ch
[EMAIL PROTECTED] cc
at.math.ethz.ch
Subject
[R] Bad optimization solution
05/07/2007 04:30
PM
Dear All
I am trying to perform the below optimization problem, but getting
(0.5,0.5) as optimal solution, which is wrong; the correct solution
should be (1,0) or (0,1).
Am I doing something wrong? I am using R 2.5.0 on Fedora Core 6 (Linux).
Thanks in advance,
Paul
--
myfunc - function(x) {
x1 - x[1]
x2 - x[2]
abs(x1-x2)
}
optim(c(0.5,0.5),myfunc,lower=c(0,0),upper=c(1,1),method=L-BFGS-B,control=list(fnscale=-1))
Yes, with (0.2,0.9), a correct solution comes out. However, how can
one be sure in general that the solution obtained by optim is correct?
In ?optim says:
Method 'L-BFGS-B' is that of Byrd _et. al._ (1995) which allows
_box constraints_, that is each variable can be given a lower
and/or upper bound. The initial value must satisfy the
constraints. This uses a limited-memory modification of the BFGS
quasi-Newton method. If non-trivial bounds are supplied, this
method will be selected, with a warning.
which only demands that the initial value must satisfy the constraints.
Furthermore, X^2 is everywhere differentiable and notwithstanding the
reported problem occurs with
myfunc - function(x) {
x1 - x[1]
x2 - x[2]
(x1-x2)^2
}
optim(c(0.2,0.2),myfunc,lower=c(0,0),upper=c(1,1),method=L-BFGS-B,control=list(fnscale=-1))
Paul
Then perhaps supply the gradient:
mygrad - function(x) {
x1 - x[1]
x2 - x[2]
c(2, -2) * c(x1, x2)
}
optim(c(0.2,0.2),myfunc,mygrad,lower=c(0,0),upper=c(1,1),
method=L-BFGS-B,control=list(fnscale=-1))
HTH,
--sundar
__
R-help@stat.math.ethz.ch mailing list
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.
Re: [R] Bad optimization solution
Your function, (x1-x2)^2, has zero gradient at all the starting values such
that x1 = x2, which means that the gradient-based search methods will
terminate there because they have found a critical point, i.e. a point at
which the gradient is zero (which can be a maximum or a minimum or a saddle
point).
However, I do not why optim converges to the boundary maximum, when analytic
gradient is supplied (as shown by Sundar).
Ravi.
---
Ravi Varadhan, Ph.D.
Assistant Professor, The Center on Aging and Health
Division of Geriatric Medicine and Gerontology
Johns Hopkins University
Ph: (410) 502-2619
Fax: (410) 614-9625
Email: [EMAIL PROTECTED]
Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html
-Original Message-
From: [EMAIL PROTECTED]
[mailto:[EMAIL PROTECTED] On Behalf Of Paul Smith
Sent: Monday, May 07, 2007 6:26 PM
To: R-help
Subject: Re: [R] Bad optimization solution
On 5/7/07, Paul Smith [EMAIL PROTECTED] wrote:
I think the problem is the starting point. I do not remember the
details
of the BFGS method, but I am almost sure the (.5, .5) starting point is
suspect, since the abs function is not differentiable at 0. If you
perturb
the starting point even slightly you will have no problem.
Paul Smith
[EMAIL PROTECTED]
To
Sent by: R-help r-help@stat.math.ethz.ch
[EMAIL PROTECTED]
cc
at.math.ethz.ch
Subject
[R] Bad optimization solution
05/07/2007 04:30
PM
Dear All
I am trying to perform the below optimization problem, but getting
(0.5,0.5) as optimal solution, which is wrong; the correct solution
should be (1,0) or (0,1).
Am I doing something wrong? I am using R 2.5.0 on Fedora Core 6 (Linux).
Thanks in advance,
Paul
--
myfunc - function(x) {
x1 - x[1]
x2 - x[2]
abs(x1-x2)
}
optim(c(0.5,0.5),myfunc,lower=c(0,0),upper=c(1,1),method=L-BFGS-B,control=
list(fnscale=-1))
Yes, with (0.2,0.9), a correct solution comes out. However, how can
one be sure in general that the solution obtained by optim is correct?
In ?optim says:
Method 'L-BFGS-B' is that of Byrd _et. al._ (1995) which allows
_box constraints_, that is each variable can be given a lower
and/or upper bound. The initial value must satisfy the
constraints. This uses a limited-memory modification of the BFGS
quasi-Newton method. If non-trivial bounds are supplied, this
method will be selected, with a warning.
which only demands that the initial value must satisfy the constraints.
Furthermore, X^2 is everywhere differentiable and notwithstanding the
reported problem occurs with
myfunc - function(x) {
x1 - x[1]
x2 - x[2]
(x1-x2)^2
}
optim(c(0.2,0.2),myfunc,lower=c(0,0),upper=c(1,1),method=L-BFGS-B,control=
list(fnscale=-1))
Paul
__
R-help@stat.math.ethz.ch mailing list
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@stat.math.ethz.ch mailing list
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] Analyzing Stacked Time Series
I have a question about pooling or stacking several time series
“samples” (sorry in advance for the long, possibly confusing, message).
I'm sure I'm revealing far more ignorance than I'm aware of, but
that's why I'm sending this...
[Example at bottom]
I have regional migration flows (“samples”) from, say, regions A to B, A
to C, B to A, …., C to B (Noted as xIJ in my example). Each of these is
a time series, but I want to use R’s time series tools to fit an ARIMA
(or ARMA) model along with some independent variables (Noted as YIJ,
ZIJ, where Y and Z are characteristics (e.g., per capita GDP), I
indicates the region, and J in {O, D} indicates whether it is the origin
or destination region).
I can do this easily enough for a single flow, say from A to B. But
what if I want to find one set of coefficients for all of the flows?
While I can construct vectors that look like what I think they should, I
can’t figure out how to coerce them into being “stacked” time series,
with the time indices repeating through each vector
(t=1,...t=10,t=1,...,t=10,...,t=1,...t=10). And since I can’t even
figure out how to do this (or whether it is possible/sensible), I
certainly can’t use R’s time series packages on this “stacked” time series.
For my example (below), it’s easy enough to do something like,
fit1 - lm(X ~ YO + YD + ZO + ZD)
But anything involving lags is not correct, since it seems—at best—to be
treating my “stacked” time series (i.e., 6 series of length 10) as one
series of length 60. For example, these give questionable, if any, output.
arima1 - arima(X, order = c(1, 0, 0))
fit.gls001 - gls(X ~ YO + YD + ZO + ZD,
correlation = corARMA(p = 2), method = ML)
fit.gls002 - gls(X ~ YO + YD + ZO + ZD +
lag(YO) + lag(YD) + lag(ZO) + lag(ZD),
correlation = corARMA(p = 1), method = ML)
ar001 - ar.ols(cbin(X, YO, YD, ZO, ZD))
Here is my example:
xAB - as.ts(round(rnorm(10, 0, 1), 2), start = c(1990, 1))
xAC - as.ts(round(rnorm(10, 0, 1), 2), start = c(1990, 1))
xBA - as.ts(round(rnorm(10, .75, 1.5), 2), start = c(1990, 1))
xBC - as.ts(round(rnorm(10, .25, 1.9), 2), start = c(1990, 1))
xCA - as.ts(round(rnorm(10, 1.5, 2.2), 2), start = c(1990, 1))
xCB - as.ts(round(rnorm(10, .5, 0.8), 2), start = c(1990, 1))
yA - as.ts(round(rnorm(10, -0.2, 0.3), 2), start = c(1990, 1))
yB - as.ts(round(runif(10, -1, 1), 2), start = c(1990, 1))
yC - as.ts(round(runif(10, -1, 1), 2), start = c(1990, 1))
zA - as.ts(round(rnorm(10, -1.5, 2), 2), start = c(1990, 1))
zB - as.ts(round(rnorm(10, 0, 0.5), 2), start = c(1990, 1))
zC - as.ts(round(runif(10, 0, 0.5), 2), start = c(1990, 1))
Orig - c(1, 1, 2, 2, 3, 3)
Dest - c(2, 3, 1, 3, 1, 2)
Xt - cbind(xAB, xAC, xBA, xBC, xCA, xCB)
Yt - cbind(yA, yB, yC)
Zt - cbind(zA, zB, zC)
X - vector()
for(i in 1:ncol(Xt)){
X - append(X, Xt[,i])
}
YO - vector()
for(i in 1:length(Orig)){
YO - append(YO, Yt[,Orig[i]])
}
ZO - vector()
for(i in 1:length(Orig)){
ZO - append(ZO, Zt[,Orig[i]])
}
YD - vector()
for(i in 1:length(Dest)){
YD - append(YD, Yt[,Dest[i]])
}
ZD - vector()
for(i in 1:length(Dest)){
ZD - append(ZD, Zt[,Dest[i]])
}
Many thanks is advance!
Misha
__
R-help@stat.math.ethz.ch mailing list
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.
Re: [R] A function for raising a matrix to a power?
Paul Gilbert wrote: I am getting a bit rusty on some of these things, but I seem to recall that there is a numerical advantage (speed and/or accuracy?) to diagonalizing: (...) I think this also works for non-integer, negative, large, and complex This is diverging into mathematics, maybe this is off-topic, but... Not all matrices can de diagonalized, but they can be transformed (by invertible matrices) into a canonical form, called Jordan. This canonical form has the eigenvalues in the main diagonal (they may be complex numbers), 1s or 0s in the diagonal just above the main diagonal, and 0 everywhere else. Based on this, it's possible to define f(M) for _any_ function f that has enough derivatives. f(J), for a matrix in the Jordan canonical form, has f(x) in the main diagonal (where x are the eigenvalues) and f'(x) in some values of the diagonal just above that, f''(x)/2! in the next, etc. It's in Jordan normal form, in the borg of all wisdom, the wikipedia. Alberto Monteiro __ R-help@stat.math.ethz.ch mailing list 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] How might I -remove- a tree from a random forest?
I see the function getTree, which is very interesting. As I'm trying to teach myself more and more about R, and dealing with lists, it occurred to me that it might be fun to remove (as in delete) a single tree from a forest...say to go from 500 to 499. I know, I know... why? Why, to play, of course! I've been doing a lot of reading on various tuning parameters, and wanted to play with thinning as a demonstration project. Problem is, each tree is a different length (naturally), and when I try to set a portion (identified by examining the getTree function) to NULL, I get an error about lengths not matching up. Wondering if there's a tool somewhere out there to do this in a straightforward fashion, or any breadcrumbs to get me going. If this is a new thing, I'll post the function when I'm done for others to play with it, too. -- --- David L. Van Brunt, Ph.D. mailto:[EMAIL PROTECTED] If Tyranny and Oppression come to this land, it will be in the guise of fighting a foreign enemy. --James Madison [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list 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] ordered logistic regression with random effects. Howto?
I'd like to estimate an ordinal logistic regression with a random effect for a grouping variable. I do not find a pre-packaged algorithm for this. I've found methods glmmML (package: glmmML) and lmer (package: lme4) both work fine with dichotomous dependent variables. I'd like a model similar to polr (package: MASS) or lrm (package: Design) that allows random effects. I was thinking there might be a trick that might allow me to use a program written for a dichotomous dependent variable with a mixed effect to estimate such a model. The proportional odds logistic regression is often written as a sequence of dichotomous comparisons. But it seems to me that, if it would work, then somebody would have proposed it already. I've found some commentary about methods of fitting ordinal logistic regression with other procedures, however, and I would like to ask for your advice and experience with this. In this article, Ching-Fan Sheu, Fitting mixed-effects models for repeated ordinal outcomes with the NLMIXED procedure Behavior Research Methods, Instruments, Computers, 2002, 34(2): 151-157. the other gives an approach that works in SAS's NLMIXED procedure. In this approach, one explicitly writes down the probability that each level will be achieved (using the linear predictor and constants for each level). I THINK I could find a way to translate this approach into a model that can be fitted with either nlme or lmer. Has someone done it? What other strategies for fitting mixed ordinal models are there in R? Finally, a definitional question. Is a multi-category logistic regression (either ordered or not) a member of the glm family? I had thought the answer is no, mainly because glm and other R functions for glms never mention multi-category qualitative dependent variables and also because the distribution does not seem to fall into the exponential family. However, some textbooks do include the multicategory models in the GLM treatment. -- Paul E. Johnson Professor, Political Science 1541 Lilac Lane, Room 504 University of Kansas __ R-help@stat.math.ethz.ch mailing list 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.
Re: [R] ordered logistic regression with random effects. Howto?
Paul-- I think the options are pretty limited for mixed-effects ordinal regression; it might be worth taking a look at Laura Thompson's R/Splus companion to Alan Agresti's text on categorical data analysis: https://home.comcast.net/~lthompson221/Splusdiscrete2.pdf She discusses some options for both GEE and random-effects approaches, though for the ordinal mixed-effects regression, I believe she writes out the likelihood function and passes it to optim() (ie, no canned functions). Hope that helps. cheers, Dave -- Dave Atkins, PhD Assistant Professor in Clinical Psychology Fuller Graduate School of Psychology Email: [EMAIL PROTECTED] Paul wrote: I'd like to estimate an ordinal logistic regression with a random effect for a grouping variable. I do not find a pre-packaged algorithm for this. I've found methods glmmML (package: glmmML) and lmer (package: lme4) both work fine with dichotomous dependent variables. I'd like a model similar to polr (package: MASS) or lrm (package: Design) that allows random effects. I was thinking there might be a trick that might allow me to use a program written for a dichotomous dependent variable with a mixed effect to estimate such a model. The proportional odds logistic regression is often written as a sequence of dichotomous comparisons. But it seems to me that, if it would work, then somebody would have proposed it already. I've found some commentary about methods of fitting ordinal logistic regression with other procedures, however, and I would like to ask for your advice and experience with this. In this article, Ching-Fan Sheu, Fitting mixed-effects models for repeated ordinal outcomes with the NLMIXED procedure Behavior Research Methods, Instruments, Computers, 2002, 34(2): 151-157. the other gives an approach that works in SAS's NLMIXED procedure. In this approach, one explicitly writes down the probability that each level will be achieved (using the linear predictor and constants for each level). I THINK I could find a way to translate this approach into a model that can be fitted with either nlme or lmer. Has someone done it? What other strategies for fitting mixed ordinal models are there in R? Finally, a definitional question. Is a multi-category logistic regression (either ordered or not) a member of the glm family? I had thought the answer is no, mainly because glm and other R functions for glms never mention multi-category qualitative dependent variables and also because the distribution does not seem to fall into the exponential family. However, some textbooks do include the multicategory models in the GLM treatment. -- Paul E. Johnson Professor, Political Science 1541 Lilac Lane, Room 504 University of Kansas __ R-help@stat.math.ethz.ch mailing list 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.
Re: [R] Bad optimization solution
G'day Paul, On Mon, 7 May 2007 22:30:32 +0100 Paul Smith [EMAIL PROTECTED] wrote: I am trying to perform the below optimization problem, but getting (0.5,0.5) as optimal solution, which is wrong; Why? As far as I can tell you are trying to minimize |x1-x2| where x1 and x2 are between 0 and 1. The minimal value that the absolute function can take is zero and any point (x1,x2)=(x,1-x) where x is between 0 and 1 will achieve this value and also respect the constraints that you have imposed. Hence, any such point, including (0.5,0.5) is a solution to your problem. the correct solution should be (1,0) or (0,1). Why? Unless there are some additional constraint that you have not told optim() (and us) about, these are two possible solutions from an infinite set of solutions. As I said, any point of the form (x, 1-x) with x between 0 and 1 is a solution to your problem, unless I am missing something Cheers, Berwin __ R-help@stat.math.ethz.ch mailing list 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.
Re: [R] Bad optimization solution
G'day Paul,
On Mon, 7 May 2007 23:25:52 +0100
Paul Smith [EMAIL PROTECTED] wrote:
[...]
Furthermore, X^2 is everywhere differentiable and notwithstanding the
reported problem occurs with
myfunc - function(x) {
x1 - x[1]
x2 - x[2]
(x1-x2)^2
}
Same argument as with abs(x1-x2) holds. (x1-x2)^2 is non-negative for
all x1, x2. All points of the form (x, 1-x) where x is between 0 and 1
satisfy the constraints and achieve a function value of 0. Hence, all
such points are solutions.
There is no problem. Except if there are further constraints that
reduce the set of possible solutions which you have not told us about.
Cheers,
Berwin
__
R-help@stat.math.ethz.ch mailing list
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.
Re: [R] Bad optimization solution
G'day all, On Tue, 8 May 2007 12:10:25 +0800 Berwin A Turlach [EMAIL PROTECTED] wrote: I am trying to perform the below optimization problem, but getting (0.5,0.5) as optimal solution, which is wrong; Why? As far as I can tell you are trying to minimize |x1-x2| It was pointed out to me that you are actually trying to maximise the function, sorry, didn't read the command carefully enough. And I also noticed that my comments were anyhow wrong, for minimisation the points (x,x) and not (x, 1-x) would be the solutions; must have also misread the function as |x1+x2|... Looks as if I need more coffee to wake up and get my brain going... Please ignore my last two mails on this issue and apologise to the list Cheers, Berwin __ R-help@stat.math.ethz.ch mailing list 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.
Re: [R] ordered logistic regression with random effects. Howto?
On the definitional question, some texts do indeed consider multi-category
logistic regression as a glm. But the original definition by Nelder does
not. I've never seen polr considered to be a glm (but it way well have
been done).
Adding random effects is a whole different ball game: you need to
integrate over the random effects to find a likelihood. That integration
is tricky, and I am not sure we yet have reliable software for it in the
binary ('dichotomous dependent variable') case: SAS's NLMIXED certainly is
not reliable. I've had students run real problems through a variety of
software, and get quite different results. (It is possible that the shape
of the likelihood is a problem but it is not the only one.)
MCMC approaches to that integration are an alternative not mentioned
below.
On Mon, 7 May 2007, Paul Johnson wrote:
I'd like to estimate an ordinal logistic regression with a random
effect for a grouping variable. I do not find a pre-packaged
algorithm for this. I've found methods glmmML (package: glmmML) and
lmer (package: lme4) both work fine with dichotomous dependent
variables. I'd like a model similar to polr (package: MASS) or lrm
(package: Design) that allows random effects.
I was thinking there might be a trick that might allow me to use a
program written for a dichotomous dependent variable with a mixed
effect to estimate such a model. The proportional odds logistic
regression is often written as a sequence of dichotomous comparisons.
But it seems to me that, if it would work, then somebody would have
proposed it already.
You need to combine all the binary comparisons to get the likelihood, and
the models have parameters in common.
I've found some commentary about methods of fitting ordinal logistic
regression with other procedures, however, and I would like to ask for
your advice and experience with this. In this article,
Ching-Fan Sheu, Fitting mixed-effects models for repeated ordinal
outcomes with the NLMIXED procedure Behavior Research Methods,
Instruments, Computers, 2002, 34(2): 151-157.
the other gives an approach that works in SAS's NLMIXED procedure. In
this approach, one explicitly writes down the probability that each
level will be achieved (using the linear predictor and constants for
each level). I THINK I could find a way to translate this approach
into a model that can be fitted with either nlme or lmer. Has someone
done it?
What other strategies for fitting mixed ordinal models are there in R?
Finally, a definitional question. Is a multi-category logistic
regression (either ordered or not) a member of the glm family? I had
thought the answer is no, mainly because glm and other R functions for
glms never mention multi-category qualitative dependent variables and
also because the distribution does not seem to fall into the
exponential family. However, some textbooks do include the
multicategory models in the GLM treatment.
--
Brian D. Ripley, [EMAIL PROTECTED]
Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/
University of Oxford, Tel: +44 1865 272861 (self)
1 South Parks Road, +44 1865 272866 (PA)
Oxford OX1 3TG, UKFax: +44 1865 272595
__
R-help@stat.math.ethz.ch mailing list
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.
