Re: [R] plot p(Y=1) vs as
I'd recommend having a look at the lrm function in the Design package. Looking at the examples should show you how you can product this type of plot using some of the other components of this package. On 26/11/06, Aimin Yan [EMAIL PROTECTED] wrote: I am trying to fit a logistic regression model for this data set. Firstly, I want to plot P(Y=1) vs As and P(Y=1) vs Aa. Does any body know how to do these in R. Thanks, Aimin p5 - read.csv(http://www.public.iastate.edu/~aiminy/data/p_5_2.csv;) str(p5) 'data.frame': 1030 obs. of 6 variables: $ P : Factor w/ 5 levels 821p,8ABP,..: 1 1 1 1 1 1 1 1 1 1 ... $ Aa : Factor w/ 19 levels ALA,ARG,ASN,..: 12 16 7 18 11 10 19 19 19 1 ... $ As : num 126.9 64.1 82.7 7.6 42.0 ... $ Ms : num 92.4 50.7 75.3 17.2 57.7 ... $ Cur: num -0.1320 -0.0977 -0.0182 0.2368 0.1306 ... $ Y : int 0 0 1 1 0 0 1 0 1 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. -- = David Barron Said Business School University of Oxford Park End Street Oxford OX1 1HP __ 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] Plot with two seperate y axis
Thorsten Muehge wrote: Hello, I would like to plot the following matrix: Wk=x achsis. Para 1 = left y-axis as a barplot para 2 right y-axis as a normal scatter plat. I could not find such a solution in any of my documentation. Can someone help me? Thanks a lot Thorsten WkPara 1 Para 2 312000 99.8 322005 99.0 332002 98.0 341090 98.5 352001 99.1 362010 97.0 372010 98.8 Hi Thorsten, Is this what you mean? par(mar=c(5,4,4,4)) xpos-barplot(mueghe.df$Para.1) axis(1,at=xpos,labels=mueghe.df$Wk) par(new=TRUE) plot(xpos,mueghe.df$Para.2,xlim=par(usr)[1:2],xaxs=i, ylim=c(0,100),type=n,xlab=,ylab=,axes=FALSE) axis(4) points(xpos,mueghe.df$Para.2) Jim __ 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] cat not evaluated before file.choose
I believe this is specific to the MacOS (and Windows) GUI, which you did
not mention you were using.
See ?flush.console for the solution. Your subject line is misleading: it
is the console that is delaying showing the result.
On Sat, 25 Nov 2006, Kevin Middleton wrote:
As part of a larger function I have code similar to the reduced
example below. The user is instructed to choose a file, which gets
read using read.csv. In this example, I just have the name of the
file print out.
When I call this function with choosefile(), the file dialog box
appears before the first cat line is printed. After I choose a file,
both cat lines are printed.
choosefile - function (){
cat(Choose the data file.\n)
filename - file.choose(new = FALSE)
cat(You chose: , filename, sep = )
}
Is there a way to force the first cat line to print before the call
to file.choose? I'm using R 2.4.0 Patched (2006-11-24 r39989) on OS
X. Session info below.
Any suggestions would be greatly appreciated.
Kevin Middleton
---
sessionInfo()
R version 2.4.0 Patched (2006-11-24 r39989)
powerpc-apple-darwin8.8.0
locale:
en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
attached base packages:
[1] stats graphics grDevices utils datasets
methods base
--
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] Nonlinear statistical modeling -- a comparison of R and AD Model Builder
Spencer, I tried the mixed effects approach you suggest using the random effects module of AD Model Builder: (http://www.otter-rsch.ca/admbre/admbre.html). What are 94 unbounded parameters in Schnute et al (1998), now become realizations of a Gaussian random variable, with the corresponding standard deviation being estimated as a parameter. The approach works, but the computation time is increased substantially. This is however understandable as the computational problem is a very different one. The likelihood function now involves an integral in dimension 94, which I believe cannot be broken into a product of lower dimensional integrals as is usual for clustered data (the reason being the recursive nature of the population dynamics). hans ___ Spencer Graves wrote: Have you considered nonlinear mixed effects models for the types of problems considered in the comparison paper you cite? Those benchmark trials consider T years of data ... for A age classes and the total number of parameters is m = T+A+5. Without knowing more about the problem, I suspect that the T year parameters and the A age class parameters might be better modeled as random effects. If this were done, the optimization problem would then involve 7 parameters, the 5 fixed-effect parameters suggested by the computation of m plus two variance parameters, one for the random year effects and another for the random age class effect. This would replace the problem of maximizing, e.g., a likelihood over T+A+5 parameters with one of maximizing a marginal likelihood over 2+5 parameters after integrating out the T and A random effects. These integrations may not be easy, and I might stick with the fixed-effects solution if I couldn't get answers in the available time using a model I thought would be theoretically more appropriate. Also, I might use the fixed-effects solution to get starting values for an attempt to maximize a more appropriate marginal likelihood. For the latter, I might first try 'nlmle'. If that failed, I might explore Markov Chain Monte Carlo (MCMC). I have not done MCMC myself, but the MCMCpack R package looks like it might make it feasible for the types of problems considered in this comparison. The CRAN summary of that package led me to an Adobe Acrobat version of a PPT slide presentation that seemed to consider just this type of problem (e.g., http://mcmcpack.wustl.edu/files/MartinQuinnMCMCpackslides.pdf). Have you considered that? Hope this helps. Spencer Graves [[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] axis in an image() plot
Ricardo Rodríguez - Your XEN ICT Team wrote: Hi all, I've found quite usefull colored-grid created by image() but I'm facing a doubt I am not able to solve. Given the following data rectangle... 1 2 3 4 5 6 7 8 9 10 11 12 13 14 1 12 22 0 7 2 1 0 2 0 2 6 -3 0 3 2 0 -1 0 9 3 -4 0 0 0 0 3 0 0 0 3 29 45 6 12 16 85 -2 0 -3 -4 89 -1 -1 1 4 2 9 3 6 17 3 -2 -9 -2 8 -1 0 0 0 5 44 16 -3 21 23 3 2 1 0 -2 13 18 -5 2 I am not able to draw x and y axis labeled 1 to 14 and 1 to 5 by 1. I've tried a number of options by using axis() to no avail. It will be also very helpfull to be able to draw the value of each x,y combination within its position in the grid. Text() seems to be the answer, but I am not able yet to get the correct position for each label. Please, could you point me in the right direction or offer some example? Hi Ricardo, This might be what you want (say your data frame is called my.df): library(plotrix) color2D.matplot(my.df,c(1,0),c(0,0),c(0,1)) text(rep(0.5:13.5,5),rep(seq(4.5,0.5,by=-1),14), unlist(my.df),col=white) and in fact it looks so neat that I might add it as an option. Jim __ 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] Questions about generating samples in R
Hello! I have a data set with 8 columns and in about 5000 rows. What I want to do is to generate samples of this data set. Samples of a special size, as example 200. What is the easiest way to do this? No special things are needed, only the random selection of 200 rows of the data set. Thanks Alex -- Alexander Geisler * Kaltenbach 151 * A-6272 Kaltenbach email: [EMAIL PROTECTED] | [EMAIL PROTECTED] phone: +43 650 / 811 61 90 | skpye: al1405ex __ 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] Questions about generating samples in R
?sample should tell you what you need to know. On 26/11/06, Alexander Geisler [EMAIL PROTECTED] wrote: Hello! I have a data set with 8 columns and in about 5000 rows. What I want to do is to generate samples of this data set. Samples of a special size, as example 200. What is the easiest way to do this? No special things are needed, only the random selection of 200 rows of the data set. Thanks Alex -- Alexander Geisler * Kaltenbach 151 * A-6272 Kaltenbach email: [EMAIL PROTECTED] | [EMAIL PROTECTED] phone: +43 650 / 811 61 90 | skpye: al1405ex __ 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. -- = David Barron Said Business School University of Oxford Park End Street Oxford OX1 1HP __ 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] Fitting mixed-effects models with lme with fixed error term variances
Hello! Douglas Bates wrote: ... Not quite. There is now a capability in lmer to specify starting estimates for the relative precision matrices which means that you can use the estimates from one fit as starting estimates for another fit. It looks like fm1 - lmer(...) fm2 - lmer(y ~ x + (...), start = [EMAIL PROTECTED]) I should say that in my experience this has not been as effective as I had hoped it would be. What I see in the optimization is that the first few iterations reduce the deviance quite quickly but the majority of the time is spent refining the estimates near the optimum. So, for example, if it took 40 iterations to converge from the rather crude starting estimates calculated within lmer it might instead take 35 iterations if you give it good starting estimates. Yes, I also have the same experience with other programs, say VCE[1]. What I was hopping for was just solutions from linear mixed model i.e. either via GLS or PLS representations and no optimization if values for variance components are passed as input - there should be a way to distinguish that from just passing starting values.. The PLS representation in lmer is in terms of the relative variance-covariance matrices of the random effects where relative means relative to s^2. (In fact, the Omega slot contains the relative What exactly is s^2? precision matrices (inverses of the relative variance matrices) but its the same idea. All variance components are expressed relative to s^2.) If it is sufficient to fix these then you can easily do that. Just follow the code in lmer until it gets to .Call(mer_ECMEsteps, mer, cv$niterEM, cv$EMverbose) LMEoptimize(mer) - cv return(new(lmer, mer, and replace the first two lines by .Call(mer_factor, mer, PACKAGE = lme4) Hmm. Sounds simple. I will try this. Do you find any value in adding this as an option in lmer or perhaps new? function, say lmeSol() or something similar? A side-effect of performing the factorization is to evaluate the ML and REML deviances and store the result in the deviance slot. The follow the code in lmer part will get a bit tricky because of the number of functions that are hidden in the lme4 namespace but I'm sure you know how to get around that. Thank you very much! Gregor __ 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] axis in an image () plot
jim holtman[EMAIL PROTECTED] 24/11/2006 23:21 If you data is a matrix, then try: image(1:5, 1:14, data.rect) text(row(data.rect), col(data.rect), data.rect) Thanks, Jim. It works great and perfectly fill my requirements. I better understand now how text() works. By the way, there is a line at the bottom of your message reading What is the problem you are trying to solve?, is this a kind of motto or are you asking me what I am trying to solve? :-) Thanks! Cheers, Ricardo -- Ricardo Rodríguez Your XEN ICT Team __ 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] axis in an image () plot
Jim Lemon[EMAIL PROTECTED] 26/11/2006 12:05 Hi Ricardo, This might be what you want (say your data frame is called my.df): library(plotrix) color2D.matplot(my.df,c(1,0),c(0,0),c(0,1)) text(rep(0.5:13.5,5),rep(seq(4.5,0.5,by=-1),14), unlist(my.df),col=white) and in fact it looks so neat that I might add it as an option. Jim Thanks, Jim! Once the original problem was solved by using image(), your plotrix() package is of major interest to keep improving this kind of graphics! Thanks for your contribution, Ricardo -- Ricardo Rodríguez Your XEN ICT Team __ 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] nonlinear regression-getting the explained variation
On 11/23/06, [EMAIL PROTECTED] [EMAIL PROTECTED] wrote: I'm trying to teach myself R, and by the way, re-learning statistics using Crawley's Statistics: an introduction using R. I've reached the regression chapter, and when it deals with non-linear regresion using the nls library I face the following problem: I follow the steps--- deer-read.table(c:\\temp\\jaws.txt,header=T) ---data available at http://www.bio.ic.ac.uk/research/crawley/statistics/data/ attach(deer) names(deer) library(nls) model2-nls(bone~a*(1-exp(-c*age)),start=list(a=120,c=0.064)) summary(model2) I've taken away those steps that regard plotting. But I hope you did actually plot the data when you tried this exercise. Also, you can, if you wish, read the data from the URL without copying the data to a local file. deer - read.table(http://www.bio.ic.ac.uk/research/crawley/statistics/data/jaws.txt;, header = TRUE) model2-nls(bone~a*(1-exp(-c*age)),deer, start=list(a=120,c=0.064)) summary(model2) Formula: bone ~ a * (1 - exp(-c * age)) Parameters: Estimate Std. Error t value Pr(|t|) a 115.580562.84365 40.645 2e-16 c 0.118820.01233 9.635 3.69e-13 Residual standard error: 13.1 on 52 degrees of freedom If you want to make things even easier you could use the self-starting model SSasympOrig as summary(fm1 - nls(bone ~ SSasympOrig(age, Asym, lrc), deer)) Formula: bone ~ SSasympOrig(age, Asym, lrc) Parameters: Estimate Std. Error t value Pr(|t|) Asym 115.5805 2.8436 40.65 2e-16 lrc -2.1302 0.1038 -20.52 2e-16 Residual standard error: 13.1 on 52 degrees of freedom The parameter lrc for that model is the natural logarithm of the rate constant 'c'. And everything goes fine, I understand the steps taken and so, but on the text Crawley says the model... and explained 84.6% of the total variation in bone lenght. I guess this is linked to the adjusted squared R for linear models, but I just can't find how to get it in this case... Actually it's just an R-squared, not an adjusted R-squared. It is being calculated as 1 - RSS/AdjSS where RSS is the residual sum of squares from the fitted model and AdjSS is the adjusted sum of squares or the sum of squares of the deviations of the observed responses from their mean. In this case the values are sum(resid(model2)^2) # RSS [1] 8929.143 with(deer, sum((bone - mean(bone))^2)) # AdjSS [1] 59007.99 so the value of R^2 from the formula is with(deer, 1 - sum(resid(model2)^2)/sum((bone - mean(bone))^2)) [1] 0.848679 I've tried in Statgraphics, and it plots the anova table and r^2 right the way, how could I do so in R? Yes, many software packages that fit nonlinear regression models do produce an anova table and an R^2 value for any model, even when that table and the R^2 value do not apply. (I'm assuming that the anova table is based on dividing the adjusted sum of squares into model and residual components so it is essentially the same calculation as above.) It would not be difficult to include these values in the summary output from an nls model in R but we don't because they could be nonsense. To see why we must examine what the R^2 should represent. First you should read what Bill Venables has to say on the general subject of analysis of variance for linear models. Use install.packages(fortunes); library(fortunes); fortune(curious) All the summaries like an anova table or an R^2 value are based on the comparison of two model fits - the model we just fit to the data and the corresponding trivial model. The trivial model is either an arbitrary constant or zero, depending on whether the model consisting of a constant only can be embedded in the model we have fit. If we can select values of the parameters that turn our model into a one-parameter model consisting of a constant then the comparison sum of squares is AdjSS as above because the parameter estimate in the trivial model is mean(response). If we can't embed the constant model in our model then the appropriate comparison sum of squares is the unadjusted sum of squares sum(response^2) Technically we can't embed the model for which the predictions are an arbitrary constant in this model because of that point at age == 0. No matter how we change the parameters 'a' and 'c' in a*(1-exp(-c*age)) we always predict zero at age == 0. Thus the only model with constant predictions that is embedded in this model is the model that predicts 0 for all ages so we should use the unadjusted sum of squares. However, if we ignore the point at age == 0 (which doesn't contribute any information to the model fit) then we can get a constant model by letting c go to +Inf. As c goes to +Inf the conditional estimate of a goes to mean(bone) so the constant model is embedded in the model we have fit. We don't include the R^2 or the anova table in the summary output for a nonlinear regression model because we can't tell if these are the appropriate formulas. Look at the model in
[R] adding elemens to a list
Hi, I have a list of 20 elements, each of them of variable length and with a structure like this: lasker[[1]][1:10,] Var1 Freq 1 1988-023 2 1988-031 3 1988-041 4 1988-052 5 1988-063 6 1988-071 7 1988-081 8 1988-091 9 1989-031 10 1989-041 How do I can insert in this list: 1988-01 0 1988-10 0 1988-11 0 1988-12 0 1989-01 0 1989-02 0 It's to say, the appropiate missing Year-Month row equal 0 Antonio PD: follows dput output of lasker[[1]] dput(lasker[[1]],control=all) structure(list(Var1 = structure(as.integer(c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 145, 146, 147, 148, 149, 150, 151, 152, 153, 154, 155, 156, 157, 158, 159, 160, 161, 162, 163, 164, 165, 166, 167, 168, 169, 170)), .Label = c(1988-02, 1988-03, 1988-04, 1988-05, 1988-06, 1988-07, 1988-08, 1988-09, 1989-03, 1989-04, 1989-05, 1989-07, 1989-08, 1989-09, 1989-11, 1990-01, 1990-02, 1990-03, 1990-04, 1990-05, 1990-06, 1990-07, 1990-08, 1990-10, 1990-11, 1991-01, 1991-03, 1991-04, 1991-05, 1991-06, 1991-09, 1991-11, 1991-12, 1992-01, 1992-03, 1992-04, 1992-05, 1992-06, 1992-07, 1992-08, 1992-09, 1992-10, 1992-11, 1992-12, 1993-01, 1993-02, 1993-03, 1993-04, 1993-05, 1993-07, 1993-08, 1993-09, 1993-10, 1993-12, 1994-02, 1994-03, 1994-04, 1994-05, 1994-07, 1994-08, 1994-09, 1994-10, 1994-11, 1994-12, 1995-04, 1995-05, 1995-06, 1995-07, 1995-08, 1995-09, 1995-10, 1995-11, 1996-02, 1996-04, 1996-05, 1996-06, 1996-07, 1996-08, 1996-09, 1996-11, 1996-12, 1997-01, 1997-02, 1997-04, 1997-06, 1997-07, 1997-08, 1997-09, 1997-10, 1998-01, 1998-02, 1998-04, 1998-05, 1998-06, 1998-08, 1998-09, 1998-10, 1998-11, 1998-12, 1999-03, 1999-05, 1999-06, 1999-07, 1999-08, 1999-09, 1999-10, 1999-11, 1999-12, 2000-01, 2000-02, 2000-05, 2000-06, 2000-07, 2000-08, 2000-09, 2000-10, 2000-11, 2000-12, 2001-02, 2001-03, 2001-04, 2001-05, 2001-07, 2001-08, 2001-09, 2001-10, 2001-11, 2001-12, 2002-01, 2002-02, 2002-04, 2002-05, 2002-06, 2002-08, 2002-09, 2002-10, 2002-11, 2002-12, 2003-01, 2003-02, 2003-03, 2003-04, 2003-05, 2003-06, 2003-07, 2003-08, 2003-09, 2003-12, 2004-01, 2004-02, 2004-03, 2004-04, 2004-05, 2004-06, 2004-07, 2004-08, 2004-09, 2004-10, 2004-11, 2004-12, 2005-01, 2005-02, 2005-03, 2005-04, 2005-06, 2005-07, 2005-08, 2005-09, 2005-10, 2005-12), class = factor), Freq = as.integer(c(3, 1, 1, 2, 3, 1, 1, 1, 1, 1, 2, 1, 2, 1, 1, 1, 2, 1, 1, 1, 3, 2, 1, 1, 3, 1, 1, 2, 1, 2, 1, 2, 1, 1, 2, 2, 3, 3, 2, 2, 1, 1, 1, 1, 1, 2, 2, 1, 2, 1, 2, 1, 2, 1, 1, 1, 2, 2, 1, 1, 2, 1, 2, 1, 2, 2, 2, 2, 1, 2, 1, 1, 2, 2, 1, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 3, 2, 2, 1, 1, 1, 1, 1, 2, 3, 2, 4, 2, 2, 1, 3, 1, 3, 2, 1, 1, 2, 1, 2, 1, 1, 2, 2, 2, 2, 3, 1, 1, 2, 1, 2, 1, 2, 2, 1, 1, 1, 1, 2, 2, 1, 1, 2, 2, 2, 1, 1, 1, 1, 1, 1, 3, 1, 2, 1, 2, 2, 1, 2, 1, 2, 3, 3, 1, 2, 1, 1, 2, 1, 1, 1, 1, 2, 2, 2, 1, 2, 2, 1, 2))), .Names = c(Var1, Freq), row.names = as.integer(c(NA, 170)), class = data.frame) __ 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] adding elemens to a list
If DF is the data frame in your post then turn it into a zoo object, convert that to class ts and then back to zoo replacing NA's with zero: library(zoo) z - zoo(DF$Freq, as.yearmon(as.Date(paste(DF$Var1, 01, sep = - zz - as.zoo(as.ts(z)) zz[is.na(zz)] - 0 Suggest you read the zoo vignette: vignette(zoo) On 11/26/06, antonio rodriguez [EMAIL PROTECTED] wrote: Hi, I have a list of 20 elements, each of them of variable length and with a structure like this: lasker[[1]][1:10,] Var1 Freq 1 1988-023 2 1988-031 3 1988-041 4 1988-052 5 1988-063 6 1988-071 7 1988-081 8 1988-091 9 1989-031 10 1989-041 How do I can insert in this list: 1988-01 0 1988-10 0 1988-11 0 1988-12 0 1989-01 0 1989-02 0 It's to say, the appropiate missing Year-Month row equal 0 Antonio PD: follows dput output of lasker[[1]] dput(lasker[[1]],control=all) structure(list(Var1 = structure(as.integer(c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 145, 146, 147, 148, 149, 150, 151, 152, 153, 154, 155, 156, 157, 158, 159, 160, 161, 162, 163, 164, 165, 166, 167, 168, 169, 170)), .Label = c(1988-02, 1988-03, 1988-04, 1988-05, 1988-06, 1988-07, 1988-08, 1988-09, 1989-03, 1989-04, 1989-05, 1989-07, 1989-08, 1989-09, 1989-11, 1990-01, 1990-02, 1990-03, 1990-04, 1990-05, 1990-06, 1990-07, 1990-08, 1990-10, 1990-11, 1991-01, 1991-03, 1991-04, 1991-05, 1991-06, 1991-09, 1991-11, 1991-12, 1992-01, 1992-03, 1992-04, 1992-05, 1992-06, 1992-07, 1992-08, 1992-09, 1992-10, 1992-11, 1992-12, 1993-01, 1993-02, 1993-03, 1993-04, 1993-05, 1993-07, 1993-08, 1993-09, 1993-10, 1993-12, 1994-02, 1994-03, 1994-04, 1994-05, 1994-07, 1994-08, 1994-09, 1994-10, 1994-11, 1994-12, 1995-04, 1995-05, 1995-06, 1995-07, 1995-08, 1995-09, 1995-10, 1995-11, 1996-02, 1996-04, 1996-05, 1996-06, 1996-07, 1996-08, 1996-09, 1996-11, 1996-12, 1997-01, 1997-02, 1997-04, 1997-06, 1997-07, 1997-08, 1997-09, 1997-10, 1998-01, 1998-02, 1998-04, 1998-05, 1998-06, 1998-08, 1998-09, 1998-10, 1998-11, 1998-12, 1999-03, 1999-05, 1999-06, 1999-07, 1999-08, 1999-09, 1999-10, 1999-11, 1999-12, 2000-01, 2000-02, 2000-05, 2000-06, 2000-07, 2000-08, 2000-09, 2000-10, 2000-11, 2000-12, 2001-02, 2001-03, 2001-04, 2001-05, 2001-07, 2001-08, 2001-09, 2001-10, 2001-11, 2001-12, 2002-01, 2002-02, 2002-04, 2002-05, 2002-06, 2002-08, 2002-09, 2002-10, 2002-11, 2002-12, 2003-01, 2003-02, 2003-03, 2003-04, 2003-05, 2003-06, 2003-07, 2003-08, 2003-09, 2003-12, 2004-01, 2004-02, 2004-03, 2004-04, 2004-05, 2004-06, 2004-07, 2004-08, 2004-09, 2004-10, 2004-11, 2004-12, 2005-01, 2005-02, 2005-03, 2005-04, 2005-06, 2005-07, 2005-08, 2005-09, 2005-10, 2005-12), class = factor), Freq = as.integer(c(3, 1, 1, 2, 3, 1, 1, 1, 1, 1, 2, 1, 2, 1, 1, 1, 2, 1, 1, 1, 3, 2, 1, 1, 3, 1, 1, 2, 1, 2, 1, 2, 1, 1, 2, 2, 3, 3, 2, 2, 1, 1, 1, 1, 1, 2, 2, 1, 2, 1, 2, 1, 2, 1, 1, 1, 2, 2, 1, 1, 2, 1, 2, 1, 2, 2, 2, 2, 1, 2, 1, 1, 2, 2, 1, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 3, 2, 2, 1, 1, 1, 1, 1, 2, 3, 2, 4, 2, 2, 1, 3, 1, 3, 2, 1, 1, 2, 1, 2, 1, 1, 2, 2, 2, 2, 3, 1, 1, 2, 1, 2, 1, 2, 2, 1, 1, 1, 1, 2, 2, 1, 1, 2, 2, 2, 1, 1, 1, 1, 1, 1, 3, 1, 2, 1, 2, 2, 1, 2, 1, 2, 3, 3, 1, 2, 1, 1, 2, 1, 1, 1, 1, 2, 2, 2, 1, 2, 2, 1, 2))), .Names = c(Var1, Freq), row.names = as.integer(c(NA, 170)), class = data.frame) __ 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] adding elemens to a list
Gabor Grothendieck escribió: If DF is the data frame in your post then turn it into a zoo object, convert that to class ts and then back to zoo replacing NA's with zero: library(zoo) z - zoo(DF$Freq, as.yearmon(as.Date(paste(DF$Var1, 01, sep = - zz - as.zoo(as.ts(z)) zz[is.na(zz)] - 0 Suggest you read the zoo vignette: vignette(zoo) Gabor, Once again, thanks! BR Antonio On 11/26/06, antonio rodriguez [EMAIL PROTECTED] wrote: Hi, I have a list of 20 elements, each of them of variable length and with a structure like this: lasker[[1]][1:10,] Var1 Freq 1 1988-023 2 1988-031 3 1988-041 4 1988-052 5 1988-063 6 1988-071 7 1988-081 8 1988-091 9 1989-031 10 1989-041 How do I can insert in this list: 1988-01 0 1988-10 0 1988-11 0 1988-12 0 1989-01 0 1989-02 0 It's to say, the appropiate missing Year-Month row equal 0 Antonio PD: follows dput output of lasker[[1]] dput(lasker[[1]],control=all) structure(list(Var1 = structure(as.integer(c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 145, 146, 147, 148, 149, 150, 151, 152, 153, 154, 155, 156, 157, 158, 159, 160, 161, 162, 163, 164, 165, 166, 167, 168, 169, 170)), .Label = c(1988-02, 1988-03, 1988-04, 1988-05, 1988-06, 1988-07, 1988-08, 1988-09, 1989-03, 1989-04, 1989-05, 1989-07, 1989-08, 1989-09, 1989-11, 1990-01, 1990-02, 1990-03, 1990-04, 1990-05, 1990-06, 1990-07, 1990-08, 1990-10, 1990-11, 1991-01, 1991-03, 1991-04, 1991-05, 1991-06, 1991-09, 1991-11, 1991-12, 1992-01, 1992-03, 1992-04, 1992-05, 1992-06, 1992-07, 1992-08, 1992-09, 1992-10, 1992-11, 1992-12, 1993-01, 1993-02, 1993-03, 1993-04, 1993-05, 1993-07, 1993-08, 1993-09, 1993-10, 1993-12, 1994-02, 1994-03, 1994-04, 1994-05, 1994-07, 1994-08, 1994-09, 1994-10, 1994-11, 1994-12, 1995-04, 1995-05, 1995-06, 1995-07, 1995-08, 1995-09, 1995-10, 1995-11, 1996-02, 1996-04, 1996-05, 1996-06, 1996-07, 1996-08, 1996-09, 1996-11, 1996-12, 1997-01, 1997-02, 1997-04, 1997-06, 1997-07, 1997-08, 1997-09, 1997-10, 1998-01, 1998-02, 1998-04, 1998-05, 1998-06, 1998-08, 1998-09, 1998-10, 1998-11, 1998-12, 1999-03, 1999-05, 1999-06, 1999-07, 1999-08, 1999-09, 1999-10, 1999-11, 1999-12, 2000-01, 2000-02, 2000-05, 2000-06, 2000-07, 2000-08, 2000-09, 2000-10, 2000-11, 2000-12, 2001-02, 2001-03, 2001-04, 2001-05, 2001-07, 2001-08, 2001-09, 2001-10, 2001-11, 2001-12, 2002-01, 2002-02, 2002-04, 2002-05, 2002-06, 2002-08, 2002-09, 2002-10, 2002-11, 2002-12, 2003-01, 2003-02, 2003-03, 2003-04, 2003-05, 2003-06, 2003-07, 2003-08, 2003-09, 2003-12, 2004-01, 2004-02, 2004-03, 2004-04, 2004-05, 2004-06, 2004-07, 2004-08, 2004-09, 2004-10, 2004-11, 2004-12, 2005-01, 2005-02, 2005-03, 2005-04, 2005-06, 2005-07, 2005-08, 2005-09, 2005-10, 2005-12), class = factor), Freq = as.integer(c(3, 1, 1, 2, 3, 1, 1, 1, 1, 1, 2, 1, 2, 1, 1, 1, 2, 1, 1, 1, 3, 2, 1, 1, 3, 1, 1, 2, 1, 2, 1, 2, 1, 1, 2, 2, 3, 3, 2, 2, 1, 1, 1, 1, 1, 2, 2, 1, 2, 1, 2, 1, 2, 1, 1, 1, 2, 2, 1, 1, 2, 1, 2, 1, 2, 2, 2, 2, 1, 2, 1, 1, 2, 2, 1, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 3, 2, 2, 1, 1, 1, 1, 1, 2, 3, 2, 4, 2, 2, 1, 3, 1, 3, 2, 1, 1, 2, 1, 2, 1, 1, 2, 2, 2, 2, 3, 1, 1, 2, 1, 2, 1, 2, 2, 1, 1, 1, 1, 2, 2, 1, 1, 2, 2, 2, 1, 1, 1, 1, 1, 1, 3, 1, 2, 1, 2, 2, 1, 2, 1, 2, 3, 3, 1, 2, 1, 1, 2, 1, 1, 1, 1, 2, 2, 2, 1, 2, 2, 1, 2))), .Names = c(Var1, Freq), row.names = as.integer(c(NA, 170)), class = data.frame) __ 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] Fitting mixed-effects models with lme with fixed error term variances
On 11/26/06, Gregor Gorjanc [EMAIL PROTECTED] wrote: Hello! Douglas Bates wrote: ... Not quite. There is now a capability in lmer to specify starting estimates for the relative precision matrices which means that you can use the estimates from one fit as starting estimates for another fit. It looks like fm1 - lmer(...) fm2 - lmer(y ~ x + (...), start = [EMAIL PROTECTED]) I should say that in my experience this has not been as effective as I had hoped it would be. What I see in the optimization is that the first few iterations reduce the deviance quite quickly but the majority of the time is spent refining the estimates near the optimum. So, for example, if it took 40 iterations to converge from the rather crude starting estimates calculated within lmer it might instead take 35 iterations if you give it good starting estimates. Yes, I also have the same experience with other programs, say VCE[1]. What I was hopping for was just solutions from linear mixed model i.e. either via GLS or PLS representations and no optimization if values for variance components are passed as input - there should be a way to distinguish that from just passing starting values.. The PLS representation in lmer is in terms of the relative variance-covariance matrices of the random effects where relative means relative to s^2. (In fact, the Omega slot contains the relative What exactly is s^2? Sorry, I confused a parameter and its estimate. The variance-covariance matrices are relative to the variance of the per-observation noise which is assumed to be Gaussian with mean 0 and variance-covariance matrix sigma^2 * I. That is, they are relative to sigma^2 which is estimated by s^2 the penalized residual sum of squares divided by the number of observations. precision matrices (inverses of the relative variance matrices) but its the same idea. All variance components are expressed relative to s^2.) If it is sufficient to fix these then you can easily do that. Just follow the code in lmer until it gets to .Call(mer_ECMEsteps, mer, cv$niterEM, cv$EMverbose) LMEoptimize(mer) - cv return(new(lmer, mer, and replace the first two lines by .Call(mer_factor, mer, PACKAGE = lme4) Hmm. Sounds simple. I will try this. Do you find any value in adding this as an option in lmer or perhaps new? function, say lmeSol() or something similar? I'll have to think about it. It is the usual problem of determining how to specify this option, documenting it, explaining it and maintaining it along with the rest of the code. A side-effect of performing the factorization is to evaluate the ML and REML deviances and store the result in the deviance slot. The follow the code in lmer part will get a bit tricky because of the number of functions that are hidden in the lme4 namespace but I'm sure you know how to get around that. Thank you very much! Gregor __ 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] Fixed zeros in tables
Hello Andrew Robinson and R-List Thanks, Andrew, but this does not work. Puttting zero weights on structural zeros, one's elsewhere in glm() does not deliver the appropriate expected cell counts loglm() provides as one would expect. If you run the code provided below, you'll see loglm() delivers zero cell counts with loglm() but using glm() with the weights you suggest, the expected cell counts are not zero. Still hoping for a resolution. Andrew Criswell Associate Professor Hedmark University Postboks 104, Rena 2510, NORWAY -- ## Fienberg, The Analysis of Cross-Classified Contingency Tables, 2nd ed., p.148. ## Results from survey of teenagers regarding their health concerns. health - data.frame(expand.grid(CONCERNS = c(sex, menstral, healthy, nothing), AGE = c(12-15, 16-17), GENDER = c(male, female)), COUNT= c(4, 0, 42, 57, 2, 0, 7, 20, 9, 4, 19, 71, 7, 8, 10, 21), WEIGHTS = c(1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1)) health.tbl - xtabs(COUNT ~ CONCERNS + AGE + GENDER, data = health) zeros - data.frame(expand.grid(CONCERNS = c(sex, menstral, healthy, nothing), AGE = c(12-15, 16-17), GENDER = c(male, female)), COUNT= c(1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1)) zeros - xtabs(COUNT ~ CONCERNS + AGE + GENDER, data = zeros) library(MASS) fm.1 - loglm(~ CONCERNS + AGE + GENDER, data = health.tbl, start = zeros, fitted = TRUE) fm.1; round(fm.1$fitted, 1) fm.3 - glm(COUNT ~ CONCERNS + AGE + GENDER, data = health, weights = health$WEIGHTS, family = poisson) fm.3.fit - data.frame(expand.grid(CONCERNS = c(sex, menstral, healthy, nothing), AGE = c(12-15, 16-17), GENDER = c(male, female)), COUNT= fm.3$fitted) round(xtabs(COUNT ~ CONCERNS + AGE + GENDER, data = fm.3.fit), 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] Fixed zeros in tables
Hello Andrew, I'm not sure how this is a problem. You can multiply the fitted values by the weights again, if you wish, or ignore the structural zeros altogether. The other predicted values all seem to me to be the same. I think that is a feature. Sometimes the goal for glm with structural zeros is to estimate those values that are missing. You are at liberty to ignore them. Cheers Andrew On Mon, Nov 27, 2006 at 02:47:01AM +0700, A.R. Criswell wrote: Hello Andrew Robinson and R-List Thanks, Andrew, but this does not work. Puttting zero weights on structural zeros, one's elsewhere in glm() does not deliver the appropriate expected cell counts loglm() provides as one would expect. If you run the code provided below, you'll see loglm() delivers zero cell counts with loglm() but using glm() with the weights you suggest, the expected cell counts are not zero. Still hoping for a resolution. Andrew Criswell Associate Professor Hedmark University Postboks 104, Rena 2510, NORWAY -- ## Fienberg, The Analysis of Cross-Classified Contingency Tables, 2nd ed., p.148. ## Results from survey of teenagers regarding their health concerns. health - data.frame(expand.grid(CONCERNS = c(sex, menstral, healthy, nothing), AGE = c(12-15, 16-17), GENDER = c(male, female)), COUNT= c(4, 0, 42, 57, 2, 0, 7, 20, 9, 4, 19, 71, 7, 8, 10, 21), WEIGHTS = c(1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1)) health.tbl - xtabs(COUNT ~ CONCERNS + AGE + GENDER, data = health) zeros - data.frame(expand.grid(CONCERNS = c(sex, menstral, healthy, nothing), AGE = c(12-15, 16-17), GENDER = c(male, female)), COUNT= c(1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1)) zeros - xtabs(COUNT ~ CONCERNS + AGE + GENDER, data = zeros) library(MASS) fm.1 - loglm(~ CONCERNS + AGE + GENDER, data = health.tbl, start = zeros, fitted = TRUE) fm.1; round(fm.1$fitted, 1) fm.3 - glm(COUNT ~ CONCERNS + AGE + GENDER, data = health, weights = health$WEIGHTS, family = poisson) fm.3.fit - data.frame(expand.grid(CONCERNS = c(sex, menstral, healthy, nothing), AGE = c(12-15, 16-17), GENDER = c(male, female)), COUNT= fm.3$fitted) round(xtabs(COUNT ~ CONCERNS + AGE + GENDER, data = fm.3.fit), 1) -- Andrew Robinson Department of Mathematics and StatisticsTel: +61-3-8344-9763 University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599 http://www.ms.unimelb.edu.au/~andrewpr http://blogs.mbs.edu/fishing-in-the-bay/ __ 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] DNA dotmatrix window argument
Hi all, I just compiled package dna v. 0.2 from source and am having some difficulty with the window argument of the dotmatrix function. Some test code; say I'm simulating some nucleotide sequences: sequence1-trunc(runif(200,0,4)) sequence2-trunc(runif(200,0,4)) I embed a nice 100 element alignment: sequence1[7:107] - sequence2[60:160] Plotting this with a window size of 1 produces a very dense plot because of the small grammar size. dotmatrix(rbind(sequence1, sequence2), window=1) However, when I plot it with a larger window (generally between 6 and 100) dotmatrix(rbind(sequence1, sequence2), window=10) I get no matches (no dots). Perhaps I'm confusing the meaning of window in this function with another type of window used to fight this random-match issue. Any help would be appreciated. Jeff. http://www.nd.edu/~jspies/ version() _ platform powerpc-apple-darwin8.7.0 arch powerpc os darwin8.7.0 system powerpc, darwin8.7.0 status major 2 minor 4.0 year 2006 month 10 day03 svn rev39566 language R version.string R version 2.4.0 (2006-10-03) [[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] axis in an image () plot
On Nov 26, 2006, at 5:23 PM, jim holtman wrote: The phrase at the bottom (what is the problem that you are trying to solve) is one that I have used for the last 20 years at work and I have it as part of my signature line since a lot of people know me from that phrase. It is one of the first questions that I ask a project when I am doing a review for them. Another way of say it is tell me what you want to do, not how you want to do it. Thanks again, Jim. Please, allow me to try to use this method from now on! Best, Ricardo __ 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] GLM and LM singularities
Hi- I'm wrestling with some of my data apparently not being called into a GLM or an LM. I'm looking at factors affecting fish annual catch rates (ie. CPUE) over 30 years. Two of the factors I'm using are sea surface temperature and sea surface temperature anomaly. A small sample of my data is below: CPUE Year Vessel_ID Base_Port Boat_Lgth Planing SST Anomaly 0.127 1976 2 BOI 40 Planing 19.4 -0.308 0.059 1976 30 BOI 40 Displacement 19.4 -0.308 0.03 1977 46 BOI 45 Displacement 18.5 -1.008 0.169231 1977 2 BOI 40 Planing 18.5 -1.008 0.044118 1977 19 BOI 34 Planing 18.5 -1.008 0.114286 1977 29 WHANGAROA 42 Displacement 18.5 -1.008 Have defined Year, Vessel_ID, Base_Port, Boat_Lgth, Planing as factors, and CPUE, SST, Anomaly are continuous variables. The formula I'm calling is: glm(formula = log(CPUE) ~ Year + Vessel_ID + SST, data = marlin). A summary of the glm includes the following output: Coefficients: (1 not defined because of singularities) Estimate Std. Error t value Pr(|t|) Vessel_ID80 -0.540930.23380 -2.314 0.021373 * Vessel_ID83 -0.364990.20844 -1.751 0.080977 . Vessel_ID84 -0.233860.19098 -1.225 0.221718 SST NA NA NA NA Can someone explain the output Coefficients: (1 not defined because of singularities)? What does this mean? I'm guessing that something to do with my SST data is causing a singularity but I don't know where to go from there? Have inspected various R discussions on this, but I haven't found the information I need. Many thanks, Tim Sippel (MSc) School of Biological Sciences Auckland University Private Bag 92019 Auckland 1020 New Zealand +64-9-373-7599 ext. 84589 (work) +64-9-373-7668 (Fax) +64-21-593-001 (mobile) [[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] Questions about generating samples in R
split - sample(2,nrow(dataframe),replace=T,prob=c(0.04,0.96)) dataframe[split==1,] # 200 dataframe[split==2,] # 4800 regards, christian ?sample should tell you what you need to know. On 26/11/06, Alexander Geisler [EMAIL PROTECTED] wrote: Hello! I have a data set with 8 columns and in about 5000 rows. What I want to do is to generate samples of this data set. Samples of a special size, as example 200. What is the easiest way to do this? No special things are needed, only the random selection of 200 rows of the data set. Thanks Alex -- Alexander Geisler * Kaltenbach 151 * A-6272 Kaltenbach email: [EMAIL PROTECTED] | [EMAIL PROTECTED] phone: +43 650 / 811 61 90 | skpye: al1405ex __ 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] Nonlinear statistical modeling -- a comparison of R and AD Model Builder
Douglas Bates wrote: snip Don't you find it somewhat disingenuous that you publish a comparison between the AD Model Builder software that you sell and R - a comparison that shows a tremendous advantage for your software - and then you write I am not proficient in R? I think there is a misunderstanding here. I did not pick this example. As I said it was undertaken by Schnute and his colleagues. I had nothing to do with it except of course to sell them my software. However as I stated at that time Splus could not run the model without crashing after a time so that no comparison was possible. I was aware of those results and decided to use R to complete the comparison. I don't see why Schnute would want to unfairly promote my software. I believe he was simply looking for the best tool for the job he had in mind. Had you been proficient in R you might have known about the symbolic differentiation capabilities, specifically the deriv function, that have been part of the S language since the late 1980s. I believe that the 'AD' in AD model builder stands for automatic differentiation, which is actually something that John Chambers and I discussed at length when we were developing nonlinear modeling methods for S. In the end we went with symbolic differentiation rather than automatic differentiation because we felt that symbolic was more flexible. Yes I am aware of the symbolic differentiation capabilities. I have checked out the deriv function and it does not seem to be capable of calculating derivatives for a model of this complexity in an efficient manner. Of course I could be wrong. There is a paper by Andreas Griewank (whose title I have forgotten but perhaps some list member recalls) around 1990 where he compares symbolic and automatic differentiation for a simple model of an oil reservoir. He demonstrates quite decisively that symbolic differentiation is not the way to go. This is not to say that automatic differentiation isn't a perfectly legitimate technique. However, my recollection is that it would have required extensive revisions to the arithmetic expression evaluator, which is already very tricky code because of the recycling rule and the desire to shield users from knowledge of the internal representations and such details as whether you are using logical or integer or double precision operands or a combination. If you want to see these details you can, of course, examine the source code. I don't believe we would have the opportunity to examine how you implemented automatic differentiation. I must also agree with Spencer Graves that when I start reading a description of a nonlinear model with over 100 parameters, the example that you chose, I immediately start thinking of nonlinear mixed effects models. In my experience the only way in which a nonlinear model ends up with that number of parameters is through applying an underlying model with a low number of parameters to various groups within the data. Table 2 in the Schnute et al. paper to which you make reference states that the number of parameters in the model is T + A + 5 where T is the number of years of data and A is the number of age classes. To me that looks a lot like a nonlinear mixed effects model. I agree that this makes a good nonlinear random effects example. Of course 10 years ago AD Model Builder did not have that capability. It now does and my colleague Hans Skaug has modified the code to incorporate random effects. I believe the model converges in a few minutes. He will report the results and hopefully they can be compared to nlme or any other software in R which can carry out the calculations. Also, your choice of subject heading for your message seems deliberatively provocative. You seem to be implying that you are discussing a comparisons of AD Model Builder and R on all aspects of nonlinear statistical modeling but you are only discussing one comparison on simulated data using a model from the applications area for which you wrote AD Model Builder. Then you follow up by saying I am not proficient in R and your results for R are from applying code that someone else gave you. It seems that ADMB had a bit of a home-field advantage in this particular comparison. I don't quite get your point. Of course I am only going to present examples where I believe ADMB is (far) superior to R. Otherwise I would just be wasting everyones time. ADMB is much more narrowly focused than R. I think that people can examine the arguments and make up their own minds. I view nonlinear statistical modeling differently. I have had a bit of experience in the area and I find that I want to examine the data carefully, usually through plots, before I embark on fitting complicated models. I like to have some assurance that the model makes sense in the context of the data. (In your example you don't need to worry about appropriateness of the model because the data were
[R] Comment on Sweave
Sorry. broken Engrish. I tried to give comment on Sweave by force. rawcode=T is necessary for an option(If you need comment). try wget http://r.nakama.ne.jp/CommentOnSweave/SweaveProfile R_PROFILE=SweaveProfile R CMD Sweave hoge.Rnw -- EI-JI Nakama [EMAIL PROTECTED] \u4e2d\u9593\u6804\u6cbb [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.
[R] Help with response CIs for lme
Hi,
Can someone please offer a procedure for going from CIs produced by
intervals.lme() to fixed-effects response CIs.
Here's a simple example:
library(mlmRev)
library(nlme)
hsb.lme - lme(mAch ~ minrty * sector, random = ~ 1 | cses, Hsb82)
(intervals(hsb.lme))
(hsb.new - data.frame
minrty = rep(c('No', 'Yes'), 2),
sector = rep(c('Public', 'Catholic'), each = 2)))
cbind(hsb.new, predict(hsb.lme, hsb.new, level = 0))
Is the following correct (I know from the previous command that the
estimate is correct)?
cbind(hsb.new, rbind(hsb.int[[1]][1,], hsb.int[[1]][1,]+hsb.int[[1]]
[2,], hsb.int[[1]][1,]+hsb.int[[1]][3,], hsb.int[[1]][1,]+hsb.int[[1]]
[2,] + hsb.int[[1]][3,] + hsb.int[[1]][4,]))
If so, is there an easier way to write it?
_
Professor Michael Kubovy
University of Virginia
Department of Psychology
USPS: P.O.Box 400400Charlottesville, VA 22904-4400
Parcels:Room 102Gilmer Hall
McCormick RoadCharlottesville, VA 22903
Office:B011+1-434-982-4729
Lab:B019+1-434-982-4751
Fax:+1-434-982-4766
WWW:http://www.people.virginia.edu/~mk9y/
__
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] x axis on sciplot
You will need a numeric axis, not a factor axis, and control of both
the at= and labels= of the axis. Here is an example using your
data and the lattice xyplot function.
try$visit.position - try$visit
levels(try$visit.position)
levels(try$visit.position)[1] - -20
levels(try$visit.position)
try$visit.position - as.numeric(as.character(try$visit.position))
try$visit.position
try$fake.ID - factor(rep(1:6, rep(5,6)))
xyplot(variable ~ visit.position,
group=interaction(group, fake.ID)[, drop=TRUE],
data=try, type=l,
scales=list(
at=unique(try$visit.position),
labels=levels(try$visit)),
auto.key=TRUE)
xyplot(variable ~ visit.position | group,
group=fake.ID,
data=try, type=l,
scales=list(
at=unique(try$visit.position),
labels=levels(try$visit)),
panel=function(...) {
panel.abline(v=-10, lty=2, col=gray50)
panel.superpose(...)
},
auto.key=TRUE)
I had some trouble with your email. It folded at the
column where you said
class=
data.frame
and therefore the class was interpreted as \ndata.frame,
and was not seen as a data.frame.
I placed a vertical boundary in the plot between the screening time
and the time=0.
Rich
__
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-sig-Geo] plot() and Jpeg() increase font size and resolution
Thanks to Edzer and Roger,
I can now plot with increased font sizes. However, jpeg still does not
reproduce these, nor does it show up in high quality. What I would like
to do is produce some highresolution jpegs.
Any help would be appreciated
Thanx
Herry
R2.4 on Mandriva 10.2 linux.
Dr Alexander Herr
Spatial and statistical analyst
CSIRO, Sustainable Ecosystems
Davies Laboratory,
University Drive, Douglas, QLD 4814
Private Mail Bag, Aitkenvale, QLD 4814
Phone/www
(07) 4753 8510; 4753 8650(fax)
Home: http://herry.ausbats.org.au
Webadmin ABS: http://ausbats.org.au
Sustainable Ecosystems: http://www.cse.csiro.au/
-Original Message-
From: Edzer J. Pebesma [mailto:[EMAIL PROTECTED]
Sent: Friday, November 24, 2006 5:37 PM
To: Herr, Alexander Herr - Herry (CSE, Townsville)
Cc: r-help@stat.math.ethz.ch; r-sig-geo@stat.math.ethz.ch
Subject: Re: [R-sig-Geo] plot() and Jpeg() increase font size and
resolution
plotting variograms in gstat is done through xyplot in lattice; you'll
find where it gets it's defaults by
library(lattice)
trellis.par.get()
?trellis.par.set
--
Edzer
[EMAIL PROTECTED] wrote:
Dear list,
I am having troubles increasing the fontize when plotting a
variogram{gstat} and its model (vgm) with plot and using jpeg(). Also
the resolution in the jpeg call does not work. I am using R2.4 on
Mandriva 10.2 linux.
I can change fontsize with cex.axis in a normal plot, so I presume it
has to do with plotting the variogram model. Any help on how to
increase the font size and resolution would be appreciated?
R-calls:
v1-variogram(log(z)~x+y, loc=coordinates(shp1),data=shp1)
m1-vgm(0.0175, Gau, 20,0.052)
jpeg(file=LOTPLAN_variogram_mod.jpg, bg=white, res=300,
pointsize=16, width=1200, height=1200, quality=100)
plot(v1,plot.number=F, model=m1, ylim=c(0.04,0.08), col=black,
cex.axis=1.5)
dev.off()
Thanks Herry
Dr Alexander Herr
Spatial and statistical analyst
CSIRO, Sustainable Ecosystems
Davies Laboratory,
University Drive, Douglas, QLD 4814
Private Mail Bag, Aitkenvale, QLD 4814
Phone/www
(07) 4753 8510; 4753 8650(fax)
Home: http://herry.ausbats.org.au
Webadmin ABS: http://ausbats.org.au
Sustainable Ecosystems: http://www.cse.csiro.au/
___
R-sig-Geo mailing list
R-sig-Geo@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
__
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-sig-Geo] plot() and Jpeg() increase font size and resolution
[EMAIL PROTECTED] wrote:
Thanks to Edzer and Roger,
I can now plot with increased font sizes. However, jpeg still does not
reproduce these, nor does it show up in high quality. What I would like
to do is produce some highresolution jpegs.
Please specify the commands that you used, preferrably such that we can
reproduce them, and specify which worked and which didn't the way you
expected. Also please give the output of sessionInfo(), and tell us the
platform you work on. Have you tried producing png, did that work?
--
Edzer
Any help would be appreciated
Thanx
Herry
R2.4 on Mandriva 10.2 linux.
Dr Alexander Herr
Spatial and statistical analyst
CSIRO, Sustainable Ecosystems
Davies Laboratory,
University Drive, Douglas, QLD 4814
Private Mail Bag, Aitkenvale, QLD 4814
Phone/www
(07) 4753 8510; 4753 8650(fax)
Home: http://herry.ausbats.org.au
Webadmin ABS: http://ausbats.org.au
Sustainable Ecosystems: http://www.cse.csiro.au/
-Original Message-
From: Edzer J. Pebesma [mailto:[EMAIL PROTECTED]
Sent: Friday, November 24, 2006 5:37 PM
To: Herr, Alexander Herr - Herry (CSE, Townsville)
Cc: r-help@stat.math.ethz.ch; r-sig-geo@stat.math.ethz.ch
Subject: Re: [R-sig-Geo] plot() and Jpeg() increase font size and
resolution
plotting variograms in gstat is done through xyplot in lattice; you'll
find where it gets it's defaults by
library(lattice)
trellis.par.get()
?trellis.par.set
--
Edzer
[EMAIL PROTECTED] wrote:
Dear list,
I am having troubles increasing the fontize when plotting a
variogram{gstat} and its model (vgm) with plot and using jpeg(). Also
the resolution in the jpeg call does not work. I am using R2.4 on
Mandriva 10.2 linux.
I can change fontsize with cex.axis in a normal plot, so I presume it
has to do with plotting the variogram model. Any help on how to
increase the font size and resolution would be appreciated?
R-calls:
v1-variogram(log(z)~x+y, loc=coordinates(shp1),data=shp1)
m1-vgm(0.0175, Gau, 20,0.052)
jpeg(file=LOTPLAN_variogram_mod.jpg, bg=white, res=300,
pointsize=16, width=1200, height=1200, quality=100)
plot(v1,plot.number=F, model=m1, ylim=c(0.04,0.08), col=black,
cex.axis=1.5)
dev.off()
Thanks Herry
Dr Alexander Herr
Spatial and statistical analyst
CSIRO, Sustainable Ecosystems
Davies Laboratory,
University Drive, Douglas, QLD 4814
Private Mail Bag, Aitkenvale, QLD 4814
Phone/www
(07) 4753 8510; 4753 8650(fax)
Home: http://herry.ausbats.org.au
Webadmin ABS: http://ausbats.org.au
Sustainable Ecosystems: http://www.cse.csiro.au/
___
R-sig-Geo mailing list
R-sig-Geo@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
___
R-sig-Geo mailing list
R-sig-Geo@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
__
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 to simulate data from copulas?
Dears, I am writing sice having question about your R package copula. I am using this package to fit my data with some of the 22 copulas mentioned in Nelsen 1999. However, when trying random number generator I have got some difficulties. For example when I have familie number A13 where cupola is equal to C(u,v) = exp(1-((1-log(u))^theta+(1-log(v))^theta-1)^(1/theta)) devuv = (exp(1 - (-1 + (1 - log(u))^theta + (1 - log(v))^theta)^theta^(-1))*(1 - log(u))^(-1 + theta)*(-1 + theta + (-1 + (1 - log(u))^theta + (1 - log(v))^theta)^theta^(-1))*(-1 + (1 - log(u))^theta + (1 - log(v))^theta)^(-2 + theta^(-1))*(1 - log(v))^(-1 + theta))/(u*v) devu = (exp(1 - (-1 + (1 - log(u))^theta + (1 - log(v))^theta)^theta^(-1))*(1 - log(u))^(-1 + theta)*(-1 + (1 - log(u))^theta + (1 - log(v))^theta)^(-1 + theta^(-1)))/u phi = (1-log(t))^theta-1 devphi = - ((theta*(1 - log(t))^(-1 + theta))/t) K(t) = t - ((1-log(t))^theta-1) / (-((theta*(1 - log(t))^(-1 + theta))/t)) I am always meeting some difficulies to random generate data. Or to make a program code to count inverse function of my phi which should help in the algorithm to generate data. It would be grateful if you could help me. Thanks a lot. -- Mgr. Rastislav Matú PhD. Student Department of Land and Water Resources Management Faculty of Civil Engineering Slovak University of Technology Blok C, 12. posch., Radlinského 11, 813 68 Bratislava 15, EUROPE BUREAU: Tel.: +421 2 59274 622, fax: +421 2 52923575 PERSONAL: +421 908 623450 E-mail: [EMAIL PROTECTED] www.kvhk.sk [[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.
