Re: [R] Help: PLSR
Shengzhe Wu writes: I have a data set with 15 variables (first one is the response) and 1200 observations. Now I use pls package to do the plsr as below. [...] Because the trainSet has been scaled before training, I think Xtotvar should be equal to 14, but unexpectedly Xtotvar = 16562, Because the Xtotvar is the total X variation, measured by sum(X^2) (where X has been centered). With 14 variables, scaled to sd == 1, and 1200 observations, you should get Xtotvar == 14*(1200-1) == 16786. (Maybe you have 1184 observations: 14*1183 == 16562.) -- Bjørn-Helge Mevik __ 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
Re: [R] Doubt about nested aov output
Hi I think the problem is related to specifying Tratamento both as a fixed factor and in the error term. I attached a script with a reproducible example that shows a similar output. I do not know the details of the original data and the questions of interest, but maybe a model including Tratamento is more what you wanted to implement. Regards, Christoph ## R-Script library(nlme) ## Generating the data set.seed(1) ziel - rep(c(-6,8,20), each = 40) + rep(rnorm(15, 0, 20), each = 4) + rep(rnorm(60, 0, 10), each = 2) + rnorm(120, 0, 3) dat - data.frame(y = ziel, fix = factor(rep(1:3, each = 40)), R1 = factor(rep(1:15, each = 8)), R2 = factor(rep(1:60, each = 2))) ## Model including fix as fixed and random effect. res2 - aov(y ~ fix + Error(fix/R1/R2), data = dat) summary(res2) reslme2 - lme(y ~ fix , data = dat, random = ~ 1|fix/R1/R2) summary(reslme2) anova(reslme2) ## Model including fix as fixed factor. res1 - aov(y ~ fix + Error(R1/R2), data = dat) summary(res1) reslme - lme(y ~ fix , data = dat, random = ~ 1|R1/R2) summary(reslme) anova(reslme) -- Christoph Buser [EMAIL PROTECTED] Seminar fuer Statistik, LEO C13 ETH (Federal Inst. Technology) 8092 Zurich SWITZERLAND phone: x-41-44-632-4673 fax: 632-1228 http://stat.ethz.ch/~buser/ -- Spencer Graves writes: Others may know the answer to your question, but I don't. However, since I have not seen a reply, I will offer a few comments: 1. What version of R are you using? I just tried superficially similar things with the examples in ?aov in R 2.1.1 patched and consistently got F and p values. 2. My preference for this kind of thing is to use lme in library(nlme) or lmer in library(lme4). Also, I highly recommend Pinheiro and Bates (2000) Mixed-Effects Models in S and S-Plus (Springer). 3. If still want to use aov and are getting this problem in R 2.1.1, could you please provide this list with a small, self contained example that displays the symptoms that concern you? And PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html;. It might increase the speed and utility of replies. spencer graves Ronaldo Reis-Jr. wrote: Hi, I have two doubts about the nested aov output. 1) I have this: anova.ratos - aov(Glicogenio~Tratamento+Error(Tratamento/Rato/Figado)) summary(anova.ratos) Error: Tratamento Df Sum Sq Mean Sq Tratamento 2 1557.56 778.78 Error: Tratamento:Rato Df Sum Sq Mean Sq F value Pr(F) Residuals 3 797.67 265.89 Error: Tratamento:Rato:Figado Df Sum Sq Mean Sq F value Pr(F) Residuals 12 594.049.5 Error: Within Df Sum Sq Mean Sq F value Pr(F) Residuals 18 381.00 21.17 R dont make the F and P automatically, it is possible? I Like an output like this: Error: Tratamento Df Sum Sq Mean Sq F value Pr(F) Tratamento 2 1557.56 778.78 2.929 0.197 Error: Tratamento:Rato Df Sum Sq Mean Sq F value Pr(F) Residuals 3 797.67 265.89 5.372 0.0141 Error: Tratamento:Rato:Figado Df Sum Sq Mean Sq F value Pr(F) Residuals 12 594.049.5 2.339 0.0503 Error: Within Df Sum Sq Mean Sq F value Pr(F) Residuals 18 381.00 21.17 Why it not make automatic calculations? It is possible? 2) I can say that Error: Tratamento:Rato means an interaction between Tratamento and Rato? Normally the : represents an interaction, but in this output I think that it dont mean the interaction. Any explanation are welcome. Thanks Ronaldo -- Spencer Graves, PhD Senior Development Engineer PDF Solutions, Inc. 333 West San Carlos Street Suite 700 San Jose, CA 95110, USA [EMAIL PROTECTED] www.pdf.com http://www.pdf.com Tel: 408-938-4420 Fax: 408-280-7915 __ 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 !DSPAM:431b8510220677348368323! __ 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
[R] numerical intergation
how does one numerically intergate the following:
A=function(x,y)
{
xy
}
over the range: 2x0 4y10
say.
ie how would one set up the integrate function?
i forgot!__
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
Re: [R] numerical intergation
Hi Allan
you need adapt() from the adapt package.
library(adapt)
?adapt
adapt(2,lower=c(-2,4),upper=c(0,10),functn=function(z){prod(z)})
value relerr minpts lenwrk ifail
-84 7.38275e-08 165 73 0
NB: untested
HTH
rksh
On 5 Sep 2005, at 10:59, Clark Allan wrote:
how does one numerically intergate the following:
A=function(x,y)
{
xy
}
over the range: 2x04y10
say.
ie how would one set up the integrate function?
i forgot!__
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
--
Robin Hankin
Uncertainty Analyst
National Oceanography Centre, Southampton
European Way, Southampton SO14 3ZH, UK
tel 023-8059-7743
__
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
Re: [R] numerical intergation
found a solution: download the rmutil package from :
http://popgen0146uns50.unimaas.nl/~jlindsey/rcode.html
for those that are interested.
note that only works for two dimensions!
/
allan
(ranges was incorrect previously: should be:-2x0 4y10)
Clark Allan wrote:
i solved the problem:
adapt(2, lo=c(2,4), up=c(0,10))
is there any package that allows one to have infinite limits?? when
integrating over more than one variable?
Clark Allan wrote:
how does one numerically intergate the following:
A=function(x,y)
{
xy
}
over the range: 2x0 4y10
say.
ie how would one set up the integrate function?
i forgot!
__
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__
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
Re: [R] glm p-values on features
mydf - data.frame(x1=rnorm(100), x2=rnorm(100), x3=rnorm(100)) mydf$y - 0.5 * mydf$x1 + 3 * mydf$x3 + rnorm(100) fit - glm( y ~ . , data=mydf ) coefficients( summary(fit) ) Estimate Std. Errort value Pr(|t|) (Intercept) -0.2385855 0.2541498 -0.9387593 3.502103e-01 x1 0.6956811 0.1330900 5.2271462 1.003295e-06 x2 0.1313823 0.1417679 0.9267420 3.563846e-01 x3 2.7986410 0.4531338 6.1761919 1.576173e-08 On Fri, 2005-09-02 at 18:36 +0200, Joost van Evert wrote: Dear list, does anyone know how to get p-values on the coefficients returned by glm? thanks+greets, Joost __ 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 __ 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
Re: [R] Help: Plsr
Dear Bjørn-Helge, Sorry, I wrote the wrong number of observation. It should be 1184. I saw on the book that variance is defined by sd^2. If variation is a different concept from variance and defined by sd^2*(n-1) ? Since I formerly took variance and variation as the same. Thank you, Shengzhe Shengzhe Wu writes: I have a data set with 15 variables (first one is the response) and 1200 observations. Now I use pls package to do the plsr as below. [...] Because the trainSet has been scaled before training, I think Xtotvar should be equal to 14, but unexpectedly Xtotvar = 16562, Because the Xtotvar is the total X variation, measured by sum(X^2) (where X has been centered). With 14 variables, scaled to sd == 1, and 1200 observations, you should get Xtotvar == 14*(1200-1) == 16786. (Maybe you have 1184 observations: 14*1183 == 16562.) -- Bj?rn-Helge Mevik __ 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
Re: [R] Multivariate Skew Normal distribution
On Thu, 1 Sep 2005 13:09:00 + (GMT), Caio Lucidius Naberezny Azevedo wrote: CLNA CLNA Could anyone tell me if there is any package (or function) that CLNA generates values from a multivariate skew normal distribution? try the following library(sn) location - c(20,4) # e.g. for a two-dimensional variable dispers - matrix(c(3,-2,-2,2), 2,2) skew - c(10,-6) rmsn(n=10, xi=location, Omega=dispers, alpha=skew) # for skew-normal data rmst(n=10, xi=location, Omega=dispers, alpha=skew, df=5) # for skew-t data see also help(rsn) and help(rst) for the univariate case for more information, see also (as already pointed out in the list): http://azzalini.stat.unipd.it/SN best wishes, Adelchi Azzalini -- Adelchi Azzalini [EMAIL PROTECTED] Dipart.Scienze Statistiche, Università di Padova, Italia tel. +39 049 8274147, http://azzalini.stat.unipd.it/ __ 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
[R] Dummy variables model
Hi, all! Anyone know an easy way to specify the following model. Panel dataset, with stock through time, by firm. I want to run a model of y on a bunch of explanatory variables, and one dummy for each firm, which is 1 for observations that come from firm i, and 0 everywhere else. I have over 200 firms (and a factor variable that contains a firm identifier). Any easy way of going about this, without having to define all these dummies? I checked lme() with random = ~ 1|firm, but the problem is that these are random effects, i.e. that there are firm-by-firm disturbance terms and overall disturbance terms, whereas I want just overall disturbance terms. This is generally called a fixed effects model, although it seems like the term fixed effects is being used somewhat differently in the context of the nlme package. Toby -- ** When Thomas Edison invented the light bulb he tried over 2000 experiments before he got it to work. A young reporter asked him how it felt to have failed so many times. He said I never failed once. I invented the light bulb. It just happened to be a 2000-step process. __ 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
[R] model comparison and Wald-tests (e.g. in lmer)
Dear expeRts, there is obviously a general trend to use model comparisons, LRT and AIC instead of Wald-test-based significance, at least in the R community. I personally like this approach. And, when using LME's, it seems to be the preferred way (concluded from postings of Brian Ripley and Douglas Bates' article in R-News 5(2005)1), esp. because of problems with the d.f. approximation. But, on the other hand I found that not all colleagues are happy with the resulting AIC/LRT tables and the comparison of multiple models. As a compromise, and after a suggestion in Crawley's Statistical computing one may consider to supply traditional ANOVA tables as an additional explanation for the reader (e.g. field biologists). An example: one has fitted 5 models m1..m5 and after: anova(m1,m2,m3,m4,m5) # giving AIC and LRT-tests he selects m3 as the most parsimonious model and calls anova with the best model (Wald-test): anova(m3) # the additional explanatory table My questions: * Do people outside the S-PLUS/R world still understand us? * Is it wise to add such an explanatory table (in particular when the results are the same) to make the results more transparent to the reader? * Are such additional ANOVA tables *really helpful* or are they (in combination with a model comparison) just another source of confusion? Thank you! Thomas P. __ 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
[R] TeX distribution on Windows
I'm looking for a Windows distribution of TeX that works with R, after a few years' absence from Windows. On Duncan Murdoch's Rtools page fptex is still recommended, but it turns out that fptex is defunct as of May 2005, see http://www.metz.supelec.fr/~popineau/xemtex-7.html So, what is suggested? TUG (tug.org) recommends something called proTeXt, which is said to be based on MiKTeX, for Windows users. Since MikTeX could be used with R, that sounds like a good alternative. Any comments to that? -- Göran Broströmtel: +46 90 786 5223 Department of Statistics fax: +46 90 786 6614 Umeå University http://www.stat.umu.se/~goran.brostrom/ SE-90187 Umeå, Sweden e-mail: [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
[R] help
Dear helpeRs,
I seem to be a little bit confused on the result I am getting from the
few codes below:
u=v=seq(0,1,length=30)
u
[1] 0. 0.03448276 0.06896552 0.10344828 0.13793103 0.17241379
[7] 0.20689655 0.24137931 0.27586207 0.31034483 0.34482759 0.37931034
[13] 0.41379310 0.44827586 0.48275862 0.51724138 0.55172414 0.58620690
[19] 0.62068966 0.65517241 0.68965517 0.72413793 0.75862069 0.79310345
[25] 0.82758621 0.86206897 0.89655172 0.93103448 0.96551724 1.
v
[1] 0. 0.03448276 0.06896552 0.10344828 0.13793103 0.17241379
[7] 0.20689655 0.24137931 0.27586207 0.31034483 0.34482759 0.37931034
[13] 0.41379310 0.44827586 0.48275862 0.51724138 0.55172414 0.58620690
[19] 0.62068966 0.65517241 0.68965517 0.72413793 0.75862069 0.79310345
[25] 0.82758621 0.86206897 0.89655172 0.93103448 0.96551724 1.
f=function(x,y){cp=max(x+y-1,0)}
z=outer(u,v,f)
z is a 30x30 matrix which is fine, but why all its entries are equal to
1? for example, the maximum between u[22]+v[22]-1 and 0 is not 1?? I
don't really know where I went wrong!
thanks,
Dominique
__
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
Re: [R] TeX distribution on Windows
With my windows installation, MikTeX works fine, without any problem in compiling packages documentation. Antonio, Fabio Di Narzo. 2005/9/5, Göran Broström [EMAIL PROTECTED]: I'm looking for a Windows distribution of TeX that works with R, after a few years' absence from Windows. On Duncan Murdoch's Rtools page fptex is still recommended, but it turns out that fptex is defunct as of May 2005, see http://www.metz.supelec.fr/~popineau/xemtex-7.html So, what is suggested? TUG (tug.org) recommends something called proTeXt, which is said to be based on MiKTeX, for Windows users. Since MikTeX could be used with R, that sounds like a good alternative. Any comments to that? -- Göran Broströmtel: +46 90 786 5223 Department of Statistics fax: +46 90 786 6614 Umeå University http://www.stat.umu.se/~goran.brostrom/ SE-90187 Umeå, Sweden e-mail: [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 __ 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
Re: [R] TeX distribution on Windows
Hi, Göran Broström wrote: So, what is suggested? TUG (tug.org) recommends something called proTeXt, which is said to be based on MiKTeX, for Windows users. Since MikTeX could be used with R, that sounds like a good alternative. Any comments to that? I've been using MikTeX on Windows for years and have never had any problems. Its Update Wizard also has a nice and intuitive user interface. I've never had any problems using it with R. Cheers, Kev -- Ko-Kang Kevin Wang PhD Student Centre for Bioinformation Science Building 27, Room 1004 Mathematical Sciences Institute (MSI) Australian National University Canberra, ACT 2601 Australia Homepage: http://wwwmaths.anu.edu.au/~wangk/ Ph (W): +61-2-6125-2431 Ph (H): +61-2-6125-7488 Ph (M): +61-40-451-8301 __ 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
[R] [R-pkgs] New package for grouped data models
Dear R-users, We'd like to announce the release of our new package grouped (available from CRAN), for fitting models for grouped or coarse data, under the Coarsened At Random assumption. This is useful in cases where the true response variable is known only up to an interval in which it lies. Features of the package include: power calculations for two-group comparisons, computation of Bayesian residuals using Multiple Imputation, capability of fitting models for index-kind data (e.g., quality of life indexes). Any kind of feedback (questions, suggestions, bug-reports, etc.) is more than welcome. Best, Dimitris Rizopoulos Roula Tsonaka Dimitris Rizopoulos Ph.D. Student Biostatistical Centre School of Public Health Catholic University of Leuven Address: Kapucijnenvoer 35, Leuven, Belgium Tel: +32/16/336899 Fax: +32/16/337015 Web: http://www.med.kuleuven.be/biostat/ http://www.student.kuleuven.be/~m0390867/dimitris.htm Disclaimer: http://www.kuleuven.be/cwis/email_disclaimer.htm ___ 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
Re: [R] TeX distribution on Windows
[EMAIL PROTECTED] wrote: I'm looking for a Windows distribution of TeX that works with R, after a few years' absence from Windows. On Duncan Murdoch's Rtools page fptex is still recommended, but it turns out that fptex is defunct as of May 2005, see http://www.metz.supelec.fr/~popineau/xemtex-7.html So, what is suggested? TUG (tug.org) recommends something called proTeXt, which is said to be based on MiKTeX, for Windows users. Since MikTeX could be used with R, that sounds like a good alternative. Any comments to that? MikTeX works very well for us: - when writing reports using R figures, - with SWeave and - in the R package building process. You can get the web installer from: http://sourceforge.net/projects/miktex Miktex can be used with WinEDT, Emacs, Texniccenter and others as editor. HTH Thomas Petzoldt __ 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
[R] convergence for proportional odds model
Hey, everyone, I am using proportional odds model for ordinal responses in dose-response experiments. For some samll data, SAS can successfully provide estimators of the parameters, but the built-in function polr() in R fails. Would you like to tell me how to make some change so I can use polr() to obtain the estimators? Or anyone can give me a hint about the conditions for the existance of MLE in such a simple case? By the way, for the variable resp which must be ordered factor, how can I do it? Thanks a lot. Guohui The following is one example I used both in SAS and R. in R: library(MASS) dose.resp = matrix( c(1,1,1,1,2,2,2,3,3,3, 2,2,3,3,4,4,5,4,5,5), ncol=2) colnames(dose.resp)= c(resp, dose) dose.resp # dose.resp # resp dose # [1,]12 # [2,]12 # [3,]13 # [4,]13 # [5,]24 # [6,]24 # [7,]25 # [8,]34 # [9,]35 #[10,]35 polr( factor(resp, ordered=T)~dose, data=dose.resp) #Error in optim(start, fmin, gmin, method = BFGS, hessian = Hess, ...) : # initial value in 'vmmin' is not finite #In addition: Warning message: #fitted probabilities numerically 0 or 1 occurred in: #glm.fit(X, y1, wt, family = binomial(), offset = offset) in SAS NOTE: PROC LOGISTIC is fitting the cumulative logit model. The probabilities modeled are summed over the responses having the lower Ordered Values in the Response Profile table. NOTE: Convergence criterion (GCONV=1E-8) satisfied. NOTE: There were 10 observations read from the data set WORK.ONE. - Click here to donate to the Hurricane Katrina relief effort. [[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
Re: [R] TeX distribution on Windows
Thomas Petzoldt wrote: Miktex can be used with WinEDT, Emacs, Texniccenter and others as editor. Slightly off topic, if you want to get MikTeX working with Emacs and ESS, the Claus Dethlefsen has a wonderful web site (in fact, best website on this topic IMHO) http://www.math.auc.dk/~dethlef/Tips/ Cheers, Kev -- Ko-Kang Kevin Wang PhD Student Centre for Bioinformation Science Building 27, Room 1004 Mathematical Sciences Institute (MSI) Australian National University Canberra, ACT 2601 Australia Homepage: http://wwwmaths.anu.edu.au/~wangk/ Ph (W): +61-2-6125-2431 Ph (H): +61-2-6125-7488 Ph (M): +61-40-451-8301 __ 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
Re: [R] Dummy variables model
Hi If you'd like to fit a fixed effect model without random effects, you can use lm() or aov() (see ?lm and ?aov). If your variable is a factor (?factor) then you can specify your model in lm() without coding all dummy variables. Regards, Christoph Buser -- Christoph Buser [EMAIL PROTECTED] Seminar fuer Statistik, LEO C13 ETH (Federal Inst. Technology) 8092 Zurich SWITZERLAND phone: x-41-44-632-4673 fax: 632-1228 http://stat.ethz.ch/~buser/ -- Tobias Muhlhofer writes: Hi, all! Anyone know an easy way to specify the following model. Panel dataset, with stock through time, by firm. I want to run a model of y on a bunch of explanatory variables, and one dummy for each firm, which is 1 for observations that come from firm i, and 0 everywhere else. I have over 200 firms (and a factor variable that contains a firm identifier). Any easy way of going about this, without having to define all these dummies? I checked lme() with random = ~ 1|firm, but the problem is that these are random effects, i.e. that there are firm-by-firm disturbance terms and overall disturbance terms, whereas I want just overall disturbance terms. This is generally called a fixed effects model, although it seems like the term fixed effects is being used somewhat differently in the context of the nlme package. Toby -- ** When Thomas Edison invented the light bulb he tried over 2000 experiments before he got it to work. A young reporter asked him how it felt to have failed so many times. He said I never failed once. I invented the light bulb. It just happened to be a 2000-step process. __ 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 !DSPAM:431c4675196241771238468! __ 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
Re: [R] Dummy variables model
You can turn the identity vector of the firms into a factor and do lm Jean On Mon, 5 Sep 2005, Tobias Muhlhofer wrote: Hi, all! Anyone know an easy way to specify the following model. Panel dataset, with stock through time, by firm. I want to run a model of y on a bunch of explanatory variables, and one dummy for each firm, which is 1 for observations that come from firm i, and 0 everywhere else. I have over 200 firms (and a factor variable that contains a firm identifier). Any easy way of going about this, without having to define all these dummies? I checked lme() with random = ~ 1|firm, but the problem is that these are random effects, i.e. that there are firm-by-firm disturbance terms and overall disturbance terms, whereas I want just overall disturbance terms. This is generally called a fixed effects model, although it seems like the term fixed effects is being used somewhat differently in the context of the nlme package. Toby -- ** When Thomas Edison invented the light bulb he tried over 2000 experiments before he got it to work. A young reporter asked him how it felt to have failed so many times. He said I never failed once. I invented the light bulb. It just happened to be a 2000-step process. __ 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 __ 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
Re: [R] Dummy variables model
So are you guys saying to me that if I have variable firm which is the factor of all firm identifiers, I could just go lm(y ~ x + firm) and that will implicitly include a dummy for each level of factor firm, thus making this a fixed effects (aka LSDV) model? T Jean Eid wrote: You can turn the identity vector of the firms into a factor and do lm Jean On Mon, 5 Sep 2005, Tobias Muhlhofer wrote: Hi, all! Anyone know an easy way to specify the following model. Panel dataset, with stock through time, by firm. I want to run a model of y on a bunch of explanatory variables, and one dummy for each firm, which is 1 for observations that come from firm i, and 0 everywhere else. I have over 200 firms (and a factor variable that contains a firm identifier). Any easy way of going about this, without having to define all these dummies? I checked lme() with random = ~ 1|firm, but the problem is that these are random effects, i.e. that there are firm-by-firm disturbance terms and overall disturbance terms, whereas I want just overall disturbance terms. This is generally called a fixed effects model, although it seems like the term fixed effects is being used somewhat differently in the context of the nlme package. Toby -- ** When Thomas Edison invented the light bulb he tried over 2000 experiments before he got it to work. A young reporter asked him how it felt to have failed so many times. He said I never failed once. I invented the light bulb. It just happened to be a 2000-step process. __ 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 -- ** When Thomas Edison invented the light bulb he tried over 2000 experiments before he got it to work. A young reporter asked him how it felt to have failed so many times. He said I never failed once. I invented the light bulb. It just happened to be a 2000-step process. __ 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
Re: [R] model comparison and Wald-tests (e.g. in lmer)
On 9/5/05, Thomas Petzoldt [EMAIL PROTECTED] wrote: Dear expeRts, there is obviously a general trend to use model comparisons, LRT and AIC instead of Wald-test-based significance, at least in the R community. I personally like this approach. And, when using LME's, it seems to be the preferred way (concluded from postings of Brian Ripley and Douglas Bates' article in R-News 5(2005)1), esp. because of problems with the d.f. approximation. But, on the other hand I found that not all colleagues are happy with the resulting AIC/LRT tables and the comparison of multiple models. As a compromise, and after a suggestion in Crawley's Statistical computing one may consider to supply traditional ANOVA tables as an additional explanation for the reader (e.g. field biologists). An example: one has fitted 5 models m1..m5 and after: anova(m1,m2,m3,m4,m5) # giving AIC and LRT-tests he selects m3 as the most parsimonious model and calls anova with the best model (Wald-test): anova(m3) # the additional explanatory table Whether or not this is a good idea will depend on what the differences in the models are. Two mixed-effects models for the same data set can differ in their random effects specification or in the fixed-effects specification or both. The anova() function applied to a single lmer model produces a sequential anova table for the fixed-effects terms. If the models differ in the random effects specification - say the full model has random effects for slope and intercept but the restricted model has a random effect for the intercept only - then a Wald test is not appropriate (and it is not provided). In those cases I would use a likelihood ratio test and, if necessary, approximate the p-value by simulating from the null hypothesis rather than assuming a chi-squared distribution of the test statistic. Recent versions of the mlmRev package have a vignette with extensive analysis of the Exam data, including MCMC samples from the posterior distribution of the parameters. The marginal posterior distribution of the variance components are quite clearly skewed (not surprisingly, they look like scaled chi-squared distributions). Testing whether such a parameter could be zero by creating a z-statistic from the estimate and its standard error is inappropriate. Changing both the fixed-effects and the random-effects specification is tricky. I would try to do such changes in stages, settling on the fixed-effects terms first then checking the random-effects specification. My questions: * Do people outside the S-PLUS/R world still understand us? * Is it wise to add such an explanatory table (in particular when the results are the same) to make the results more transparent to the reader? * Are such additional ANOVA tables *really helpful* or are they (in combination with a model comparison) just another source of confusion? Thank you! Thomas P. __ 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 __ 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
[R] problems with outer ( was Re: help )
[[ Please use a meaningful subject line. Thank you. ]]
Dominique,
I am unable to reproduce your example in R-2.1.1. Can you reproduce this
in a fresh session of R and state which version are you using ?
u - v - seq(0, 1, length=30)
f - function(x, y){ cp=max(x+y-1,0) }
z - outer(u, v, f)
Error in outer(u, v, f) : dim- : dims [product 900] do not match the
length ofobject [1]
As shown in the R FAQ 7.17 (http://tinyurl.com/amofj), you will need to
write a wrapper function. This works for me.
wrapper - function(x, y, my.fun, ...) {
sapply(seq(along = x), FUN = function(i) f(x[i], y[i], ...))
}
z - outer(u, v, wrapper)
dim(z)
[1] 30 30
range(z)
[1] 0 1
Regards, Adai
On Mon, 2005-09-05 at 14:06 +0200, Dominique Katshunga wrote:
Dear helpeRs,
I seem to be a little bit confused on the result I am getting from the
few codes below:
u=v=seq(0,1,length=30)
u
[1] 0. 0.03448276 0.06896552 0.10344828 0.13793103 0.17241379
[7] 0.20689655 0.24137931 0.27586207 0.31034483 0.34482759 0.37931034
[13] 0.41379310 0.44827586 0.48275862 0.51724138 0.55172414 0.58620690
[19] 0.62068966 0.65517241 0.68965517 0.72413793 0.75862069 0.79310345
[25] 0.82758621 0.86206897 0.89655172 0.93103448 0.96551724 1.
v
[1] 0. 0.03448276 0.06896552 0.10344828 0.13793103 0.17241379
[7] 0.20689655 0.24137931 0.27586207 0.31034483 0.34482759 0.37931034
[13] 0.41379310 0.44827586 0.48275862 0.51724138 0.55172414 0.58620690
[19] 0.62068966 0.65517241 0.68965517 0.72413793 0.75862069 0.79310345
[25] 0.82758621 0.86206897 0.89655172 0.93103448 0.96551724 1.
f=function(x,y){cp=max(x+y-1,0)}
z=outer(u,v,f)
z is a 30x30 matrix which is fine, but why all its entries are equal to
1? for example, the maximum between u[22]+v[22]-1 and 0 is not 1?? I
don't really know where I went wrong!
thanks,
Dominique
__
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
__
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
[R] Cox model with a frailty term
Dear all, I am interested in performing survival analyses with an area-level random effect (frailty term for each area of residence of individuals). I am using coxph, and include a frailty term in the model. Here is my model: mdcox-coxph(Surv(time,censor)~ gender + age + frailty(area, dist='gauss'), data) The model provides the variance of the frailty term (between-area variance), along with a p-value. Maybe anyone knows if it is possible to obtain a standard error for the variance parameter, in order to construct a 95% confidence interval for this parameter? Second, I would like to obtain those shrunken estimates of the frailty term in each area. Do you know if there is a command allowing one to ask the frailty term for each group to be computed? I guess it should be possible to do it with predict(mdcox, type= ...), but I don't know exactly how to do. Thank you so much for your help. Best regards, Basile Chaix French National Institute of Health and Medical Research [[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
Re: [R] help
On 05-Sep-05 Dominique Katshunga wrote:
Dear helpeRs,
I seem to be a little bit confused on the result I am getting from the
few codes below:
u=v=seq(0,1,length=30)
u
[1] 0. 0.03448276 0.06896552 0.10344828 0.13793103 0.17241379
[7] 0.20689655 0.24137931 0.27586207 0.31034483 0.34482759 0.37931034
[13] 0.41379310 0.44827586 0.48275862 0.51724138 0.55172414 0.58620690
[19] 0.62068966 0.65517241 0.68965517 0.72413793 0.75862069 0.79310345
[25] 0.82758621 0.86206897 0.89655172 0.93103448 0.96551724 1.
v
[1] 0. 0.03448276 0.06896552 0.10344828 0.13793103 0.17241379
[7] 0.20689655 0.24137931 0.27586207 0.31034483 0.34482759 0.37931034
[13] 0.41379310 0.44827586 0.48275862 0.51724138 0.55172414 0.58620690
[19] 0.62068966 0.65517241 0.68965517 0.72413793 0.75862069 0.79310345
[25] 0.82758621 0.86206897 0.89655172 0.93103448 0.96551724 1.
f=function(x,y){cp=max(x+y-1,0)}
z=outer(u,v,f)
z is a 30x30 matrix which is fine, but why all its entries are equal to
1? for example, the maximum between u[22]+v[22]-1 and 0 is not 1?? I
don't really know where I went wrong!
thanks,
Dominique
A trap easy to fall into!
First, note from ?outer that FUN must be a function which
operates elementwise on arrays. In other words, given two
array arguments (x,y), for each element of x and the corresponding
element of y, it returns a value.
However, 'max' is not such a function (despite intuitive expectations).
Specifically:
max(u+v-1,0)
[1] 1
In other words, 'max' is like 'sum': it evaluates a single result
for the whole array. However, if you look at ?max and read far
enough, you will find 'pmax' -- parallel maximum -- which does
it element by element.
Look at
pmax(u+v-1,0)
So try
f-function(x,y){pmax(x+y-1,0)}
outer(u,v,f)
The reason this is a trap is that without looking at the code for
'outer' one might perhaps expect that 'outer' picks in turn each
possible pair of values (u[i],v[j]) and passes this pair as a
set of two numbers (x,y) to f(x,y). This is not so: it passes the
entire array of pairs all at once. When FUN is 'max' it returns the
maximum of this entire array for each position in the array!
Look at the code for 'outer' to see how this happens: just enter
outer
The clue is in the line
robj - array(FUN(X, Y, ...), c(dX, dY))
Best wishes,
Ted.
E-Mail: (Ted Harding) [EMAIL PROTECTED]
Fax-to-email: +44 (0)870 094 0861
Date: 05-Sep-05 Time: 16:11:22
-- XFMail --
__
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
Re: [R] TeX distribution on Windows
On 9/5/05, Göran Broström [EMAIL PROTECTED] wrote: I'm looking for a Windows distribution of TeX that works with R, after a few years' absence from Windows. On Duncan Murdoch's Rtools page fptex is still recommended, but it turns out that fptex is defunct as of May 2005, see http://www.metz.supelec.fr/~popineau/xemtex-7.html So, what is suggested? TUG (tug.org) recommends something called proTeXt, which is said to be based on MiKTeX, for Windows users. Since MikTeX could be used with R, that sounds like a good alternative. Any comments to that? -- Göran Broströmtel: +46 90 786 5223 Department of Statistics fax: +46 90 786 6614 Umeå University http://www.stat.umu.se/~goran.brostrom/ SE-90187 Umeå, Sweden e-mail: [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 As others have commented MiKTeX works with R. Duncan Murdoch's site mentions some ways to circumvent its one limitation. This is automated in the miktex-update.bat file in http://cran.r-project.org/contrib/extra/batchfiles/ for Windows XP systems. It will locate the current version of R on your system using the registry and then copy the appropriate .fd and .sty files from your R distribution to the appropriate MiKTeX directory. Assuming R and MiKTeX are installed, just issue this command without arguments: miktex-refresh at the Windows console. Using that command, all you have to do is install MiKTeX without any customizations and then run the above command to copy the mentioned files from R to MiKTeX. (If any of the TeX files in the R distribution change then you should rerun the above command after each install of R. If they have not changed you can run it or not, it does not matter.) Also in batchfiles, is the following Windows XP command: Rfind which when issued without arguments will list the location on your system of a number of programs used with R including MiKTeX. It does not actually modify any environment variables or any other aspect of your system. __ 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
Re: [R] TeX distribution on Windows
On Tue, Sep 06, 2005 at 12:07:21AM +1000, Ko-Kang Kevin Wang wrote: Slightly off topic, if you want to get MikTeX working with Emacs and ESS, the Claus Dethlefsen has a wonderful web site (in fact, best website on this topic IMHO) http://www.math.auc.dk/~dethlef/Tips/ Thanks for that tips, it looks promising (btw, it should be www.math.aau.dk). The consensus (so far) seems to be to stick to MikTeX and skip proTeXt. Thanks for all the input. Göran __ 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
Re: [R] Dummy variables model
You will need to ensure that firm is a factor and not numerical (i.e. continuous). Here is an example firm - factor( sample(1:3, 20, replace=T) ) x1 - runif(20) y- rnorm(20) summary( fit - lm( y ~ -1 + x1 + firm ) ) ... Coefficients: Estimate Std. Error t value Pr(|t|) x1-0.049640.74861 -0.0660.948 firm1 0.107320.48269 0.2220.827 firm2 0.275480.48781 0.5650.580 firm3 -0.076510.53384 -0.1430.888 NB : The -1 in the formula forces each firm to have its own intercept. Use model.matrix, you will see the dummy variables created within lm(). model.matrix( fit ) x1 firm1 firm2 firm3 1 0.6641647 0 1 0 2 0.5142712 1 0 0 3 0.2197956 1 0 0 4 0.3211675 0 1 0 5 0.1892449 1 0 0 6 0.7740754 0 0 1 7 0.3486932 0 1 0 8 0.2116816 0 0 1 9 0.2426825 0 1 0 10 0.2219768 1 0 0 11 0.9328514 1 0 0 12 0.7880405 0 0 1 13 0.8673492 0 1 0 14 0.1777998 0 1 0 15 0.3178498 1 0 0 16 0.3379726 0 0 1 17 0.9193359 1 0 0 18 0.6998152 0 1 0 19 0.2825702 0 0 1 20 0.6139586 1 0 0 Regards, Adai On Mon, 2005-09-05 at 15:53 +0100, Tobias Muhlhofer wrote: So are you guys saying to me that if I have variable firm which is the factor of all firm identifiers, I could just go lm(y ~ x + firm) and that will implicitly include a dummy for each level of factor firm, thus making this a fixed effects (aka LSDV) model? T Jean Eid wrote: You can turn the identity vector of the firms into a factor and do lm Jean On Mon, 5 Sep 2005, Tobias Muhlhofer wrote: Hi, all! Anyone know an easy way to specify the following model. Panel dataset, with stock through time, by firm. I want to run a model of y on a bunch of explanatory variables, and one dummy for each firm, which is 1 for observations that come from firm i, and 0 everywhere else. I have over 200 firms (and a factor variable that contains a firm identifier). Any easy way of going about this, without having to define all these dummies? I checked lme() with random = ~ 1|firm, but the problem is that these are random effects, i.e. that there are firm-by-firm disturbance terms and overall disturbance terms, whereas I want just overall disturbance terms. This is generally called a fixed effects model, although it seems like the term fixed effects is being used somewhat differently in the context of the nlme package. Toby -- ** When Thomas Edison invented the light bulb he tried over 2000 experiments before he got it to work. A young reporter asked him how it felt to have failed so many times. He said I never failed once. I invented the light bulb. It just happened to be a 2000-step process. __ 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 __ 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
Re: [R] Dummy variables model
here's an example data(iris) iris1-iris iris1$setosa-0 iris1[iris1$Species%in%setosa, setosa]-1 iris1$versicolor-0 iris1$virginica-0 iris1[iris1$Species%in%virginica, virginica]-1 iris1[iris1$Species%in%versicolor, versicolor]-1 iris1-iris1[, !colnames(iris1)%in%Species] summary(lm(Sepal.Length~.-1, data=iris1)) Call: lm(formula = Sepal.Length ~ . - 1, data = iris1) Residuals: Min1QMedian3Q Max -0.794236 -0.218743 0.008987 0.202546 0.731034 Coefficients: Estimate Std. Error t value Pr(|t|) Sepal.Width 0.495890.08607 5.761 4.87e-08 *** Petal.Length 0.829240.06853 12.101 2e-16 *** Petal.Width -0.315160.15120 -2.084 0.03889 * setosa2.171270.27979 7.760 1.43e-12 *** versicolor1.447700.28149 5.143 8.68e-07 *** virginica 1.147770.35356 3.246 0.00145 ** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.3068 on 144 degrees of freedom Multiple R-Squared: 0.9974, Adjusted R-squared: 0.9973 F-statistic: 9224 on 6 and 144 DF, p-value: 2.2e-16 summary(lm(Sepal.Length~.-1, data=iris)) Call: lm(formula = Sepal.Length ~ . - 1, data = iris) Residuals: Min1QMedian3Q Max -0.794236 -0.218743 0.008987 0.202546 0.731034 Coefficients: Estimate Std. Error t value Pr(|t|) Sepal.Width0.495890.08607 5.761 4.87e-08 *** Petal.Length 0.829240.06853 12.101 2e-16 *** Petal.Width -0.315160.15120 -2.084 0.03889 * Speciessetosa 2.171270.27979 7.760 1.43e-12 *** Speciesversicolor 1.447700.28149 5.143 8.68e-07 *** Speciesvirginica 1.147770.35356 3.246 0.00145 ** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.3068 on 144 degrees of freedom Multiple R-Squared: 0.9974, Adjusted R-squared: 0.9973 F-statistic: 9224 on 6 and 144 DF, p-value: 2.2e-16 On Mon, 5 Sep 2005, Tobias Muhlhofer wrote: So are you guys saying to me that if I have variable firm which is the factor of all firm identifiers, I could just go lm(y ~ x + firm) and that will implicitly include a dummy for each level of factor firm, thus making this a fixed effects (aka LSDV) model? T Jean Eid wrote: You can turn the identity vector of the firms into a factor and do lm Jean On Mon, 5 Sep 2005, Tobias Muhlhofer wrote: Hi, all! Anyone know an easy way to specify the following model. Panel dataset, with stock through time, by firm. I want to run a model of y on a bunch of explanatory variables, and one dummy for each firm, which is 1 for observations that come from firm i, and 0 everywhere else. I have over 200 firms (and a factor variable that contains a firm identifier). Any easy way of going about this, without having to define all these dummies? I checked lme() with random = ~ 1|firm, but the problem is that these are random effects, i.e. that there are firm-by-firm disturbance terms and overall disturbance terms, whereas I want just overall disturbance terms. This is generally called a fixed effects model, although it seems like the term fixed effects is being used somewhat differently in the context of the nlme package. Toby -- ** When Thomas Edison invented the light bulb he tried over 2000 experiments before he got it to work. A young reporter asked him how it felt to have failed so many times. He said I never failed once. I invented the light bulb. It just happened to be a 2000-step process. __ 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 -- ** When Thomas Edison invented the light bulb he tried over 2000 experiments before he got it to work. A young reporter asked him how it felt to have failed so many times. He said I never failed once. I invented the light bulb. It just happened to be a 2000-step process. __ 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 __ 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
Re: [R] numerical intergation
On 9/5/05, Clark Allan [EMAIL PROTECTED] wrote:
how does one numerically intergate the following:
A=function(x,y)
{
xy
}
over the range: 2x0 4y10
say.
ie how would one set up the integrate function?
i forgot!
In this particular case its separable so you could just
integrate each factor and multiply the two results
but if you want to do it in terms of A then see:
http://finzi.psych.upenn.edu/R/Rhelp02a/archive/43836.html
__
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
Re: [R] Doubt about nested aov output
That's great. In your example, the aov fit is following by a warning: res2 - aov(y ~ fix + Error(fix/R1/R2), data = dat) Warning message: Error() model is singular in: aov(y ~ fix + Error(fix/R1/R2), data = dat) Moreover, to test whether your change was statistically significant, library(nlme) supports anova with two nested models. In this case, it produced the following: anova(reslme, reslme2) Model df AIC BIClogLik Test L.Ratio p-value reslme 1 6 864.4973 881.0704 -426.2487 reslme2 2 7 866.4973 885.8325 -426.2487 1 vs 2 1.818989e-12 1 I don't know if this relates to Ronaldo Reis' question, but it certainly seems plausible. spencer graves Christoph Buser wrote: Hi I think the problem is related to specifying Tratamento both as a fixed factor and in the error term. I attached a script with a reproducible example that shows a similar output. I do not know the details of the original data and the questions of interest, but maybe a model including Tratamento is more what you wanted to implement. Regards, Christoph ## R-Script library(nlme) ## Generating the data set.seed(1) ziel - rep(c(-6,8,20), each = 40) + rep(rnorm(15, 0, 20), each = 4) + rep(rnorm(60, 0, 10), each = 2) + rnorm(120, 0, 3) dat - data.frame(y = ziel, fix = factor(rep(1:3, each = 40)), R1 = factor(rep(1:15, each = 8)), R2 = factor(rep(1:60, each = 2))) ## Model including fix as fixed and random effect. res2 - aov(y ~ fix + Error(fix/R1/R2), data = dat) summary(res2) reslme2 - lme(y ~ fix , data = dat, random = ~ 1|fix/R1/R2) summary(reslme2) anova(reslme2) ## Model including fix as fixed factor. res1 - aov(y ~ fix + Error(R1/R2), data = dat) summary(res1) reslme - lme(y ~ fix , data = dat, random = ~ 1|R1/R2) summary(reslme) anova(reslme) -- Christoph Buser [EMAIL PROTECTED] Seminar fuer Statistik, LEO C13 ETH (Federal Inst. Technology)8092 Zurich SWITZERLAND phone: x-41-44-632-4673 fax: 632-1228 http://stat.ethz.ch/~buser/ -- Spencer Graves writes: Others may know the answer to your question, but I don't. However, since I have not seen a reply, I will offer a few comments: 1. What version of R are you using? I just tried superficially similar things with the examples in ?aov in R 2.1.1 patched and consistently got F and p values. 2. My preference for this kind of thing is to use lme in library(nlme) or lmer in library(lme4). Also, I highly recommend Pinheiro and Bates (2000) Mixed-Effects Models in S and S-Plus (Springer). 3. If still want to use aov and are getting this problem in R 2.1.1, could you please provide this list with a small, self contained example that displays the symptoms that concern you? And PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html;. It might increase the speed and utility of replies. spencer graves Ronaldo Reis-Jr. wrote: Hi, I have two doubts about the nested aov output. 1) I have this: anova.ratos - aov(Glicogenio~Tratamento+Error(Tratamento/Rato/Figado)) summary(anova.ratos) Error: Tratamento Df Sum Sq Mean Sq Tratamento 2 1557.56 778.78 Error: Tratamento:Rato Df Sum Sq Mean Sq F value Pr(F) Residuals 3 797.67 265.89 Error: Tratamento:Rato:Figado Df Sum Sq Mean Sq F value Pr(F) Residuals 12 594.049.5 Error: Within Df Sum Sq Mean Sq F value Pr(F) Residuals 18 381.00 21.17 R dont make the F and P automatically, it is possible? I Like an output like this: Error: Tratamento Df Sum Sq Mean Sq F value Pr(F) Tratamento 2 1557.56 778.78 2.929 0.197 Error: Tratamento:Rato Df Sum Sq Mean Sq F value Pr(F) Residuals 3 797.67 265.89 5.372 0.0141 Error: Tratamento:Rato:Figado Df Sum Sq Mean Sq F value Pr(F) Residuals 12 594.049.5 2.339 0.0503 Error: Within Df Sum Sq Mean Sq F value Pr(F) Residuals 18 381.00 21.17 Why it not make automatic calculations? It is possible? 2) I can say that Error: Tratamento:Rato means an interaction between Tratamento and Rato? Normally the : represents an interaction, but in this output I think that it dont mean the interaction. Any explanation are welcome. Thanks Ronaldo -- Spencer Graves, PhD Senior Development Engineer PDF Solutions, Inc.
[R] Fisher's method in discriminant analysis
Hi, I'm using mda library to solve a discriminant analysis. I get results, but the thing is that I want to use Fisher's method to obtain the classification functions and I'm lost in what I should do: libraries to use, ... Can anybody give me a clue?? Thanks. Carlos Niharra López __ 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
Re: [R] help
On Mon, 5 Sep 2005, [EMAIL PROTECTED] wrote: A trap easy to fall into! So easy, in fact, that it's a FAQ. -thomas __ 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
Re: [R] Fisher's method in discriminant analysis
See this thread http://finzi.psych.upenn.edu/R/Rhelp02a/archive/34951.html Sorry, my memory has not improved since then but there are others on this list who know better this than myself. Regards, Adai On Mon, 2005-09-05 at 18:15 +0200, C NL wrote: Hi, I'm using mda library to solve a discriminant analysis. I get results, but the thing is that I want to use Fisher's method to obtain the classification functions and I'm lost in what I should do: libraries to use, ... Can anybody give me a clue?? Thanks. Carlos Niharra López __ 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 __ 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
Re: [R] TeX distribution on Windows
Goran wrote The consensus (so far) seems to be to stick to MikTeX and skip proTeXt. Thanks for all the input. Well, nothing against MikTeX, but ProteXt also works fine with R. Peter Peter L. Flom, PhD Assistant Director, Statistics and Data Analysis Core Center for Drug Use and HIV Research National Development and Research Institutes 71 W. 23rd St www.peterflom.com New York, NY 10010 (212) 845-4485 (voice) (917) 438-0894 (fax) __ 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
[R] simple line plots?
I've spent quite a lot of the day trying to construct a fairly standard line plot in R, and I have the feeling there is a simple way that I haven't discovered. I have a large vector of measurements (TIME), and each measurement falls into one of three categories (PHASE). For each PHASE value, I want a mean of the corresponding TIME measurements plotted as a point along with standard error bars. I'd like the three resulting point connected with line segments. I'd like to have two data series like this plotted on the same graph -- one in red, one in blue. Excel, as awful as it is, does this kind of graph quite easily. After sifting through the scattered documentation, the best I could do was to store the mean values of the three points, plot those three numbers against the values 1,2, and 3, then use the arrows() function to draw error bars on each one. This is a LOT of manual effort, as you can imagine (in addition to the means I have to calculate the standard errors for each point, and I still don't know how to draw each of the three line segments I need). Any suggestions? Thanks, --Ashish. - Ashish Ranpura Institute of Cognitive Neuroscience University College London 17 Queen Square LONDON WC1N 3AR __ 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
Re: [R] help
howzit dom
i just saw your mail now!
i tried the following and got an error:
u=v=seq(1,3)
u
[1] 1 2 3
v
[1] 1 2 3
f=function(x,y){max(x+y-1,0)}
z=outer(u,v,f)
Error in outer(u, v, f) : dim- : dims [product 9] do not match the
length of object [1]
it seems s if r is not performing elementwise maximisation!!!
all you need is:
f=function(x,y){pmax(x+y-1,0)}
z=outer(u,v,f)
ok lets try it!!! i used length =5 to save space
u=v=seq(0,1,length=5)
f=function(x,y){pmax(x+y-1,0)}
z=outer(u,v,f)
z
[,1] [,2] [,3] [,4] [,5]
[1,]0 0.00 0.00 0.00 0.00
[2,]0 0.00 0.00 0.00 0.25
[3,]0 0.00 0.00 0.25 0.50
[4,]0 0.00 0.25 0.50 0.75
[5,]0 0.25 0.50 0.75 1.00
think this works!
check ?pmax
see ya
allan
Dominique Katshunga wrote:
Dear helpeRs,
I seem to be a little bit confused on the result I am getting from the
few codes below:
u=v=seq(0,1,length=30)
u
[1] 0. 0.03448276 0.06896552 0.10344828 0.13793103 0.17241379
[7] 0.20689655 0.24137931 0.27586207 0.31034483 0.34482759 0.37931034
[13] 0.41379310 0.44827586 0.48275862 0.51724138 0.55172414 0.58620690
[19] 0.62068966 0.65517241 0.68965517 0.72413793 0.75862069 0.79310345
[25] 0.82758621 0.86206897 0.89655172 0.93103448 0.96551724 1.
v
[1] 0. 0.03448276 0.06896552 0.10344828 0.13793103 0.17241379
[7] 0.20689655 0.24137931 0.27586207 0.31034483 0.34482759 0.37931034
[13] 0.41379310 0.44827586 0.48275862 0.51724138 0.55172414 0.58620690
[19] 0.62068966 0.65517241 0.68965517 0.72413793 0.75862069 0.79310345
[25] 0.82758621 0.86206897 0.89655172 0.93103448 0.96551724 1.
f=function(x,y){cp=max(x+y-1,0)}
z=outer(u,v,f)
z is a 30x30 matrix which is fine, but why all its entries are equal to
1? for example, the maximum between u[22]+v[22]-1 and 0 is not 1?? I
don't really know where I went wrong!
thanks,
Dominique
__
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__
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
[R] bivirate k-function
Hi Is there anyone who can help a reporter (me) with the following: I want to use the bivirate k-function to determin if burgerbars are clustered around schools. To that purpose Ive been told to use the bivirate k-function, but my program Arcview doesnt contain that. Im stuck on a dead end how to I for a start get the data into R? Michael Holm DICAR Center for Anaytical Reporting HYPERLINK http://www.dicar.org/www.dicar.org -- No virus found in this outgoing message. Checked by AVG Anti-Virus. [[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
Re: [R] TeX distribution on Windows
Goran wrote The consensus (so far) seems to be to stick to MikTeX and skip proTeXt. Thanks for all the input. Just to clarify, I have never used proTeXt and my comment was based on having used fptex and MiKTeX and the fact that the batchfiles distribution provides specific additional R support for MiKTeX. __ 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
Re: [R] Dummy variables model
Dang! That's awesome! Being at the end of an empirical PhD in which all the econometrics was done in R, I was already a longtime R enthusiast, but you never stop learning more neat features!!! YAY to everyone involved in R's development Toby Adaikalavan Ramasamy wrote: You will need to ensure that firm is a factor and not numerical (i.e. continuous). Here is an example firm - factor( sample(1:3, 20, replace=T) ) x1 - runif(20) y- rnorm(20) summary( fit - lm( y ~ -1 + x1 + firm ) ) ... Coefficients: Estimate Std. Error t value Pr(|t|) x1-0.049640.74861 -0.0660.948 firm1 0.107320.48269 0.2220.827 firm2 0.275480.48781 0.5650.580 firm3 -0.076510.53384 -0.1430.888 NB : The -1 in the formula forces each firm to have its own intercept. Use model.matrix, you will see the dummy variables created within lm(). model.matrix( fit ) x1 firm1 firm2 firm3 1 0.6641647 0 1 0 2 0.5142712 1 0 0 3 0.2197956 1 0 0 4 0.3211675 0 1 0 5 0.1892449 1 0 0 6 0.7740754 0 0 1 7 0.3486932 0 1 0 8 0.2116816 0 0 1 9 0.2426825 0 1 0 10 0.2219768 1 0 0 11 0.9328514 1 0 0 12 0.7880405 0 0 1 13 0.8673492 0 1 0 14 0.1777998 0 1 0 15 0.3178498 1 0 0 16 0.3379726 0 0 1 17 0.9193359 1 0 0 18 0.6998152 0 1 0 19 0.2825702 0 0 1 20 0.6139586 1 0 0 Regards, Adai On Mon, 2005-09-05 at 15:53 +0100, Tobias Muhlhofer wrote: So are you guys saying to me that if I have variable firm which is the factor of all firm identifiers, I could just go lm(y ~ x + firm) and that will implicitly include a dummy for each level of factor firm, thus making this a fixed effects (aka LSDV) model? T Jean Eid wrote: You can turn the identity vector of the firms into a factor and do lm Jean On Mon, 5 Sep 2005, Tobias Muhlhofer wrote: Hi, all! Anyone know an easy way to specify the following model. Panel dataset, with stock through time, by firm. I want to run a model of y on a bunch of explanatory variables, and one dummy for each firm, which is 1 for observations that come from firm i, and 0 everywhere else. I have over 200 firms (and a factor variable that contains a firm identifier). Any easy way of going about this, without having to define all these dummies? I checked lme() with random = ~ 1|firm, but the problem is that these are random effects, i.e. that there are firm-by-firm disturbance terms and overall disturbance terms, whereas I want just overall disturbance terms. This is generally called a fixed effects model, although it seems like the term fixed effects is being used somewhat differently in the context of the nlme package. Toby -- ** When Thomas Edison invented the light bulb he tried over 2000 experiments before he got it to work. A young reporter asked him how it felt to have failed so many times. He said I never failed once. I invented the light bulb. It just happened to be a 2000-step process. __ 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 -- ** When Thomas Edison invented the light bulb he tried over 2000 experiments before he got it to work. A young reporter asked him how it felt to have failed so many times. He said I never failed once. I invented the light bulb. It just happened to be a 2000-step process. __ 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
Re: [R] TeX distribution on Windows
On 9/5/05, Gabor Grothendieck [EMAIL PROTECTED] wrote: Goran wrote The consensus (so far) seems to be to stick to MikTeX and skip proTeXt. Thanks for all the input. Just to clarify, I have never used proTeXt and my comment was based on having used fptex and MiKTeX and the fact that the batchfiles distribution provides specific additional R support for MiKTeX. One other comment. If anyone is using proTeXt or other TeX distribution and wants to modify miktex-update.bat and Rfind.bat to support proTeXt or other TeX distribution I would be happy to provide a link to it from the batchfiles documentation or even include it in batchfiles, if desired. __ 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
Re: [R] How to set starting values for lmer?
Dear Professor Bates,
Thanks, that will probably do the job.
By the way, how to cite lme4 in my work?
Shige
On 8/31/05, Douglas Bates [EMAIL PROTECTED] wrote:
On 8/30/05, Shige Song [EMAIL PROTECTED] wrote:
Dear All,
Can anyone give me some hints about how to set starting values for a
lmer
model? For complicated models like this, good starting values can help
the
numerical computation and make the model converge faster. Thanks!
Shige
I agree but I haven't gotten around to designing how that could be
done. It could be easy or difficult depending on how you want to
represent the starting values.
If you look at the (only) lmer method function you will see that it
has a section
if (lmm) { ## linear mixed model
.Call(lmer_initial, mer, PACKAGE=Matrix)
.Call(lmer_ECMEsteps, mer, cv$niterEM, cv$EMverbose,
PACKAGE = Matrix)
LMEoptimize(mer) - cv
for linear mixed models. The object mer is a mixed-effects
representation and the list cv is the control values. The only
thing that the C function lmer_initial does is set the initial
values of the relative precision matrices for the random effects.
These are the inverses of the variance-covariance matrices relative to
the variance of the per-observation noise term. They are stored
(upper triangle only) in a slot called Omega of the mer class (which
is contained in the lmer class).
There is no purpose in setting initial values for the fixed-effects
parameters or the variance of the per-observation noise term because
these are profiled out of the optimization. The optimization is only
over the values in the Omega slot.
I can allow those values to be set from an argument and only call
lmer_initialize if that argument is missing. Will that be
sufficient for you?
[[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
Re: [R] TeX distribution on Windows
On Mon, Sep 05, 2005 at 01:55:39PM -0400, Peter Flom wrote: Goran wrote The consensus (so far) seems to be to stick to MikTeX and skip proTeXt. Thanks for all the input. Well, nothing against MikTeX, but ProteXt also works fine with R. I have made some studying; it seems as if proTeXt is nothing but MikTeX together with WinEdt or TeXnicsCenter and ghostscript/ghostview, so the choice is really between WinEdt and (X)Emacs; you get MikTeX in both cases. Göran __ 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
Re: [R] TeX distribution on Windows
On 9/5/05, Göran Broström [EMAIL PROTECTED] wrote: On Mon, Sep 05, 2005 at 01:55:39PM -0400, Peter Flom wrote: Goran wrote The consensus (so far) seems to be to stick to MikTeX and skip proTeXt. Thanks for all the input. Well, nothing against MikTeX, but ProteXt also works fine with R. I have made some studying; it seems as if proTeXt is nothing but MikTeX together with WinEdt or TeXnicsCenter and ghostscript/ghostview, so the choice is really between WinEdt and (X)Emacs; you get MikTeX in both cases. In that case it probably uses the same registry entries in which case the batch utilities I mentioned would likely work with it without any change. __ 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
Re: [R] bivirate k-function
On Mon, 5 Sep 2005, Michael Holm wrote: Hi Is there anyone who can help a reporter (me) with the following: I want to use the bivirate k-function to determin if burgerbars are clustered around schools. To that purpose Ive been told to use the bivirate k-function, but my program Arcview doesnt contain that. Im stuck on a dead end how to I for a start get the data into R? The file format prefered by ArcView for point data is shapefile - see packages maptools or shapefiles on CRAN. Bivariate K-hat functions are available in splancs (k12hat), or in spatstat (Kcross and Kdot). A direct reference to the splancs functions is given by Bailey and Gatrell (1995) Interactive Spatial Data Analysis, Ch. 4. Spatstat is described clearly in: http://www.jstatsoft.org/counter.php?id=113url=v12/i06/v12i06.pdfct=1 You may need to import a bounding polygon too. Please feel free to follow this up on the R-sig-geo list, access easiest from the CRAN Spatial task view: http://cran.r-project.org/src/contrib/Views/Spatial.html (Note though, that schools may be close to other confounding locations attracting burger-vores, so you may need a school just by itself to draw any sensible conclusions) Michael Holm DICAR Center for Anaytical Reporting HYPERLINK http://www.dicar.org/www.dicar.org -- Roger Bivand Economic Geography Section, Department of Economics, Norwegian School of Economics and Business Administration, Helleveien 30, N-5045 Bergen, Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43 e-mail: [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
Re: [R] How to set starting values for lmer?
On 9/5/05, Shige Song [EMAIL PROTECTED] wrote:
Dear Professor Bates,
Thanks, that will probably do the job.
OK, I will add that capability.
By the way, how to cite lme4 in my work?
For now I suggest citing either the package itself or the R News
article that I wrote about it.
@Article{Rnews:Bates:2005,
author = {Douglas Bates},
title= {Fitting Linear Mixed Models in {R}},
journal = {R News},
year = 2005,
volume = 5,
number = 1,
pages= {27--30},
month= {May},
url = {http://CRAN.R-project.org/doc/Rnews/},
}
Eventually there will be a book that you can cite but I have to finish
writing it first :-)
Shige
On 8/31/05, Douglas Bates [EMAIL PROTECTED] wrote:
On 8/30/05, Shige Song [EMAIL PROTECTED] wrote:
Dear All,
Can anyone give me some hints about how to set starting values for a
lmer
model? For complicated models like this, good starting values can help
the
numerical computation and make the model converge faster. Thanks!
Shige
I agree but I haven't gotten around to designing how that could be
done. It could be easy or difficult depending on how you want to
represent the starting values.
If you look at the (only) lmer method function you will see that it
has a section
if (lmm) { ## linear mixed model
.Call(lmer_initial, mer, PACKAGE=Matrix)
.Call(lmer_ECMEsteps, mer, cv$niterEM, cv$EMverbose,
PACKAGE = Matrix)
LMEoptimize(mer) - cv
for linear mixed models. The object mer is a mixed-effects
representation and the list cv is the control values. The only
thing that the C function lmer_initial does is set the initial
values of the relative precision matrices for the random effects.
These are the inverses of the variance-covariance matrices relative to
the variance of the per-observation noise term. They are stored
(upper triangle only) in a slot called Omega of the mer class (which
is contained in the lmer class).
There is no purpose in setting initial values for the fixed-effects
parameters or the variance of the per-observation noise term because
these are profiled out of the optimization. The optimization is only
over the values in the Omega slot.
I can allow those values to be set from an argument and only call
lmer_initialize if that argument is missing. Will that be
sufficient for you?
[[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
__
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
Re: [R] partial association model
I'm not familiar with the term partial association model, and I don't have the Sloane and Morgan paper you cite. Neither Google nor RSiteSearch(partial association model) produced anything that looked to me like what you were asking. Part of the R culture is that one can fit anything with R; the only question is how. To determine that and to answer your other question on the advantages of log-linear models, I believe you could best approach those questions by consulting the references in the help file for glm. In particular, I recommend you start with Venables and Ripley (2002) Modern Applied Statistics with S, 4th ed. (Springer). It is probably the best known book on this issue. If you are not already familiar with this book, I highly recommend it. I don't have a copy at my fingertips now, but with luck, you may find answers to both your questions. It may not discuss partial association models by that name, but if you know partial association models, you might be able to find something relevant in that book. McCullagh P. and Nelder, J. A. (1989) Generalized Linear Models. (London: Chapman and Hall) is the second and substantially enlarged edition of the initial defining book on this topic. Hastie, T. J. and Pregibon, D. (1992) Generalized linear models, Chapter 6 of _Statistical Models in S_ eds J. M. Chambers and T. J. Hastie, (Wadsworth Brooks/Cole) may provide useful information not available in Venables and Ripley. I doubt if this answered your question, but I hope it at least will help you in some way. Feel free to post another question. However, I think you will help yourself if you read and follow the posting guide! http://www.R-project.org/posting-guide.html;. I suspect that on average, posts prepared following this guide get more useful answers quicker. spencer graves [EMAIL PROTECTED] wrote: my last post was filtered,so I post it again with another title. If I do not make a mistake,the partial association model is an extension of log-linear model.I read a papers which gives an example of it.(Sloane and Morgan,1996,An Introduction to Categorical Data Analysis,Annual Review of Sociology.22:351-375) Can R fit such partial association model? ps:Another somewhat off-topic question.What's the motivations to use log-linear model?Or why use log-linear model?I learn the log-linear model butI still do not master the the advantage of the model.thank you! __ 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 -- Spencer Graves, PhD Senior Development Engineer PDF Solutions, Inc. 333 West San Carlos Street Suite 700 San Jose, CA 95110, USA [EMAIL PROTECTED] www.pdf.com http://www.pdf.com Tel: 408-938-4420 Fax: 408-280-7915 __ 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
[R] tcltk, X11 protocol error: Bug?
Hi,
I am having trouble debugging this one. The code is attached below, but
it seems to be a problem at the
C-tk interface. If I run this 1 time there are no problems if I run it
more than once I start to get warnings
that increase in multiples of 11 everytime I run it. Here is a sample
session
source(clrramp2.r)
Loading required package: tcltk
Loading Tcl/Tk interface ... done
clrRamp()
tt-clrRamp()
tt
function (n)
{
x - ramp(seq(0, 1, length = n))
rgb(x[, 1], x[, 2], x[, 3], max = 255)
}
environment: 0x8b8674c
image(matrix(1:10),col=tt(10))
tt-clrRamp()
There were 22 warnings (use warnings() to see them)
image(matrix(1:10),col=tt(10))
There were 11 warnings (use warnings() to see them)
warnings()
Warning messages:
1: X11 protocol error: BadWindow (invalid Window parameter)
2: X11 protocol error: BadWindow (invalid Window parameter)
3: X11 protocol error: BadWindow (invalid Window parameter)
4: X11 protocol error: BadWindow (invalid Window parameter)
5: X11 protocol error: BadWindow (invalid Window parameter)
6: X11 protocol error: BadWindow (invalid Window parameter)
7: X11 protocol error: BadWindow (invalid Window parameter)
8: X11 protocol error: BadWindow (invalid Window parameter)
9: X11 protocol error: BadWindow (invalid Window parameter)
10: X11 protocol error: BadWindow (invalid Window parameter)
11: X11 protocol error: BadWindow (invalid Window parameter)
I am running R-2.1.1 on ubuntu linux 5.04, compiled from source (not the
deb package).
My version of tcl/tk is 8.4. The code is below. If anyone sees something
I am doing foolish
let me know, otherwise I will file a bug report.
Nicholas
# File clrramp2.r ##
require(tcltk)
clrRamp - function(n.col, b.color=NULL,e.color=NULL){
B.ChangeColor - function()
{
b.color - tclvalue(tkcmd(tk_chooseColor,initialcolor=e.color,
title=Choose a color))
if (nchar(b.color)0){
tkconfigure(canvas.b,bg=b.color)
Rmp.Draw()
}
}
E.ChangeColor - function()
{
e.color - tclvalue(tkcmd(tk_chooseColor,initialcolor=e.color,
title=Choose a color))
##cat(e.color)
if (nchar(e.color)0){
tkconfigure(canvas.e,bg=e.color)
Rmp.Draw()
}
}
Rmp.Draw -function(){
cr-colorRampPalette(c(b.color,e.color),space=Lab,interpolate=spline)
rmpcol - cr(n.col)
#rmpcol-rgb( rmpcol[,1],rmpcol[,2],rmpcol[,3])
inc - 300/n.col
xl - 0
for(i in 1:n.col){
##item -
tkitemconfigure(canvas.r,barlst[[i]],fill=rmpcol[i],outline=rmpcol[i])
#xl - xl+inc
}
}
save.ramp - function(){
cr-colorRampPalette(c(b.color,e.color),space=Lab,interpolate=spline)
tkdestroy(tt)
##invisible(cr)
}
tt - tktoplevel()
tkwm.title(tt,Color Ramp Tool)
frame - tkframe(tt)
bframe - tkframe(frame,relief=groove,borderwidth=1)
if(is.null(b.color)) b.color - blue
if(is.null(e.color)) e.color - yellow
if(missing(n.col)) n.col - 100
canvas.b - tkcanvas(bframe,width=50,height=25,bg=b.color)
canvas.e - tkcanvas(bframe,width=50,height=25,bg=e.color)
canvas.r - tkcanvas(tt,width=300,height=50,bg=white)
BColor.button - tkbutton(bframe,text=Begin
Color,command=B.ChangeColor)
##tkgrid(canvas.b,BColor.button)
EColor.button - tkbutton(bframe,text=End
Color,command=E.ChangeColor)
killbutton - tkbutton(bframe,text=Save,command=save.ramp)
tkgrid(canvas.b,BColor.button,canvas.e,EColor.button)
tkgrid(bframe)
tkgrid(frame)
tkgrid(canvas.r)
tkgrid(killbutton)
cr-colorRampPalette(c(b.color,e.color),space=Lab,interpolate=spline)
##rmpcol - hex(mixcolor(alpha,bc,ec,where=LUV))
rmpcol - cr(n.col)
inc - 300/n.col
xl - 0
#barlst - vector(length=n.col,mode=list)
barlst - tclArray()
for(i in 1:n.col){
item-tkcreate(canvas.r,rect,xl,0,xl+inc,50,
fill=rmpcol[i],outline=rmpcol[i])
##tkaddtag(canvas.r, point, withtag, item)
barlst[[i]]-item
xl - xl+inc
}
tkgrab.set(tt)
tkwait.window(tt)
##tkdestroy(tt)
invisible(cr)
}
__
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
Re: [R] question sur R
Thanks for providing such a concise example. I ran your example and got the same error message with R 2.1.1 patched and with the MASS package version 7.2-19. Moreover, I got the same error with the following simplification: test1 = polr(rating ~ X1) summary(test1) There are two possible next steps from this point: (1) Copy the code for polr into a script file and work through it line by line until you understand where the error message arises and why. (2) Report the problem to the package maintainer, who in this case is Prof. Brian Ripley [EMAIL PROTECTED]. In reporting your problem, please preceed your script with something like set.seed(1). Because you use pseudo-random numbers, someone else may get a different result unless you both start with the same seed. Also, please report which version of R and the MASS package you used. And you can certainly add that I get essentially the same result. If it were my problem, I would try these two steps in this order. However, a polite note to Prof. Ripley asking about this without tracing the error yourself would also be acceptable. I'm sorry I was not of more help, but I don't have time to trace this problem myself. spencer graves Abdelhafid BERRICHI wrote: hello I've tried to simulate a normal law, like that : X1 = c(rnorm(90,50,5583),rnorm(160,1198,13034597),rnorm(40,13,125)) then, I've regressed my ordinal polytomic variable rating rating=c(rep(2,90),rep(3,160), rep(4,40)) rating = as.factor(rating) rating = as.ordered(rating) (ratins is an ordered factor) on the continuous variable X1 like that library(MASS) test2 = polr(rating ~ X1, method = c(probit)) summary(test2) but R indicates me the following error whereas X does not have infinite or missing values? Re-fitting to get Hessian Error in svd(X) : infinite or missing values in x THANKS good by abdelhafid berrichi [[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 -- Spencer Graves, PhD Senior Development Engineer PDF Solutions, Inc. 333 West San Carlos Street Suite 700 San Jose, CA 95110, USA [EMAIL PROTECTED] www.pdf.com http://www.pdf.com Tel: 408-938-4420 Fax: 408-280-7915 __ 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
Re: [R] TeX distribution on Windows
Göran Broström wrote: I'm looking for a Windows distribution of TeX that works with R, after a few years' absence from Windows. On Duncan Murdoch's Rtools page fptex is still recommended, but it turns out that fptex is defunct as of May 2005, see http://www.metz.supelec.fr/~popineau/xemtex-7.html So, what is suggested? TUG (tug.org) recommends something called proTeXt, which is said to be based on MiKTeX, for Windows users. Since MikTeX could be used with R, that sounds like a good alternative. I use MikTeX, with one or another of the workarounds listed on my page. I've never tried proTeXt; I did a little googling, but I still don't see the point of it exactly. fptex is still available in various repositories, and is likely to keep working for quite a long time: R doesn't demand the latest and greatest innovations from TeX/eTeX. Duncan Murdoch __ 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
[R] MASS: rlm, MM and errors in observations AND regressors
Hello, I need to perform a robust regression on data which contains errors in BOTH observations and regressors. Right now I am using rlm from the MASS package with 'method=MM' and get visually very nice results. MASS is quite clear, however, that the described methodologies are only applicable to observation-error only data (p. 157, 4th Ed.). So here's the questions now: a) is there methodology for robust and resistant regression on data with errors in BOTH observations and regressors? b) if yes at a): does somebody know of an implementation in R? c) if no at b): would somebody reading this be open to implement this? If yes: please get in contact with me off list. Please excuse my unfamiliarity with terminology and subject - I'm very new to statistics. Joh -- +--+ | Johannes Graumann, Dipl. Biol. | | | | Graduate StudentTel.: ++1 (626) 395 6602| | Deshaies LabFax.: ++1 (626) 395 5739| | Department of Biology | | CALTECH, M/C 156-29 | | 1200 E. California Blvd.| | Pasadena, CA 91125 | | USA | +--+ __ 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
Re: [R] How to set starting values for lmer?
Dear Professor Bates,
Thank you very much for the help. I cannot waiting to see your new book!
Best,
Shige
On 9/6/05, Douglas Bates [EMAIL PROTECTED] wrote:
On 9/5/05, Shige Song [EMAIL PROTECTED] wrote:
Dear Professor Bates,
Thanks, that will probably do the job.
OK, I will add that capability.
By the way, how to cite lme4 in my work?
For now I suggest citing either the package itself or the R News
article that I wrote about it.
@Article{Rnews:Bates:2005,
author = {Douglas Bates},
title = {Fitting Linear Mixed Models in {R}},
journal = {R News},
year = 2005,
volume = 5,
number = 1,
pages = {27--30},
month = {May},
url = {http://CRAN.R-project.org/doc/Rnews/},
}
Eventually there will be a book that you can cite but I have to finish
writing it first :-)
Shige
On 8/31/05, Douglas Bates [EMAIL PROTECTED] wrote:
On 8/30/05, Shige Song [EMAIL PROTECTED] wrote:
Dear All,
Can anyone give me some hints about how to set starting values for a
lmer
model? For complicated models like this, good starting values can
help
the
numerical computation and make the model converge faster. Thanks!
Shige
I agree but I haven't gotten around to designing how that could be
done. It could be easy or difficult depending on how you want to
represent the starting values.
If you look at the (only) lmer method function you will see that it
has a section
if (lmm) { ## linear mixed model
.Call(lmer_initial, mer, PACKAGE=Matrix)
.Call(lmer_ECMEsteps, mer, cv$niterEM, cv$EMverbose,
PACKAGE = Matrix)
LMEoptimize(mer) - cv
for linear mixed models. The object mer is a mixed-effects
representation and the list cv is the control values. The only
thing that the C function lmer_initial does is set the initial
values of the relative precision matrices for the random effects.
These are the inverses of the variance-covariance matrices relative to
the variance of the per-observation noise term. They are stored
(upper triangle only) in a slot called Omega of the mer class (which
is contained in the lmer class).
There is no purpose in setting initial values for the fixed-effects
parameters or the variance of the per-observation noise term because
these are profiled out of the optimization. The optimization is only
over the values in the Omega slot.
I can allow those values to be set from an argument and only call
lmer_initialize if that argument is missing. Will that be
sufficient for you?
[[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
[[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
