[R] Applying ME in spdep to phylogenetic autocorrelation
Hi all, I want to fit a model with the form: trait 1 ~ trait 2 + environmental variable while controlling for the phylogenetic relatedness in my study species. I would use phylogenetic least squares regression, but I want to apply an independent effects analysis to partition the effects of trait 2 and environment on trait 1. (Partial correlation doesn't work here since all 3 variables are correlated). It looks like the ME function in the spdep package does exactly what I would like to do for phylogenetic autocorrelation, but for spatial autocorrelation. Is there a way to apply ME to phylogenetically structured instead of spatially structured data? I can use the mat2listw function to convert the phylogenetic distance matrix to a listw object, but I'm not sure I'm representing the neighbor relationships correctly. Should each taxa be connected to every other taxa by a weight that reflects the phylogenetic distance, or should only the most closely related tips be counted as neighbors? Should the phylogenetic distance matrix itself go into the mat2listw function, or should the values be converted to 1/distance, to be analogous to spatial weights? How can I use the S0, S1, and S2 values from the summary() of the listw object to figure out whether I've represented the neighbors correctly? Thanks very much for your help! Best, Megan Bartlett [[alternative HTML version deleted]] __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Correlating multiple effect sizes within a study to study-level predictors: metafor package
Hi Michael, I think you're right, it would be a good idea for me to get a better grip on mixed effects modeling before I charge ahead with metafor. Thank you very much for all of your help with this! Much appreciated! Best, Megan On Thu, Jul 17, 2014 at 3:50 AM, Michael Dewey i...@aghmed.fsnet.co.uk wrote: At 17:49 16/07/2014, Megan Bartlett wrote: Hi Michael, Thank you! Just to clarify, in my question, I was thinking that in this regression each study should be treated as one point, instead of each species, so that each effect size x value has a unique climate y value. Is that what the random= list(~1|Species, ~1|Site) argument is doing? No. Since the climate variable is per study (I assume) you are assuming that it has the same effect on each species. If that is not true you need to add species as another moderator and then add the interaction between climate and species. The random parameter is saying that each site has its own intercept but you are only estimating its variance and each species also has its own intercept drawn from another distribution whose variance is being estimated. I think you probably need to get local statistical help now from someone who understands the science of what you are doing and the statistics of mixed effects models. I am a bit concerned that without that knowledge we on the list may end up giving you misleading advice, Thanks, Megan On Wed, Jul 16, 2014 at 1:53 AM, Michael Dewey i...@aghmed.fsnet.co.uk wrote: At 23:19 14/07/2014, Megan Bartlett wrote: Thanks very much, Wolfgang and Michael! I feel like I understand rma much more clearly. But just to make sure, is there any way to do this kind of analysis for a continuous predictor variable? Yes, just put it in as a moderator. I am not sure I fully understand the rest of your question but the answer may be that the weights are a property of the individual effect sizes For each site level, I have a value for a climate variable, and it would be great to see whether the average effect size for each site is correlated with that climate variable. But I'm not sure what variance would produce the appropriate weighting for each site-level average- would it be the variance in effect sizes across species within each site? Or does this analysis not really make any sense for effect sizes? Thanks again! Best, Megan On Mon, Jul 14, 2014 at 6:06 AM, Viechtbauer Wolfgang (STAT) wolfgang.viechtba...@maastrichtuniversity.nl wrote: Somehow that initial post slipped under the radar for me ... Yes, I would give the same suggestion as Michael. Besides random effects for 'site', I would also suggest to add random effects for each estimates (as in a regular random-effects model). So, if you have an 'id' variable that is unique to each observed d-value, you would use: random = list(~ 1 | site, ~ 1 | id) with the rma.mv() function. This is in essence the model given by equation (6) in: Nakagawa, S., Santos, E. S. A. (2012). Methodological issues and advances in biological meta-analysis. Evolutionary Ecology, 26(5), 1253-1274. (at the time of publication, this model could not be fitted with metafor, but it can now). Same model is described with a bit more detail in: Konstantopoulos, S. (2011). Fixed effects and variance components estimation in three-level meta-analysis. Research Synthesis Methods, 2(1), 61-76. Best, Wolfgang -- Wolfgang Viechtbauer, Ph.D., Statistician Department of Psychiatry and Psychology School for Mental Health and Neuroscience Faculty of Health, Medicine, and Life Sciences Maastricht University, P.O. Box 616 (VIJV1) 6200 MD Maastricht, The Netherlands +31 (43) 388-4170 | http://www.wvbauer.com -Original Message- From: r-help-boun...@r-project.org [mailto:r-help-bounces@r- project.org] On Behalf Of Michael Dewey Sent: Monday, July 14, 2014 14:42 To: Megan Bartlett; r-help@r-project.org Subject: Re: [R] Correlating multiple effect sizes within a study to study-level predictors: metafor package At 23:18 11/07/2014, Megan Bartlett wrote: Hi everyone, Since metafor doesn't have its own list, I hope this is the correct place for this posting- my apologies if there is a more appropriate list. metafor questions welcome here, Megan Wolfgang seems to be off-list so while we wait for the definitive answer here are some hints. I'm conducting a meta-analysis where I would like to determine the correlation between plasticity in leaf traits and climate. I'm calculating effect sizes as Hedge's d. My data is structured so that each study collected data from one forest site, so there is one set of climate variable values for that study, and there are one or more species in each
Re: [R] Correlating multiple effect sizes within a study to study-level predictors: metafor package
Hi Michael, Thank you! Just to clarify, in my question, I was thinking that in this regression each study should be treated as one point, instead of each species, so that each effect size x value has a unique climate y value. Is that what the random= list(~1|Species, ~1|Site) argument is doing? Thanks, Megan On Wed, Jul 16, 2014 at 1:53 AM, Michael Dewey i...@aghmed.fsnet.co.uk wrote: At 23:19 14/07/2014, Megan Bartlett wrote: Thanks very much, Wolfgang and Michael! I feel like I understand rma much more clearly. But just to make sure, is there any way to do this kind of analysis for a continuous predictor variable? Yes, just put it in as a moderator. I am not sure I fully understand the rest of your question but the answer may be that the weights are a property of the individual effect sizes For each site level, I have a value for a climate variable, and it would be great to see whether the average effect size for each site is correlated with that climate variable. But I'm not sure what variance would produce the appropriate weighting for each site-level average- would it be the variance in effect sizes across species within each site? Or does this analysis not really make any sense for effect sizes? Thanks again! Best, Megan On Mon, Jul 14, 2014 at 6:06 AM, Viechtbauer Wolfgang (STAT) wolfgang.viechtba...@maastrichtuniversity.nl wrote: Somehow that initial post slipped under the radar for me ... Yes, I would give the same suggestion as Michael. Besides random effects for 'site', I would also suggest to add random effects for each estimates (as in a regular random-effects model). So, if you have an 'id' variable that is unique to each observed d-value, you would use: random = list(~ 1 | site, ~ 1 | id) with the rma.mv() function. This is in essence the model given by equation (6) in: Nakagawa, S., Santos, E. S. A. (2012). Methodological issues and advances in biological meta-analysis. Evolutionary Ecology, 26(5), 1253-1274. (at the time of publication, this model could not be fitted with metafor, but it can now). Same model is described with a bit more detail in: Konstantopoulos, S. (2011). Fixed effects and variance components estimation in three-level meta-analysis. Research Synthesis Methods, 2(1), 61-76. Best, Wolfgang -- Wolfgang Viechtbauer, Ph.D., Statistician Department of Psychiatry and Psychology School for Mental Health and Neuroscience Faculty of Health, Medicine, and Life Sciences Maastricht University, P.O. Box 616 (VIJV1) 6200 MD Maastricht, The Netherlands +31 (43) 388-4170 | http://www.wvbauer.com -Original Message- From: r-help-boun...@r-project.org [mailto:r-help-bounces@r- project.org] On Behalf Of Michael Dewey Sent: Monday, July 14, 2014 14:42 To: Megan Bartlett; r-help@r-project.org Subject: Re: [R] Correlating multiple effect sizes within a study to study-level predictors: metafor package At 23:18 11/07/2014, Megan Bartlett wrote: Hi everyone, Since metafor doesn't have its own list, I hope this is the correct place for this posting- my apologies if there is a more appropriate list. metafor questions welcome here, Megan Wolfgang seems to be off-list so while we wait for the definitive answer here are some hints. I'm conducting a meta-analysis where I would like to determine the correlation between plasticity in leaf traits and climate. I'm calculating effect sizes as Hedge's d. My data is structured so that each study collected data from one forest site, so there is one set of climate variable values for that study, and there are one or more species in each study, so all the species in a study have the same values for the climate variables. I'm not sure how to account for this structure in modeling the relationship between plasticity and climate. I think you need rma.mv for your situation and you need to specify a random effect for site. Try going ?rma.mv and looking for the section entitled Specifying random effects You will need to set up your dataframe with one row per species and an indicator variable for site and then use random = ~ 1 | site Not tested obviously and Wolfgang may have other suggestions My first thought was to calculate mean effect size and variance across species for every study with multiple species and correlate that with the climate variable values for those study with the rma() function, but trying to do that returns an error message: rma(yi = EffectSize, vi = Var, data = sitestable, mod = Precip) returns: Error in wi * (yi - X %*% b)^2 : non-conformable arrays This leaves me with two questions: 1) Am I even accounting for the data structure correctly with this approach, and 2) am I fundamentally misunderstanding how to use metafor to do so? Thanks
Re: [R] Correlating multiple effect sizes within a study to study-level predictors: metafor package
Thanks very much, Wolfgang and Michael! I feel like I understand rma much more clearly. But just to make sure, is there any way to do this kind of analysis for a continuous predictor variable? For each site level, I have a value for a climate variable, and it would be great to see whether the average effect size for each site is correlated with that climate variable. But I'm not sure what variance would produce the appropriate weighting for each site-level average- would it be the variance in effect sizes across species within each site? Or does this analysis not really make any sense for effect sizes? Thanks again! Best, Megan On Mon, Jul 14, 2014 at 6:06 AM, Viechtbauer Wolfgang (STAT) wolfgang.viechtba...@maastrichtuniversity.nl wrote: Somehow that initial post slipped under the radar for me ... Yes, I would give the same suggestion as Michael. Besides random effects for 'site', I would also suggest to add random effects for each estimates (as in a regular random-effects model). So, if you have an 'id' variable that is unique to each observed d-value, you would use: random = list(~ 1 | site, ~ 1 | id) with the rma.mv() function. This is in essence the model given by equation (6) in: Nakagawa, S., Santos, E. S. A. (2012). Methodological issues and advances in biological meta-analysis. Evolutionary Ecology, 26(5), 1253-1274. (at the time of publication, this model could not be fitted with metafor, but it can now). Same model is described with a bit more detail in: Konstantopoulos, S. (2011). Fixed effects and variance components estimation in three-level meta-analysis. Research Synthesis Methods, 2(1), 61-76. Best, Wolfgang -- Wolfgang Viechtbauer, Ph.D., Statistician Department of Psychiatry and Psychology School for Mental Health and Neuroscience Faculty of Health, Medicine, and Life Sciences Maastricht University, P.O. Box 616 (VIJV1) 6200 MD Maastricht, The Netherlands +31 (43) 388-4170 | http://www.wvbauer.com -Original Message- From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On Behalf Of Michael Dewey Sent: Monday, July 14, 2014 14:42 To: Megan Bartlett; r-help@r-project.org Subject: Re: [R] Correlating multiple effect sizes within a study to study-level predictors: metafor package At 23:18 11/07/2014, Megan Bartlett wrote: Hi everyone, Since metafor doesn't have its own list, I hope this is the correct place for this posting- my apologies if there is a more appropriate list. metafor questions welcome here, Megan Wolfgang seems to be off-list so while we wait for the definitive answer here are some hints. I'm conducting a meta-analysis where I would like to determine the correlation between plasticity in leaf traits and climate. I'm calculating effect sizes as Hedge's d. My data is structured so that each study collected data from one forest site, so there is one set of climate variable values for that study, and there are one or more species in each study, so all the species in a study have the same values for the climate variables. I'm not sure how to account for this structure in modeling the relationship between plasticity and climate. I think you need rma.mv for your situation and you need to specify a random effect for site. Try going ?rma.mv and looking for the section entitled Specifying random effects You will need to set up your dataframe with one row per species and an indicator variable for site and then use random = ~ 1 | site Not tested obviously and Wolfgang may have other suggestions My first thought was to calculate mean effect size and variance across species for every study with multiple species and correlate thatwith the climate variable values for those study with the rma() function, but trying to do that returns an error message: rma(yi = EffectSize, vi = Var, data = sitestable, mod = Precip) returns: Error in wi * (yi - X %*% b)^2 : non-conformable arrays This leaves me with two questions: 1) Am I even accounting for the data structure correctly with this approach, and 2) am I fundamentally misunderstanding how to use metafor to do so? Thanks very much for your help! Best, Megan [[alternative HTML version deleted]] __ R-help@r-project.org 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] Correlating multiple effect sizes within a study to study-level predictors: metafor package
Hi everyone, Since metafor doesn't have its own list, I hope this is the correct place for this posting- my apologies if there is a more appropriate list. I'm conducting a meta-analysis where I would like to determine the correlation between plasticity in leaf traits and climate. I'm calculating effect sizes as Hedge's d. My data is structured so that each study collected data from one forest site, so there is one set of climate variable values for that study, and there are one or more species in each study, so all the species in a study have the same values for the climate variables. I'm not sure how to account for this structure in modeling the relationship between plasticity and climate. My first thought was to calculate mean effect size and variance across species for every study with multiple species and correlate thatwith the climate variable values for those study with the rma() function, but trying to do that returns an error message: rma(yi = EffectSize, vi = Var, data = sitestable, mod = Precip) returns: Error in wi * (yi - X %*% b)^2 : non-conformable arrays This leaves me with two questions: 1) Am I even accounting for the data structure correctly with this approach, and 2) am I fundamentally misunderstanding how to use metafor to do so? Thanks very much for your help! Best, Megan [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Getting a root edge error when trying to read the phylocom megatree into R
Hi, I'm trying to use the picante package in R to build phylogenetic trees, based on a list of taxa I have data for and the phylocom APG3 megatree (version R20091110; http://www.phylodiversity.net/phylomatic/). However, trying to read the tree in R using: read.tree(R20091110.new.txt) Gives the error message: Error in if (sum(obj[[i]]$edge[, 1] == ROOT) == 1 dim(obj[[i]]$edge)[1] : missing value where TRUE/FALSE needed This seems like such a basic question, but I don't know what the problem with the root edge is. Are there any modifications you have to make to the megatree in order to read it in R? Thanks so much for your help! Best, Megan Bartlett [[alternative HTML version deleted]] __ R-help@r-project.org 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] Adding points to a wireframe: 'x and units must have length 0' error
Hi, I'm trying to follow the suggestions given by Deepayan Sarkar in this message: http://tolstoy.newcastle.edu.au/R/help/05/11/16135.html to plot 3-D points on a wireframe plot. The problem is that I keep getting a partly formed plot- with the colored lattice visible but no axis labels or additional points- with the error message error using packet 1, 'x' and 'units' must have length 0. Does anyone know what I'm doing wrong? The code I've been using is: # generate some fake data data.frame(x = seq(-4, 0, 0.5), y = seq(0, 40, 5))- df expand.grid(x = df$x, y = df$y) - gridd (gridd$y* gridd$x) - gridd$z data.frame(x = runif(10, -4, 0), y = runif(10, 0, 40))- pts1 pts1$z - pts1$x*pts1$y # plot wireframe(z ~ y*x, gridd, drape = TRUE, colorkey= TRUE, scales = list(arrows = FALSE), pts = pts1, panel.3d.wireframe = function(x, y, z,xlim, ylim, zlim, xlim.scaled, ylim.scaled, zlim.scaled, pts1, ...) { panel.3dwire(x = x, y = y, z = z, xlim = c(-4, 0), ylim = c(0, 40), zlim = c(-160, 0), xlim.scaled = c(-0.5, 0.5), ylim.scaled = c(-0.5, 0.5), zlim.scaled = c(-0.5, 0.5), ...) xx - xlim.scaled[1] + diff(xlim.scaled) * (pts1$x - xlim[1]) / diff(xlim) yy - ylim.scaled[1] + diff(ylim.scaled) * (pts1$y - ylim[1]) / diff(ylim) zz - zlim.scaled[1] + diff(zlim.scaled) * (pts1$z - zlim[1]) / diff(zlim) panel.3dscatter(x = xx, y = yy, z = zz, xlim = xlim, ylim = ylim, zlim = zlim, xlim.scaled = xlim.scaled, ylim.scaled = ylim.scaled, zlim.scaled = zlim.scaled, ...) }) Thanks so much! -Megan Bartlett [[alternative HTML version deleted]] __ R-help@r-project.org 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 points to a wireframe: 'x and units must have length 0' error
Hi David, Thanks for the suggestion - I changed pts1 to pts, but I still got the same error as before. Do you know what else I'm doing wrong? Thanks, Megan PS. New code is: wireframe(z ~ y*x, gridd, drape = TRUE, colorkey= TRUE, scales = list(arrows = FALSE), pts = pts1, panel.3d.wireframe = function(x, y, z,xlim, ylim, zlim, xlim.scaled, ylim.scaled, zlim.scaled, pts, ...) { panel.3dwire(x = x, y = y, z = z, xlim = c(-4, 0), ylim = c(0, 40), zlim = c(-160, 0), xlim.scaled = c(-0.5, 0.5), ylim.scaled = c(-0.5, 0.5), zlim.scaled = c(-0.5, 0.5), ...) xx - xlim.scaled[1] + diff(xlim.scaled) * (pts$x - xlim[1]) / diff(xlim) yy - ylim.scaled[1] + diff(ylim.scaled) * (pts$y - ylim[1]) / diff(ylim) zz - zlim.scaled[1] + diff(zlim.scaled) * (pts$z - zlim[1]) / diff(zlim) panel.3dscatter(x = xx, y = yy, z = zz, xlim = xlim, ylim = ylim, zlim = zlim, xlim.scaled = xlim.scaled, ylim.scaled = ylim.scaled, zlim.scaled = zlim.scaled, ...) }) On Mon, Oct 24, 2011 at 12:04 PM, David Winsemius dwinsem...@comcast.netwrote: On Oct 24, 2011, at 1:23 PM, Megan Bartlett wrote: Hi, I'm trying to follow the suggestions given by Deepayan Sarkar in this message: http://tolstoy.newcastle.edu.**au/R/help/05/11/16135.htmlhttp://tolstoy.newcastle.edu.au/R/help/05/11/16135.html to plot 3-D points on a wireframe plot. The problem is that I keep getting a partly formed plot- with the colored lattice visible but no axis labels or additional points- with the error message error using packet 1, 'x' and 'units' must have length 0. Does anyone know what I'm doing wrong? The code I've been using is: # generate some fake data data.frame(x = seq(-4, 0, 0.5), y = seq(0, 40, 5))- df expand.grid(x = df$x, y = df$y) - gridd (gridd$y* gridd$x) - gridd$z data.frame(x = runif(10, -4, 0), y = runif(10, 0, 40))- pts1 pts1$z - pts1$x*pts1$y # plot wireframe(z ~ y*x, gridd, drape = TRUE, colorkey= TRUE, scales = list(arrows = FALSE), pts = pts1, panel.3d.wireframe = function(x, y, z,xlim, ylim, zlim, xlim.scaled, ylim.scaled, zlim.scaled, pts1, ...) { panel.3dwire(x = x, y = y, z = z, xlim = c(-4, 0), ylim = c(0, 40), zlim = c(-160, 0), xlim.scaled = c(-0.5, 0.5), ylim.scaled = c(-0.5, 0.5), zlim.scaled = c(-0.5, 0.5), ...) xx - xlim.scaled[1] + diff(xlim.scaled) * (pts1$x - xlim[1]) / diff(xlim) yy - ylim.scaled[1] + diff(ylim.scaled) * (pts1$y - ylim[1]) / diff(ylim) zz - zlim.scaled[1] + diff(zlim.scaled) * (pts1$z - zlim[1]) / diff(zlim) panel.3dscatter(x = xx, You assigned 'pts1' to 'pts', so you should have used 'pts' above: y = yy, z = zz, xlim = xlim, ylim = ylim, zlim = zlim, xlim.scaled = xlim.scaled, ylim.scaled = ylim.scaled, zlim.scaled = zlim.scaled, ...) }) Thanks so much! -Megan Bartlett [[alternative HTML version deleted]] __** R-help@r-project.org mailing list https://stat.ethz.ch/mailman/**listinfo/r-helphttps://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/** posting-guide.html http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. David Winsemius, MD West Hartford, CT [[alternative HTML version deleted]] __ R-help@r-project.org 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 points to a wireframe: 'x and units must have length 0' error
Hi, I got the same result with David's version, but defining the *lim values prior to the function actually made it work. I also renamed the points pts so it wouldn't matter whether it was pts1 or pts. The new code is: data.frame(x = seq(-4, 0, 0.5), y = seq(0, 40, 5))- df expand.grid(x = df$x, y = df$y) - gridd (gridd$y* gridd$x) - gridd$z data.frame(x = runif(10, -4, 0), y = runif(10, 0, 40))- pts pts$z - pts$x*pts$y wireframe(z ~ x*y, gridd, drape = TRUE, colorkey= TRUE, scales = list(arrows = FALSE), xlim = c(-4, 0), ylim = c(0, 40), zlim = c(-160, 0), pts = pts, panel.3d.wireframe = function(x, y, z,xlim, ylim, zlim, xlim.scaled, ylim.scaled, zlim.scaled, pts, ...) { panel.3dwire(x = x, y = y, z = z, xlim.scaled = c(-0.5, 0.5), ylim.scaled = c(-0.5, 0.5), zlim.scaled = c(-0.5, 0.5), ...) xx - xlim.scaled[1] + diff(xlim.scaled) * (pts$x - xlim[1]) / diff(xlim) yy - ylim.scaled[1] + diff(ylim.scaled) * (pts$y - ylim[1]) / diff(ylim) zz - zlim.scaled[1] + diff(zlim.scaled) * (pts$z - zlim[1]) / diff(zlim) panel.3dscatter(x = xx, y = yy, z = zz, xlim = xlim, ylim = ylim, zlim = zlim, xlim.scaled = xlim.scaled, ylim.scaled = ylim.scaled, zlim.scaled = zlim.scaled, ...) }) Thanks so much!! -Megan On Mon, Oct 24, 2011 at 2:25 PM, Dennis Murphy djmu...@gmail.com wrote: Hi David: When I try your code, I get the wireframe with the x, y, z axes sans bounding cube and points, along with the error message Error using packet 1 object 'pts' not found sessionInfo() R version 2.13.1 (2011-07-08) Platform: x86_64-pc-mingw32/x64 (64-bit) locale: [1] LC_COLLATE=English_United States.1252 [2] LC_CTYPE=English_United States.1252 [3] LC_MONETARY=English_United States.1252 [4] LC_NUMERIC=C [5] LC_TIME=English_United States.1252 attached base packages: [1] stats graphics grDevices utils datasets methods base other attached packages: [1] sos_1.3-1 brew_1.0-6 lattice_0.19-33 loaded via a namespace (and not attached): [1] grid_2.13.1 tools_2.13.1 Dennis On Mon, Oct 24, 2011 at 1:45 PM, David Winsemius dwinsem...@comcast.net wrote: On Oct 24, 2011, at 4:12 PM, Megan Bartlett wrote: Hi David, Thanks for the suggestion - I changed pts1 to pts, but I still got the same error as before. Do you know what else I'm doing wrong? I do not. The code that works for me is: wireframe(z ~ y*x, gridd, drape = TRUE, colorkey= TRUE, scales = list(arrows = FALSE), pts = pts1, panel.3d.wireframe = function(x, y, z,xlim, ylim, zlim, xlim.scaled, ylim.scaled, zlim.scaled, pts1, ...) { panel.3dwire(x = x, y = y, z = z, xlim = c(-4, 0), ylim = c(0, 40), zlim = c(-160, 0), xlim.scaled = c(-0.5, 0.5), ylim.scaled = c(-0.5, 0.5), zlim.scaled = c(-0.5, 0.5), ...) xx - xlim.scaled[1] + diff(xlim.scaled) * (pts$x - xlim[1]) / diff(xlim) yy - ylim.scaled[1] + diff(ylim.scaled) * (pts$y - ylim[1]) / diff(ylim) zz - zlim.scaled[1] + diff(zlim.scaled) * (pts$z - zlim[1]) / diff(zlim) panel.3dscatter(x = xx, y = yy, z = zz, xlim = xlim, ylim = ylim, zlim = zlim, xlim.scaled = xlim.scaled, ylim.scaled = ylim.scaled, zlim.scaled = zlim.scaled, ...) }) -- David. Thanks, Megan PS. New code is: wireframe(z ~ y*x, gridd, drape = TRUE, colorkey= TRUE, scales = list(arrows = FALSE), pts = pts1, panel.3d.wireframe = function(x, y, z,xlim, ylim, zlim, xlim.scaled, ylim.scaled, zlim.scaled, pts, ...) { panel.3dwire(x = x, y = y, z = z, xlim = c(-4, 0), ylim = c(0, 40), zlim = c(-160, 0), xlim.scaled = c(-0.5, 0.5), ylim.scaled = c(-0.5, 0.5), zlim.scaled = c(-0.5, 0.5), ...) xx - xlim.scaled[1] + diff(xlim.scaled) * (pts$x - xlim[1]) / diff(xlim) yy - ylim.scaled[1] + diff(ylim.scaled) * (pts$y - ylim[1]) / diff(ylim) zz - zlim.scaled[1] + diff(zlim.scaled) * (pts$z - zlim[1]) / diff(zlim) panel.3dscatter(x = xx, y = yy, z = zz, xlim = xlim, ylim = ylim