[R] Piecewise Regression with a known slope
Hey, all. I'm working on a data set with a broken stick linear regression where I know one of the two slopes. It is a negative linear function until the line intersects with the x-axis, at which point it becomes 0. It is not a nonlinear asymptotic function, and, indeed, using negative exponential or logistic types of fits as an approximation has tended to lead to an under or overestimation of values. I am also very interested to know just what the breakpoint in the data is. Essentially if xpsi y = a + bx + error, where b is negative else y=0+error and I want to know the value of psi, as well as a and b. The error is also most likely gamma distributed, as values0 are not possible (nutrient concentrations). I have attempted to use the segmented package, specifying something close to the visually estimated breakpoint as the value of psi, but it continues to return fits with a breakpoint that occurs somewhere in the middle of the linear portion of the line, and the slope and intercept of the second half of the regression is not 0. Is there either a package that exists that will allow me to estimate such a model, or a function that I can use for optim or nlm (I admit, I am a rank novice at coding such functions). Thanks so much! -Jarrett __ 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] Going from Coda Files to R2WinBUGS
I'm currently using JAGS as my Bayesian program of choice due to
working off of an older mac running OSX. I'd like to utilize some of
the functions from R2WinBUGS, however. As such, I'm attempting to
write something to go from coda output dumped by JAGS into the bugs
object format - I've looked for functions that will convert from an
mcmc object to a bugs object, but have had no luck as of yet. I've
attempted to adapt Yu-Sung Su's method over at http://
yusung.blogspot.com/2007/01/analyzing-coda-statistics-produced-
by.html . However, whenever I run it for a coda file set generated
by jags, I get the following error:
Error in if (trans[i] == log) { : missing value where TRUE/FALSE
needed
This is for a run whose jags.ind is as follows
cypraea.effect 1 1
intercept 10001 2
sponge.sd 20001 3
deviance30001 4
When I debuged R2WinBUGS:::monitor, which is where the whole thing
borks, I found that trans is as follows
cypraea.effect intercept sponge.sd deviance
NA NA
NA
And the error comes when first looking at intercept, which is NA, not
. I am somewhat unclear as to why this is so.
The code for the method is as follows. Any thoughts would be greatly
appreciated, and if this works out, feel free to use it yourself!
Could be quite useful!
-Jarrett
#note, the test run was something along the lines of c.bugs-coda2bugs
(n.burnin=1000)
coda2bugs-function(codafile=jags.out, indexfile=jags.ind,
n.chains=1,
n.iter=NA, n.burnin=NA, n.thin=1,
DIC=FALSE, file.rm=T, ...){
require(R2WinBUGS)
#first, split up the coda file for R2WinBUGS
codaSplit(codafile, indexfile)
#get the parameter names
index.table-read.table(indexfile)
varNames-as.vector(index.table[,1])
#determine the n.iter
if(is.na(n.iter)){n.iter-index.table[1,3]}
#you will need to put the n.burnin in yourself
#for the cypraea example, it is 1000
bugs.fit - R2WinBUGS:::bugs.sims(varNames, n.chains=n.chains,
n.iter=n.iter, n.burnin=n.burnin,
n.thin=n.thin, DIC = DIC)
class(bugs.fit)-bugs
bugs.fit$isDIC - FALSE
#clean up the new coda files
if(file.rm==T){
file.remove(codaIndex.txt)
for(i in rownames(index.table)){
file.remove(paste(coda,i,.txt,sep=))
}
}
return(bugs.fit)
}
codaSplit-function(codafile=jags.out, indexfile=jags.ind){
index.table-read.table(indexfile)
write.table(index.table, codaIndex.txt, quote=F, row.names=F,
col.names=F, sep=\t)
coda.table-read.table(codafile)
#write the new coda files
for(i in rownames(index.table)){
new.file=paste(coda,i,.txt, sep=)
new.out-coda.table[index.table[i,2]:index.table[i,3],]
write.table(new.out, new.file, row.names=F, col.names=F,
sep=\t)
}
}
__
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] glms with poisson and negative binomial errors
A reviewer recently remarked to me that, due to my data being constrained to not fall below zero, a generalized linear model with a negative binomial error (or poisson) with a log link would be more appropriate for fitting my model. I ran it in R with glm.nb() and got results that matched just using lm on log transformed data pretty well. However, R indicated some warnings. I checked warnings(), and saw a list of warnings as follows: Warning messages: 1: non-integer x = 0.254825 I got the same error when trying to use the poisson family. My data is indeed continuous, not discrete (lots of non-integers). Does this mean that the model was not fit properly? Was data dropped when fitting the model? Is there an option to deal with this that I have overlooked? It would seem all is in order, but i just wanted to make sure. Thanks! Thanks. -Jarrett __ 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] Advice on visual graph packages
Hey, all. I'm looking for packages that are good at two things 1) Drawing directed graphs (i.e nodes and edges), both with single and double headed arrows, as well as allowing for differences in line width and solid versus dashed. Note: I've tried Rgraphviz here, but have run into some problems (which seem fixable and I may go with it in the end), and it doesn't satisfy need # 2 (which would be ideal if there is a package that does both). 2) Allowing a user to create a directed graph, and have some text object created that can be reprocessed easily reprocessed into a matrix representation, or other representation of my choosing. I've tried dynamicGraph, but it seems buggy, and continually either crashes, behaves very erratically (nodes disappearing when I modify edges), nor is it clear from the UI how one outputs a new graph, nor how one even accesses many graph attributes. This may be my own ignorance on the latter. Do you have any suggestions? Thanks! -Jarrett -- A Quick and (Very) Dirty Guide to Stats in R http://didemnid.ucdavis.edu/rtutorial.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 and provide commented, minimal, self-contained, reproducible code.
[R] comments in scan
I had a question about scan in R. For better code readability, I would like to have lines in the block of data to be scanned that are commented - not just lines that have a comment at the end. For example #age, weight, height 33,128,65 34,56,155 instead of having to do something like 33,128,65 #age, weight, height 34,56,155 Is this at all possible? __ 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] comments in scan
Right, but what if it is not the first line? For example #data block 1 34,54,23 23,53,12 35,23,23 #data block 2 64,24,13 354,24,12 324,13,3 On Nov 27, 2006, at 4:39 PM, Duncan Murdoch wrote: On 11/27/2006 7:25 PM, Jarrett Byrnes wrote: I had a question about scan in R. For better code readability, I would like to have lines in the block of data to be scanned that are commented - not just lines that have a comment at the end. For example #age, weight, height 33,128,65 34,56,155 instead of having to do something like 33,128,65 #age, weight, height 34,56,155 Is this at all possible? If it's on the first line, it's easy: just use skip=1 to skip the first line. The general case #age, weight, height 33,128,65 # and now for something completely different 34,56,155 probably needs something like this: scan(foo, what=list(0,0,0), comment.char=#, sep=,, multi.line=T) i.e. you need to tell it how many objects are in a record, and allow records to span multiple lines. Watch out for typos in the file that put different numbers of entries on each line, because scan() won't complain about it. 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 and provide commented, minimal, self-contained, reproducible code.
[R] Updating lmer - object is not subsettable?
I'm attempting to write a general function to implement Faraway's bootstrapping algorithm for mixed models with lmer, but have run into a curious problem. I'm comparing two models model.1-lmer(Response ~ Treatment + (1|Trial), data=exp.data, method=ML) model.2-lmer(Response ~ 1 + (1|Trial), data=exp.data, method=ML) When I attempt to update model.2 with simulated data, however, I get the following error: sim.data-unlist(simulate(model.1)) sim.model.2-update(model.2, sim.data~.) Error in x[[3]] : object is not subsettable Now, the following sim.model.1-update(model.1, sim.data~.) appears to work just fine. Does anyone know why update won't work, and is there something I can do about this? -Jarrett __ 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] Updating lmer - object is not subsettable?
Of course. The traceback is as follows. If you wish, I can
privately email you the data, as well as the function I'm working on.
17: all.names(x)
16: inherits(x, factor)
15: is.factor(table)
14: match(x, table, nomatch = 0)
13: / %in% all.names(x)
12: slashTerms(x[[3]])
11: FUN(X[[1]], ...)
10: lapply(bb, function(x) {
if (is.list(trms - slashTerms(x[[3]])))
return(lapply(unlist(makeInteraction(trms)), function
(trm) substitute(foo |
bar, list(foo = x[[2]], bar = trm
x
})
9: unlist(lapply(bb, function(x) {
if (is.list(trms - slashTerms(x[[3]])))
return(lapply(unlist(makeInteraction(trms)), function
(trm) substitute(foo |
bar, list(foo = x[[2]], bar = trm
x
}))
8: expandSlash(findbars(formula[[3]]))
7: lmer(formula = sim.data ~ (1 | Trial), data = exp.data, method =
ML)
6: lmer(formula = sim.data ~ (1 | Trial), data = exp.data, method =
ML)
5: eval(expr, envir, enclos)
4: eval(call, parent.frame())
3: .local(object, ...)
2: update(model.2, sim.data ~ .)
1: update(model.2, sim.data ~ .)
On Sep 12, 2006, at 8:16 PM, Douglas Bates wrote:
Can you show a traceback on this example? It may be related to a
problem that I just fixed in the development version of the lme4
package.
Alternatively if you can make the data available I can generate a
traceback myself.
On 9/12/06, Jarrett Byrnes [EMAIL PROTECTED] wrote:
I'm attempting to write a general function to implement Faraway's
bootstrapping algorithm for mixed models with lmer, but have run into
a curious problem. I'm comparing two models
model.1-lmer(Response ~ Treatment + (1|Trial), data=exp.data,
method=ML)
model.2-lmer(Response ~ 1 + (1|Trial), data=exp.data, method=ML)
When I attempt to update model.2 with simulated data, however, I get
the following error:
sim.data-unlist(simulate(model.1))
sim.model.2-update(model.2, sim.data~.)
Error in x[[3]] : object is not subsettable
Now, the following
sim.model.1-update(model.1, sim.data~.)
appears to work just fine. Does anyone know why update won't work,
and is there something I can do about this?
-Jarrett
__
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] Autoregressive Model with Independent Variable
On Mar 1, 2006, at 8:35 PM, Dirk Eddelbuettel wrote: On 1 March 2006 at 20:06, Jarrett Byrnes wrote: | Hey, all, I may just be missing something, but I'm trying to construct | a temporal autoregression with an independant variable other than just | what is happened at a previous point in time. So, the model structure | would be something like | | y(t)=b0+b1*y(t-1)+b2*y(t-2)...+a*x(t) | Yes: arima(), see in particular the xreg argument. Thanks so much! arima() seems to mostly fit the bill. I have data from multiple sites to use, as well. e.g. Timey1 x1 y2 x2 1 4 6 7 10 2 5 10 5 20 3 10 1 7 15 etc. I would like to use all of the sites in creating a model - I realize that the structure of the model would now be along the lines of: y(t)=b0+b1*y(t-1)+b2*y(t-2)...+a1*x(t)+a2*x(t-1)...+c Where c is the site effect - I know this can get all wrapped up in the intercept, but, how does one pass this data to arima() to make it work? I know that arima() takes a vector of y values - can it take a matrix of y values and a corresponding matrix of x values, or is there some other function that does this? -Jarrett __ 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] Autoregressive Model with Independent Variable
On Mar 1, 2006, at 8:35 PM, Dirk Eddelbuettel wrote: On 1 March 2006 at 20:06, Jarrett Byrnes wrote: | Hey, all, I may just be missing something, but I'm trying to construct | a temporal autoregression with an independant variable other than just | what is happened at a previous point in time. So, the model structure | would be something like | | y(t)=b0+b1*y(t-1)+b2*y(t-2)...+a*x(t) | Yes: arima(), see in particular the xreg argument. Thanks so much! arima() seems to mostly fit the bill. I have data from multiple sites to use, as well. e.g. Timey1 x1 y2 x2 1 4 6 7 10 2 5 10 5 20 3 10 1 7 15 etc. I would like to use all of the sites in creating a model - I realize that the structure of the model would now be along the lines of: y(t)=b0+b1*y(t-1)+b2*y(t-2)...+a1*x(t)+a2*x(t-1)...+c Where c is the site effect - I know this can get all wrapped up in the intercept, but, how does one pass this data to arima() to make it work? I know that arima() takes a vector of y values - can it take a matrix of y values and a corresponding matrix of x values, or is there some other function that does this? -Jarrett __ 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] Autoregressive Model with Independent Variable
Hey, all, I may just be missing something, but I'm trying to construct a temporal autoregression with an independant variable other than just what is happened at a previous point in time. So, the model structure would be something like y(t)=b0+b1*y(t-1)+b2*y(t-2)...+a*x(t) I'm even considering a model of y(t)=b0+b1*y(t-1)+b2*y(t-2)...+a1*x(t)+a2*x(t-1)... So, my data looks like Timey x 1 4 6 2 5 10 3 10 1 etc. When looking at ar() and similar methods, however, it seemed that the input was a single vector - say, in this case, the value y. Is there a method that allows me to specify an explicit model that would then incorporate x? __ 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] Canonical Values and Centroids for MANOVA plots
Hey, all, I'm trying to construct a centroid plot using canonical values from a MANOVA. I know that from the summary.manova object you can get Eigenvalues, and the H and E matrices (from SS$Treatment and SS$Residuals), but I am at a loss to get the loadings for the canonical values, nor values for the centroid centers and radii. Is there a package that does this that I am just missing, or any helpful code out there that I have been unable to locate? Thanks so much! -Jarrett p.s. There are a bunch of ecologists and evolutionary biologists here at Davis who, as a group, have just taken it upon themselves to learn R, and it's going swimmingly so far! Jarrett Byrnes Population Biology Graduate Group, UC Davis Bodega Marine Lab 707-875-1969 http://www-eve.ucdavis.edu/stachowicz/byrnes.shtml __ 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] Repeates Measures MANOVA for Time*Treatment Interactions
Dear R folk, First off I want to thank those of you who responded with comments for my R quick and dirty stats tutorial. They've been quite helpful, and I'm in the process of revising them. When it comes to repeated measures MANOVA, I'm in a bit of a bind, however. I'm beginning to see that all of the documentation is written for psychologists, who have a slightly different mind-set behind their experiments than, say, an ecologist, who is interested in the effects of time per se, and not just the effects of a treatment. For example, here's my dataset, say, looking at plant height in cm with and without fertilizer Treatment, Time1, Time2, Time3, Time4, Time5 Fertilizer, 1, 4, 8, 10, 12 Control,1,2,3,4,5 Fertilizer,1,8,10,12,20 Control,1,3,5,6,6 Fertilizer,2,5,10,20,25 Control,1,2,4,4,4 Clearly there is a time*treatment interaction (just eyeballing the dataset) My question is, how does one set this up using the anova.mlm approach so that in the end I can write up a table that says Treatment Time Time*Treatment I can see from ?anova.mlm how one would get the Treatment effect using something like response-with(my.data, rbind(Time1, Time2, Time3, Time4, Time5)) mlmfit-lm(response~1) mlmfit0-update(mlmfit, ~0) anova(mlmfit, mlmfit0, X = ~ Treatment, idata=my.data, test=Spherical) Although this yields the result that, after correction, it's not significant - perhaps due to the low DF from this simple example -- Analysis of Variance Table Model 1: response ~ 1 Model 2: response ~ 1 - 1 Contrasts orthogonal to ~Treatment Greenhouse-Geisser epsilon: 0.3565 Huynh-Feldt epsilon:0.4982 Res.Df Df Gen.var. F num Df den Df Pr(F) G-G Pr H-F Pr 1 4 0.43167 2 5 1 0.50937 3.7939 4 16 0.02356 0.09620 0.06966 -- But, I still want to get my time and time*treatment interactions - what would be the appropriate anova statements here? Thanks so much, and hopefully this will resolve the confusion both for myself and for LOTS of other ecology types! -Jarrett Jarrett Byrnes Population Biology Graduate Group, UC Davis Bodega Marine Lab 707-875-1969 http://www-eve.ucdavis.edu/stachowicz/byrnes.shtml __ 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] Repeates Measures MANOVA for Time*Treatment Interactions
Ah, so, if, say, I sampled on the 1st, 4th, 19th, 20th, and 30th day of a month, I should use: idata=data.frame(time=c(1,4,19,20,30)) On Nov 15, 2005, at 11:08 AM, Peter Dalgaard wrote: Peter Dalgaard [EMAIL PROTECTED] writes: No... idata is the *intra*-block structure, so it should just be your five times. I'm somewhat baffled that you're getting away with To clarify: Your five times, if anything. the default is a dataframe containing the single vector 1:p where p is the number of columns. This allows you to fit simple polynomials in time and test them for interaction with treatment, etc., but of course it requires that the times are equispaced, so if they are not, then you should use (say) idata=data.frame(time=c(0,3,6,9,12,18,24)) snip rest -- O__ Peter Dalgaard Øster Farimagsgade 5, Entr.B c/ /'_ --- Dept. of Biostatistics PO Box 2099, 1014 Cph. K (*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918 ~~ - ([EMAIL PROTECTED]) FAX: (+45) 35327907 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ 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] Type II and III sums of squares with Error in AOV
While my original design was balanced, I lost several replicates due to a storm, making the whole thing unbalanced. Ah, the realities of ecology. So, how does one look at individual strata, and then how would one report an aggregate test of the effect in general? On Nov 9, 2005, at 4:38 AM, Peter Dalgaard wrote: Prof Brian Ripley [EMAIL PROTECTED] writes: A multistratum aov() fit is just a list of aov() fits, so you can apply functions such as Anova to the individual strata. However, why do you want types II and III sums of squares? It is usual to do this type of analysis only with balanced designs. In the cases I can envisage that these make any sense, they are the same as type I (and in cases with only one treatment effect, they always are). I was about to make a similar comment. A possible exception is ANCOVA where you likely want to test both the within-stratum effect of a covariate and the effect of design factors adjusted for the covariate. On Tue, 8 Nov 2005, Jarrett Byrnes wrote: I've recently run into the problem of using aov with nested factors, and wanting to get the type II and III sums of squares. Normally Anova from the car package would do fine, but it doesn't like having an Error included, so my.aov -aov(Response ~ Treatment + Error(Treatment:Replicate)) Anova(my.aov, type=II) yields Error in Anova(nested.anova) : no applicable method for Anova And lm does not take Error as an argument. Nor does log()! Is that relevant? Error in eval(expr, envir, enclos) : couldn't find function Error Is there a way to get these additional types of sums of squares when Error is included in an aov model? Have you read the reference on the aov help page? That might be a good step to deepening your understanding of what Error() does. -- 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 -- O__ Peter Dalgaard Øster Farimagsgade 5, Entr.B c/ /'_ --- Dept. of Biostatistics PO Box 2099, 1014 Cph. K (*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918 ~~ - ([EMAIL PROTECTED]) FAX: (+45) 35327907 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ 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] Simple Nesting question/Odd error message
That works, but I am still puzzled as to why lme(fixed = X.open ~ Dock, random=~Slip%in%Dock, data=my.data) Does not work. It yields the error Error in getGroups.data.frame(dataMix, groups) : Invalid formula for groups Whereas were I to treat slip as a fixed factor aov(X.open ~ Dock + slipfactor%in%Dock, data=my.data) works just fine. Is there something fundamental that I am missing in the syntax? On Nov 8, 2005, at 2:06 AM, Dieter Menne wrote: Jarrett Byrnes redbeard at arrr.net writes: I'm having syntactical issues, however. When I try dock.lme-lme(X.open ~ Dock, random=Slip|~Dock, data=my.data) I get Error in inherits(object, reStruct) : Object Slip not found Syntax is wrong, should be something like dock.lme-lme(X.open ~ Dock, random=~Slip|Dock, data=my.data) __ 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] Type II and III sums of squares with Error in AOV
I've recently run into the problem of using aov with nested factors, and wanting to get the type II and III sums of squares. Normally Anova from the car package would do fine, but it doesn't like having an Error included, so my.aov -aov(Response ~ Treatment + Error(Treatment:Replicate)) Anova(my.aov, type=II) yields Error in Anova(nested.anova) : no applicable method for Anova And lm does not take Error as an argument. Error in eval(expr, envir, enclos) : couldn't find function Error Is there a way to get these additional types of sums of squares when Error is included in an aov model? __ 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] A Quick and (Very) Dirty Intro to Stats in R
Greetings to all, First off, I want to thank you all for answering any nagging questions I've had over the past few days. I've been in the process of putting together A Quick and (Very) Dirty Intro to Doing Your Statistics in R (which I have posted to http://didemnid.ucdavis.edu/rtutorial.html ) in order to teach an R workshop for the graduate students in my department. This is a guide for your everyday stats crunchers who want to free themselves from the cycle of SAS updates, have more flexibility than JMP or Statview will allow, but are not hardcore programming/think-about-stats-allday types. These are people who get data from the natural world, and then find out what it's telling them. So, to that end, I've put the guide together, and would be very interested in any comments you all would have. Are there any statistical methods that you think I really should have included for this type of audience that I didn't (and if it's over my head, would you be interested in contributing)? Is there anything just blatantly wrong or is unclear to a casual reader? Most importantly, there are still a few holes that need to be filled - if they can be 1) A SIMPLE explanation for how to do mixed models using lme. I am quite unsatisfied with most of what I've seen on the net, and if it even comes close to going over my head, it really won't fly with most folk I know. I've done the best I can, but I know if falls short. 2) A method of looking at type II and III sums of squares for aov if there is a different error term included. 3) How does one plot canonical values and centroid groupings for a MANOVA? 4) How does one use manova to do repeated measures? I've got the univariate method down, but would like to use manova a la the repeated statement in SAS. 5) Better output for post-hocs, and a Ryan's Q implementation. Thanks in advance for any input, and I hope this can be a resource to a lot of people! Jarrett Byrnes Population Biology Graduate Group, UC Davis Bodega Marine Lab 707-875-1969 http://www-eve.ucdavis.edu/stachowicz/byrnes.shtml __ 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 Nesting question/Odd error message
I'm attempting to analyze some survey data comparing multiple docks. I surveyed all of the slips within each dock, but as slips are nested within docks, getting multiple samples per slip, and don't really represent any meaningful gradient, slip is a random effect. There are also an unequal number of slips at each dock. I'm having syntactical issues, however. When I try dock.lme-lme(X.open ~ Dock, random=Slip|~Dock, data=my.data) I get Error in inherits(object, reStruct) : Object Slip not found I'm uncertain as to what this means, as Slip most certainly is there (yes, I checked) - any pointers, or pointers about syntax for this as a general topic. And while I'm on it, is there a way to fit models with random effects using least squares instead of REML? Thanks! -Jarrett __ 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] Plotting Factorial GLMs
Thanks to all who replied (although next time, reply to the list!) The simplest answer is to replace A in the curve statement with the following: rep(A, times=length(x)) Now if I could just figure out the workaround to plotting multiple sets of points onto my graph in the first place (i.e. break up the x,y values by block, then color them differently, and plot them on top of each other as add=TRUE doesn't seem to work for the plot statement - any pointers?) On Nov 3, 2005, at 10:31 PM, Jarrett Byrnes wrote: Hello all, I'm attempting to plot the functions from a generalized linear model while iterating over multiple levels of a factor in the model. In other words, I have a data set Block, Treatment.Level, Response.Level So, the glm and code to plot should be logit.reg-glm(formula = Response.Level ~ Treatment.Level + Block, family=quasibinomial(link=logit))) plot( Response.Level ~ Treatment.Level) logit.reg.function - function (trt, blk) predict(logit.reg, data.frame(Treatment.Level=trt, Block=blk) curve(logit.reg.function(x, A), add=TRUE) But I get the error: Error in xy.coords(x, y) : 'x' and 'y' lengths differ Now, if I set Block=A in the function, and take blk out, as well as taking the A out of the curve statement, it plots just fine. What am I doing wrong, as this would be a nice, quick, and easy way to whip up multiple curves from a factorial dataset! __ 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] Plotting Factorial GLMs
Hello all, I'm attempting to plot the functions from a generalized linear model while iterating over multiple levels of a factor in the model. In other words, I have a data set Block, Treatment.Level, Response.Level So, the glm and code to plot should be logit.reg-glm(formula = Response.Level ~ Treatment.Level + Block, family=quasibinomial(link=logit))) plot( Response.Level ~ Treatment.Level) logit.reg.function - function (trt, blk) predict(logit.reg, data.frame(Treatment.Level=trt, Block=blk) curve(logit.reg.function(x, A), add=TRUE) But I get the error: Error in xy.coords(x, y) : 'x' and 'y' lengths differ Now, if I set Block=A in the function, and take blk out, as well as taking the A out of the curve statement, it plots just fine. What am I doing wrong, as this would be a nice, quick, and easy way to whip up multiple curves from a factorial dataset! __ 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 Graphs in Powerpoint
Hey, all. Quick question. I'm attempting to use some of the great graphs generated in R for an upcoming talk that I'm writing in Powerpoint. Copying and pasting (I'm using OSX) yields graphs that look great in Powerpoint - until I resize them. Then fonts, points, and lines all become quite pixelated and blurry. Even if I size the window properly first, and then copy and paste in the graph, when I then view the slideshow, the graphs come out pixelated and blurry. Is there any good solution to this, or is this some fundamental incompatibility that I can't get around? -Jarrett __ 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] lm type of sums of squares
I'm curious, I realize there are methods for Type II and III sums of squares, and yet, when I've been constructing models with lm, I've noticed that position of the term of the model has not mattered in terms of its p-value. Does lm use sequential Type I sums of squares, or something else? Thanks! -Jarrett __ 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] Post Hoc Groupings
Quick question, as I attempt to learn R. For post-hoc tests 1) Is there an easy function that will take, say the results of tukeyHSD and create a grouping table. e.g., if I have treatments 1, 2, and 3, with 1 and 2 being statistically the same and 3 being different from both Group Treatment A 1 A 2 B 3 2) I've been stumbling over the proper syntax for simple effects for a tukeyHSD test. Is it TukeyHSD(model.aov, Treatment1, Treatment2) or TukeyHSD(model, c(Treatment1, Treatment2)) or something else, as neither of those seem to really work. __ 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] Post Hoc Groupings
Indeed, the following works as well On Oct 26, 2005, at 5:23 PM, P Ehlers wrote: fm1 - aov(breaks ~ wool*tension, data = warpbreaks) TukeyHSD(fm1, c(wool,tension, wool:tension)) However, when working with my own dataset, I get the following errors. I have some inkling this may be due to a slightly unbalanced sample size, but am uncertain of this. rich.aov - aov(X.open ~ Dock*Slip, data=dock_2004_data) TukeyHSD(rich.aov, c(Dock, Slip)) Error in rep.int(n, length(means)) : unimplemented type 'NULL' in 'rep' In addition: Warning messages: 1: non-factors ignored: Slip in: replications(paste(~, xx), data = mf) 2: non-factors ignored: Dock, Slip in: replications(paste(~, xx), data = mf) I am pleased to know that these errors are not quite the __ 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] Post Hoc Groupings
Indeed, that does it. Odd. I guess as Slip was a number, it needed to be categorized. Interesting... Not if only I can figure out how to get the TukeyHSD grouping table On Oct 26, 2005, at 7:35 PM, Mark Lyman wrote: Looking at the errors your code produces, it looks like you need to make Dock and Slip factors. dock_2004_data$Dockf-factor(dock_2004_data$Dock) dock_2004_data$Slipf-factor(dock_2004_data$Slip) rich.aov - aov(X.open ~ Dockf*Slipf, data=dock_2004_data) TukeyHSD(rich.aov, c(Dockf, Slipf)) Jarrett Byrnes wrote: Indeed, the following works as well On Oct 26, 2005, at 5:23 PM, P Ehlers wrote: fm1 - aov(breaks ~ wool*tension, data = warpbreaks) TukeyHSD(fm1, c(wool,tension, wool:tension)) However, when working with my own dataset, I get the following errors. I have some inkling this may be due to a slightly unbalanced sample size, but am uncertain of this. rich.aov - aov(X.open ~ Dock*Slip, data=dock_2004_data) TukeyHSD(rich.aov, c(Dock, Slip)) Error in rep.int(n, length(means)) : unimplemented type 'NULL' in 'rep' In addition: Warning messages: 1: non-factors ignored: Slip in: replications(paste(~, xx), data = mf) 2: non-factors ignored: Dock, Slip in: replications(paste(~, xx), data = mf) I am pleased to know that these errors are not quite the __ 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] Ryan's Q Post-Hoc for ANOVA
I'm using lm to run an ANOVA, and would like to use Ryan's Q as my post-hoc (as recommended by Day and Quinn, 1989, Ecological Monographs). I can't seem to find any methods in the base stats package that implement this post-hoc. Is there a good package of post-hoc methods out there, or has someone written a method for Ryan's Q previously? Thanks! -Jarrett __ 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] Getting univariate information from a multivariate data set
A quick question that I've had only partial success in answering. I have a multivariate dataset, and would like to extract some simple univariate information from it grouped by treatments, etc. I am encountering two problems however Note: I am importing my data with my_data - read.csv(/path/data.csv) 1) Scoping of unstack If I attempt sorted_data - unstack(response, response ~ treatment, data=my_data) I get Error in unstack(response, response ~ treatment, data=my_data): Object response not found However, if I first use attatch(my_data) and drop the data statement in unstack, I'm fine. This is all well and good, but I often work with multiple data sets with similar column headings, and would rather not have to attach and detach over and over again. 2) getting the univariate data Once, I attach, I find that I cannot then easily just use mean, sd, or anysuch, as mean(sorted_data) yields argument is not numeric or logical: returning NA in: mean.default(a) Nor does something like mean(sorted_data[1]) work - I thought perhaps breaking it down by treatment might help. Thanks in advance for any help! -Jarrett Jarrett Byrnes Population Biology Graduate Group, UC Davis Bodega Marine Lab 707-875-1969 http://www-eve.ucdavis.edu/stachowicz/byrnes.shtml __ 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
