[R] incompatible dimensions error
Hi Listers Using mvpart to run a MV regression tree with PCA= TRUE to get a PCA plotted with sites coloured according to the tree output. Unfortunately it wont produce the PCA, instead giving the error message.. Error in cor(xall, xx[order(tree$where), ]) : incompatible dimensions. However, when I run a PCA on the data using the rda command I have no problems producing a PCA. data is attached as a text file my code is thus... fish05.hel <- decostand(fish05,"hellinger") fish05.mrt <- mvpart(data.matrix(fish05.hel)~.,env,margin=0.08,cp=0,rsq=TRUE,xv="pick",xval=nrow(fish05),xvmult=100,which=4,pca=TRUE) The tree is produced no problem but it wont produce a PCA. I am just keen to understand what this error means as I dont see anything unusual about the dataset used, notwithstanding the data is rather messy. Andy -- Andrew Halford Ph.D Senior Coastal Fisheries Scientist Pacific Community | Communauté du Pacifique CPS – B.P. D5 | 98848 Noumea, New Caledonia | Nouméa, Nouvelle-Calédonie CTE.STRICHL.SORDCHA.LTUSLAB.DIMI ACA.NANSMON.GRANZAN.CORNGOM.VARIHAL.HORT THA.LUTECHA.AURICEN.FLAVZEB.SCOPBAL.UNDU CHA.RETINAS.LITUCHA.EPHICHE.FASCTHA.HARD HEM.MELACEP.URODCHL.BLEESCA.SCHLMEL.VIDU CHE.CHLOEPB.INSISCA.PSITZEB.VELICHA.TLIS CHA.VAGALAB.BICOHIP.LONGEPI.MERRACA.NUDA CHA.LUNUCHA.ULIETHA.QUINCEP.ARGUCHL.MICR SCA.NIGEHAL.TRIMLUT.FULVLUT.GIBBSCA.DIMI CHA.CITRMAC.MELEGND.AUREAPH.FURCLUT.MONO PYG.DIACACA.TRIORHI.ACULCHE.TRILTHA.LUNA ACA.NCUSHEM.FASCTHA.PURPSCA.GLOBSCA.SPP SIG.ARGEACA.NRISHEN.CHRYHAL.CHRYOXY.UNIF SCA.FORSSCA.OVICACA.PYRONAS.UNICMON.HETE CET.OCELSCA.FRENACA.THOMFOR.FLAVPSC.HEXA STJ.BANDSCA.FLAVPSE.PASCSIG.VULPACA.LINE ACA.OLIVMEL.NIGESUF.BURSCHA.KLEICHE.UNDU COR.GAIMHAL.RICHMCR.NIGESCA.GHOBCHA.MERT FOR.LONGHAL.MELALUT.BOHACEN.LORICTE.CYAN SUF.CHRYHTY.POLYLBY.UNILPSC.TETRTHA.AMBL ACA.BLOCCHA.MELACHA.MEYECHE.SPP HAL.MTUS NOV.TAENLET.HARALUT.SEMISCO.LINECEN.BICO SCA.SPINACA.MACUACA.SPP CHA.LINECHA.PUNC KYP.VAIGANA.MELECOR.BATULAB.PECTOXY.DIGR STJ.STRIMCR.MACUSCA.ALTISCA.RIVUPMS.AREO CTE.STRGNAS.BREVNAS.HEXACHA.SEMEHEN.MONO ANA.MELAANA.TWISCOR.AYGUHAL.SPP SCO.BILI CEN.VROLPCA.IMPESCA.FESTSIG.FUSCACA.MATA NAS.CAESHAL.PROSHAL.SCAPLBS.MICRTHA.TRIL SCA.RUBRPSE.BARTSIG.PUEL DrenmeoMPA_Tr4_2011 75 41 5.5 38.50 3 0 1 0 2 3.5 9.5 26 0 0.5 1 0.5 2 7.5 1.5 1 0 0 0 0 0 0 0.5 3.5 0 1 2.5 1.5 0 0.5 1.5 0 0.5 0 0.5 3 0 0 0 0 0 0 0 1.5 0.5 0 0 0 0 0 0 0 21.32.5 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0.5 0 0 0 0 0 0 0 0.5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 0 0 0 0 0 0 0 0 4.5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 14 0 0 0 0 0 0 0 0 0 0 0 0 0 0 DrenmeoMPA_Tr10_20114.5 9 2 12.50 0 0 1 1 2.5 0 0.5 0 0 1 0 0 0 0.5 0 0 0 0 0 1.5 0 0 0 2 0 0 0 1 0 0 0 0 0 0 0 28 0 0 0 0 0
[R] Fwd: defining group colours in a call to rda
-- Forwarded message -
From: Abby Spurdle
Date: Wed, Aug 5, 2020 at 3:07 PM
Subject: Re: [R] defining group colours in a call to rda
To: Andrew Halford
Hi Andrew,
Perhaps you want this:
cols <- rep_len (c ("red", "green", "blue", "aquamarine", "magenta"), 9)
cols
Or this:
cols = rep ("", 9)
cols [unique (MI_fish_all.mrt$where)] = c ("red", "green", "blue",
"aquamarine", "magenta")
cols
Then you can substitute either into your original example:
plotcolor <- cols [MI_fish_all.mrt$where]
On Wed, Aug 5, 2020 at 1:02 PM Andrew Halford
wrote:
>
> Hi Abby,
>
> Apologies for not providing more info but you have worked out what I was
on about anyways.
>
> I thought it would scroll through and allocate the colours to each unique
number sequentially. I will add more colours to my vector but I would like
to know if it is possible to do what I originally hoped for.
>
> cheers
>
> Andy
>
> On Wed, Aug 5, 2020 at 8:40 AM Abby Spurdle wrote:
>>
>> Hi,
>>
>> Your example is not reproducible.
>> However, I suspect that the following is the problem:
>>
>> c("red","green","blue","aquamarine","magenta")[MI_fish_all.mrt$where]
>>
>> Here's my version:
>>
>> where = c (3, 3, 8, 6, 6, 9, 5, 5, 9, 3, 8, 6, 9, 6, 5, 9, 5, 3, 8, 6,
>> 9, 6, 5, 9, 5, 3, 3, 8, 6, 6, 9, 5, 5, 9, 6, 9, 5, 9)
>> unique (where)
>>
>> c("red", "green", "blue", "aquamarine", "magenta")[where]
>>
>> There's five colors.
>> But only two of the indices are within one to five.
>> So, the resulting color vector contains missing values.
>>
>> In the base graphics system, if you set colors to NA, it usually means
no color.
>>
>> I'm not sure exactly what you want to do, but I'm assuming you can fix
>> it from here.
>>
>> On Tue, Aug 4, 2020 at 9:49 PM Andrew Halford
wrote:
>> >
>> > Hi,
>> >
>> > I've been trying to use the output on group membership of the final
leaves
>> > in a MRT analysis to define my own colours, however I am not getting
the
>> > result I'm after.
>> >
>> > Here is the code
>> > fish.pca <-rda(fish_all.hel,scale=TRUE)
>> > fish.site <- scores(fish.pca,display="sites",scaling=3)
>> > fish.spp <-
>> >
scores(fish.pca,display="species",scaling=3)[fish.MRT.indval$pval<=0.05,]
>> > plot(fish.pca,display=c("sites","species"),type="n",scaling=3)
>> > points(fish.site,pch=21,bg=MI_fish_all.mrt$where,cex=1.2)
>> > plotcolor <-
>> > c("red","green","blue","aquamarine","magenta")[MI_fish_all.mrt$where]
>> > fish.pca <-rda(fish_all.hel,scale=TRUE)
>> > plot(fish.pca,display=c("sites","species"),type="n",scaling=3)
>> > points(fish.site,pch=21,bg=plotcolor,cex=1.2)
>> > MI_fish_all.mrt$where
>> >
>> > If I run the points command and insert the group membership direct
from the
>> > MRT analysis e.g. bg=MI_fish_all.mrt$where , then the subsequent
points
>> > plot up correctly with a different colour for each group.However if I
try
>> > to impose my own colour combo with plotcolor.It prints colours for
2
>> > groups and leaves the rest uncoloured.
>> >
>> > The call to MI_fish_all.mrt$where gives...
>> > [1] 3 3 8 6 6 9 5 5 9 3 8 6 9 6 5 9 5 3 8 6 9 6 5 9 5 3 3 8 6 6 9 5 5
9 6
>> > 9 5 9.
>> >
>> > These are the split groupings for all 39 sites in the analysis and
there
>> > are 5 numbers corresponding to 5 final leaves in the tree.
>> >
>> > I cant see why my colour scheme isnt being recognised.
>> >
>> > All help accepted.
>> >
>> > Andy
>> >
>> >
>> > --
>> > Andrew Halford Ph.D
>> > Senior Coastal Fisheries Scientist
>> > Pacific Community | Communauté du Pacifique CPS – B.P. D5 | 98848
Noumea,
>> > New Caledonia | Nouméa, Nouvelle-Calédonie
>> >
>> > [[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.
>
>
>
> --
> Andrew Halford Ph.D
> Senior Coastal Fisheries Scientist
> Pacific Community | Communauté du Pacifique CPS – B.P. D5 | 98848 Noumea,
> New Caledonia | Nouméa, Nouvelle-Calédonie
--
Andrew Halford Ph.D
Senior Coastal Fisheries Scientist
Pacific Community | Communauté du Pacifique CPS – B.P. D5 | 98848 Noumea,
New Caledonia | Nouméa, Nouvelle-Calédonie
[[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.
[R] defining group colours in a call to rda
Hi,
I've been trying to use the output on group membership of the final leaves
in a MRT analysis to define my own colours, however I am not getting the
result I'm after.
Here is the code
fish.pca <-rda(fish_all.hel,scale=TRUE)
fish.site <- scores(fish.pca,display="sites",scaling=3)
fish.spp <-
scores(fish.pca,display="species",scaling=3)[fish.MRT.indval$pval<=0.05,]
plot(fish.pca,display=c("sites","species"),type="n",scaling=3)
points(fish.site,pch=21,bg=MI_fish_all.mrt$where,cex=1.2)
plotcolor <-
c("red","green","blue","aquamarine","magenta")[MI_fish_all.mrt$where]
fish.pca <-rda(fish_all.hel,scale=TRUE)
plot(fish.pca,display=c("sites","species"),type="n",scaling=3)
points(fish.site,pch=21,bg=plotcolor,cex=1.2)
MI_fish_all.mrt$where
If I run the points command and insert the group membership direct from the
MRT analysis e.g. bg=MI_fish_all.mrt$where , then the subsequent points
plot up correctly with a different colour for each group.However if I try
to impose my own colour combo with plotcolor.It prints colours for 2
groups and leaves the rest uncoloured.
The call to MI_fish_all.mrt$where gives...
[1] 3 3 8 6 6 9 5 5 9 3 8 6 9 6 5 9 5 3 8 6 9 6 5 9 5 3 3 8 6 6 9 5 5 9 6
9 5 9.
These are the split groupings for all 39 sites in the analysis and there
are 5 numbers corresponding to 5 final leaves in the tree.
I cant see why my colour scheme isnt being recognised.
All help accepted.
Andy
--
Andrew Halford Ph.D
Senior Coastal Fisheries Scientist
Pacific Community | Communauté du Pacifique CPS – B.P. D5 | 98848 Noumea,
New Caledonia | Nouméa, Nouvelle-Calédonie
[[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] multiple graphs on one plot
Hi Jim,
Many thanks for your help and yes CW is carapace width. Here is the final
coding I used...I set the peak of the curves at max frequency bin for each
sex. I also added the means and SD's from my data. According to my visual
diagnostics (qqplots, density plots) the frequency distributions do appear
normally distributed in this case. To be honest Ive struggled to find
whether CW is normally distributed in mud crab populations or not.
f<- lf_crabs$cw[lf_crabs$sex=='female']
m<- lf_crabs$cw[lf_crabs$sex=='male']
# find mean and sd to determine normal curve dimensions
m_m <-mean(m)
sd_m <-sqrt(var(m))
m_f <-mean(f)
sd_f <-sqrt(varf(f))
mf <- list(f,m)
multhist(mf, xlab="CW", ylab="Frequency", ylim=c(0,100),main="All Measured
Crabs", col=c("dark gray", "light gray"),
breaks=seq(90,210, by=10),beside=TRUE,space=c(0,0.5))
legend("topright", c("Females", "Males"), fill=c("dark gray", "light gray"))
lines(seq(0,32,length.out=121),rescale(dnorm(90:210,145.4867,20.99906),c(0,50)),col="dark
gray",lwd=2)
lines(seq(0,32,length.out=121),rescale(dnorm(90:210,151.0783,21.88299),c(0,80)),col="light
gray",lwd=2)
abline(v=145.4867,lwd=2,col="red")
On Mon, 13 May 2019 at 17:30, Jim Lemon wrote:
> Hi Andrew,
> First, a little mind reading. My crystal ball says that "cw" can be
> interpreted as "carapace width". It didn't tell me the parameters of
> the distribution, so:
>
> set.seed(1234)
> mf<-list(rnorm(400,145,15),rnorm(400,160,15))
> library(plotrix)
> multhist(mf, xlab="CW", ylab="Frequency", ylim=c(0,100),main="All Measured
> Crabs", col=c("dark gray", "light gray"),
> breaks=seq(90,210, by=10),beside=TRUE,space=c(0,0.5))
> legend("topright", c("Females", "Males"), fill=c("dark gray", "light
> gray"))
> lines(seq(0,32,length.out=121),rescale(dnorm(90:210,145,15),c(0,100)))
>
> This produces what I think you are after. Note that it may be
> misleading as the distribution of carapace width in real mud crabs
> doesn't look normal to me.
>
> Jim
>
> On Mon, May 13, 2019 at 3:00 PM Andrew Halford
> wrote:
> >
> > Hi Listers
> >
> > I've been trying to make a single graphic that has frequency histograms
> for
> > male and female mud crabs displayed side by side (such as when using the
> > beside=TRUE command for barplots). I then want to display a normal
> > distribution on top of the male and female histograms.
> >
> > I have been using the multhist command in Plotrix to generate the
> > histograms without too much problem, but I cannot get the normal
> > distributions to plot up on the same graph.
> >
> > Histograms plot
> >
> > mf <-
> >
> list(lf_crabs$cw[lf_crabs$sex=='female'],lf_crabs$cw[lf_crabs$sex=='male'])
> > multhist(mf, xlab="CW", ylab="Frequency", ylim=c(0,100),main="All
> Measured
> > Crabs", col=c("dark gray", "light gray"),
> > breaks=seq(90,210, by=10),beside=TRUE,space=c(0,0.5))
> > legend("topright", c("Females", "Males"), fill=c("dark gray", "light
> gray"))
> >
> > Then I try to add a normal distribution curve just to the female data
> but I
> > cant get the output to plot
> >
> > points(seq(min(mf[[1]]), max(mf[[1]]), length.out=300),
> >dnorm(seq(min(mf[[1]]), max(mf[[1]]), length.out=300),
> > mean(mf[[1]]), sd(mf[[1]])),type="l", col="dark gray")
> >
> > Even trying to add an abline to the plot doesn't work.
> >
> > What am I missing?
> >
> > cheers
> >
> > Andy
> >
> > --
> > Andrew Halford Ph.D
> > Senior Coastal Fisheries Scientist
> > Pacific Community | Communauté du Pacifique CPS – B.P. D5 | 98848 Noumea,
> > New Caledonia | Nouméa, Nouvelle-Calédonie
> >
> > [[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.
>
--
Andrew Halford Ph.D
Senior Coastal Fisheries Scientist
Pacific Community | Communauté du Pacifique CPS – B.P. D5 | 98848 Noumea,
New Caledonia | Nouméa, Nouvelle-Calédonie
[[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.
[R] multiple graphs on one plot
Hi Listers
I've been trying to make a single graphic that has frequency histograms for
male and female mud crabs displayed side by side (such as when using the
beside=TRUE command for barplots). I then want to display a normal
distribution on top of the male and female histograms.
I have been using the multhist command in Plotrix to generate the
histograms without too much problem, but I cannot get the normal
distributions to plot up on the same graph.
Histograms plot
mf <-
list(lf_crabs$cw[lf_crabs$sex=='female'],lf_crabs$cw[lf_crabs$sex=='male'])
multhist(mf, xlab="CW", ylab="Frequency", ylim=c(0,100),main="All Measured
Crabs", col=c("dark gray", "light gray"),
breaks=seq(90,210, by=10),beside=TRUE,space=c(0,0.5))
legend("topright", c("Females", "Males"), fill=c("dark gray", "light gray"))
Then I try to add a normal distribution curve just to the female data but I
cant get the output to plot
points(seq(min(mf[[1]]), max(mf[[1]]), length.out=300),
dnorm(seq(min(mf[[1]]), max(mf[[1]]), length.out=300),
mean(mf[[1]]), sd(mf[[1]])),type="l", col="dark gray")
Even trying to add an abline to the plot doesn't work.
What am I missing?
cheers
Andy
--
Andrew Halford Ph.D
Senior Coastal Fisheries Scientist
Pacific Community | Communauté du Pacifique CPS – B.P. D5 | 98848 Noumea,
New Caledonia | Nouméa, Nouvelle-Calédonie
[[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.
[R] spatial analysis using quickPCNM
Hi Listers,
I posted this message to the R-sig-ecology group last Friday but have not
had a response hence my post here.
I have been trying to run spatial analyses on a fish community dataset.
My fish dataset has 114 species(variables) x 45 sites
My spatial dataset has the Lat and Long values for each site, converted to
cartesian coordinates
> fish <- read.table(file.choose())> latlong <- read.table(file.choose()) >
> fish.h <- decostand (fish, "hellinger")> fish.PCNM.quick <-
> quickPCNM(fish.h,latlong)
Truncation level = 639.5348
Time to compute PCNMs = 0.82 sec Error in if (temp2.test[1, 5] <=
alpha) { : argument is of length zeroTiming stopped at: 1.06 0.05 1.19
I do not understand the error message coming up and would appreciate
some advice.
Andy
--
Andrew Halford Ph.D
Research Scientist (Kimberley Marine Parks)
Dept. Parks and Wildlife
Western Australia
Ph: +61 8 9219 9795
Mobile: +61 (0) 468 419 473
[[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] grImport/ghostscript problems
Hi Again, I still cannot get the PostscriptTrace to work. The error message remains the same. I have run the following with no problems Sys.setenv(R_GSCMD=c:/Program Files(x86)/gs/gs9.07/bin/gswin32) but the call to Postscript Trace still fails PostScriptTrace(fish.ps) Error in PostScriptTrace(fish.ps) : status 127 in running command 'c:/Program Files(x86)/gs/gs9.07/bin/gswin32 -q -dBATCH -dNOPAUSE -sDEVICE=pswrite -sOutputFile=C:\Users\ahalford\AppData\Local\Temp\RtmpiIf7Zx\file8507fcd9d1 -sstdout=fish.ps.xml capturefish.ps' What am I doing wrong here? Andy On 14 August 2013 12:10, Andrew Halford andrew.half...@gmail.com wrote: Hi Listers I have been trying to import a .ps graphic file into R using the grImport package but I keep getting the following error message Error in PostScriptTrace(fish.ps) : status 127 in running command 'gswin32c.exe -q -dBATCH -dNOPAUSE -sDEVICE=pswrite -sOutputFile=C:\Users\ahalford\AppData\Local\Temp\Rtmp6BOVDe\fileffc30613d6 -sstdout=fish.ps.xml capturefish.ps' Any advice appreciated. Andy -- Andrew Halford Ph.D Adjunct Research Scientist University of Guam Curtin University Ph: +61 (0) 468 419 473 -- Andrew Halford Ph.D Adjunct Research Scientist University of Guam Curtin University Ph: +61 (0) 468 419 473 [[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] grImport/ghostscript problems
Hi Listers I have been trying to import a .ps graphic file into R using the grImport package but I keep getting the following error message Error in PostScriptTrace(fish.ps) : status 127 in running command 'gswin32c.exe -q -dBATCH -dNOPAUSE -sDEVICE=pswrite -sOutputFile=C:\Users\ahalford\AppData\Local\Temp\Rtmp6BOVDe\fileffc30613d6 -sstdout=fish.ps.xml capturefish.ps' Any advice appreciated. Andy -- Andrew Halford Ph.D Adjunct Research Scientist University of Guam Curtin University Ph: +61 (0) 468 419 473 [[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] extracting x and y coordinates from rnorm
Hi Listers, I have generated a random normal distribution using rnorm y - rnorm(180,116.7,5.4) How do I get a list of the x values that correspond to each y value that was generated. If I do a call to y, then I get the y-values but I also need the x values so I can generate a plot in another package I am using. Andy -- Andrew Halford Ph.D Adjunct Research Scientist University of Guam Curtin University Ph: +61 (0) 468 419 473 [[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] loop for calculating 1-se in rpart
Hi Listers
I need to calculate and then plot a frequency histogram of the best tree
calculated using the 1-se rule. I have included some code that has worked
well for me in the past but it was only for selecting the minimum
cross-validation error. I include the code for my model, some relevant
output and the code for selecting and plotting the frequency histogram of
minimum xerror.
Here is the output that is being referenced in the code below
Regression tree:
rpart(formula = chbiomsq ~ HC + BC + POC + RUG + Depth + Exp +
DFP + FI + LAT, data = ch, method = anova, control =
rpart.control(minsplit = 10,
cp = 0.01, xval = 10))
Variables actually used in tree construction:
[1] BCDepth DFP Exp
Root node error: 47456/99 = 479.35
n= 99
CP nsplit rel error xerror xstd
1 0.344626 0 1.0 1.02074 0.139585
2 0.179054 1 0.65537 0.76522 0.107470
3 0.072037 2 0.47632 0.68115 0.092627
4 0.063469 3 0.40428 0.67320 0.094830
5 0.036190 4 0.34081 0.58516 0.096726
6 0.034677 5 0.30462 0.56747 0.074953
7 0.018219 6 0.26995 0.52036 0.069447
8 0.017033 7 0.25173 0.54695 0.074249
9 0.010672 8 0.23469 0.55741 0.075119
10 0.01 9 0.22402 0.55756 0.074017
#Here is the code to run the model 50 times and take the minimum xerror
each time, followed by the histogram. However to calculate the appropriate
tree size which satisfies the 1-se rule I need to select minimum xerror
xerror (minimum xerror + xstd error) and which also has the smallest
nsplit.
cp50 - replicate(50,{ fit1 -
rpart(chbiomsq~HC+BC+POC+RUG+Depth+Exp+DFP+FI+LAT, data=ch,method=anova,
control=rpart.control(minsplit=10,cp=0.01,xval=10));x=printcp(fit1);x[which.min(x[,xerror]),nsplit]})
hist(cp50,main=optimal splits for tree,xlab= no. of optimal tree
splits, ylab= frequency)
Any help appreciated.
Andy
--
Andrew Halford Ph.D
Adjunct Research Scientist
University of Guam Curtin University
Ph: +61 (0) 468 419 473
[[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] plotting issues with PCA
Hi Listers, This has a simple answer but it has been eluding me nonetheless. I have been building a PCA plot from scratch with the ability to plot predefined groups in different colors. This has worked fine but when I try to get a polygon drawn around each of the groups it is not recognising my colour file correctly and is only printing the first colour in the filecode is below site.codings - c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,2,5,4,4,4,4,4,3,3,3,6,6,6,6,6,6,6,6,5,5,5,5) names(site.codings) - c( WM1, WM2, WM3, NM1, NM2, NM3, KH1, KH2, KH3, LM1 ,LM2 ,LM3, DB1 ,DB2 ,DB3, DM1 , DM2 , DM3 , FI1, FI2, BKI1, BKI2, BKO1, BKO2, BKO3, SUR1,MI1,MI2,MI3,BHE1,BHE2,BHE3,BHW1,BHW2,BHW3,HAL1,HAL2,HAL3,HAL4,HAL5,HAL6,HAL7,DOH1,DOH2,DOH3,DOH4,DOH5) fish.pca -rda(fish.sqrt.h) fish.site - scores(fish.pca,display=sites,scaling=3) fish.spp - scores(fish.pca,display=species,scaling=3)[omanfish.mrt.indval$pval=0.05,] graph - plot(fish.pca,display=c(sites,species),type=n,scaling=3) plotcolor - c(red,green,blue,aquamarine,magenta,yellow)[site.codings] points(fish.site,pch=21,bg=plotcolor,cex=1.2) #up to this point all works well but when I try to draw the polygons I cant get the lines to colour code the same way as the points did ordiellipse(graph,site.codings,kind=sd,conf=0.90,draw=polygon) I see there is a command called show.groups but I cant work out how to use it to access the plotcolor file. Any help appreciated. -- Andrew Halford Ph.D Associate Research Scientist Marine Laboratory University of Guam Ph: +1 671 734 2948 [[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] plot 3d info in 2d
Hi Listers, Is it possible to produce an ordination plot in 2d, where bubbles represent the location of sites (this part is easy enough) and the size of the bubbles is proportional to the sites location in 3d space (I am stuck on this option). So sites that are very near the 2d plane of the xy axes would be larger while sites that are actually further away in 3 d space would be proportionally smaller. any help/advice appreciated Andy -- Andrew Halford Ph.D Associate Research Scientist Marine Laboratory University of Guam Ph: +1 671 734 2948 [[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] reflecting a PCA biplot
Hi Listers, I am trying to reflect a PCA biplot in the x-axis (i.e. PC1) but am not having much success. In theory I believe all I need to do is multiply the site and species scores for the PC1 by -1, which would effectively flip the biplot. I am creating a blank plot using the plot command and accessing the results from a call to rda. I then use the calls to scores to obtain separate site and species coordinates and I have worked out how to multiply the appropriate PC1 scores by -1 to create the site and species scores I want. However I am not sure how to change the call to plot which accesses the results of the call to rda to draw the blank plot. The coordinates it is accessing are for the unreflected ordination and this does not match the new site and species scores that I have. fish.pca - rda(fish.hel) fish.site - scores(fish.pca,display=sites,scaling=3) fish.spp - scores(fish.pca,display=species,scaling=3) fish.site[,PC1] - -1*(fish.site[,PC1]) fish.spp[,PC1] - -1*(fish.spp[,PC1]) graph - plot(fish.pca,display=c(sites,species),type=n,scaling=3) # how do I get the plot to draw up the blank display based on the reversed site and species scores? Any help appreciated. cheers Andy -- Andrew Halford Ph.D Associate Research Scientist Marine Laboratory University of Guam Ph: +1 671 734 2948 [[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] ordiellipse
Hi Users, I am trying to plot 95% confidence ellipses around the 4 groupings found in a PCA ordination and identified a-priori by clustering. I am using ordiellipse to plot the ellipses and it works fine except one of my groups contains only 2 sites and I cannot get an ellipse around them. I'm assuming that 2 points is not enough to perform the relevant calculations here, however I would like to plot one if I could, if only for the sake of pictorial consistency. Andy -- Andrew Halford Ph.D Associate Research Scientist Marine Laboratory University of Guam Ph: +1 671 734 2948 [[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] which multivariate regression?
Hi R-Users, I have a student doing work with lionfish and she has been trying to analyse a multivariate dataset to see what variables/factors are influencing the behaviour of lionfish. We have attempted a number of analyses, including rpart, relimpo and standard linear regression but we are not having much luck with quality output. The data is very non-normal and we would appreciate some advice on the best way to go about analysing it. Kathy has provided a synopsis below along with part of the dataset below. Any help/advice appreciated. I am stuck in a problem with a dataset on a behavior study on Indo-Pacific lionfish *Pterois volitans*. The idea is to find out whether lionfish behave differently at different locations and times of day and whether these differences can be accounted for by any of the explanatory variables measured. My response variable is a series of behavior categories: (1) rest, (2) passive hunting and (3) active hunting. I have chosen to treat them individually because each one has a different biological importance, so basically I am trying to come up with an answer for 3 response variables. Measurement for these behavior categories is proportion of time (10 minute observation) spent at the activity described and values range from 0 to 1. Explanatory variables are a mix of categorical and continuous variables and are six: Region (Guam and Philippines), Hours after Sunrise, Habitat (5 categories), Weather (3 categories), Current (3 categories) and Lionfish Size (cm). The following is an example of the dataset for response variable Rest (R) R REG HAS HAB WE CU SI 0.05 0 11.017 Artificial 2 0 10 0.05 0 0.5667 Rock_boulder_cave 1 1 11 0.05 0 9.1333 Artificial 1 1 18 0.1 0 4.2 Sand_rubble 1 2 20 0.1 0 9.1333 Rock_boulder_cave 1 2 10 0.1 0 9.6 Sand_rubble 0 0 7 0.1 0 0.7833 Rock_boulder_cave 1 0 31 0.1 0 1.2833 Artificial 1 0 20 0.1 0 10.867 Coral 1 0 22 0.15 0 10.417 Coral 0 1 27 0.2 0 3.4667 Rock_boulder_cave 0 0 8 0.2 0 1.2333 Rock_boulder_cave 1 0 25 0.45 1 11.683 Coral 2 0 15 0.5 1 11.017 Artificial 1 2 14 0.5 1 11.917 Artificial 0 0 14 0.5 1 9.5333 Artificial 1 0 24 0.5 1 9.8333 Artificial 1 0 15 0.5 1 11.583 Rock_boulder_cave 1 1 29 0.53 1 5.9167 Coral 1 1 15 0.6 1 11.017 Artificial 1 2 17 0.6 1 9.7833 Rock_boulder_cave 0 0 12 0.6 1 4.6833 Sand_rubble 2 0 14 0.6 1 5.0167 Rock_boulder_cave 2 0 16 0.6 1 3.1833 Artificial 2 1 19 0.65 1 5.25 Coral 2 0 15 0.65 1 9.6333 Sand_rubble 1 1 17 As you can see here I have converted categorical variables region, current and weather to numerical; region because it can be expressed in binary form and the other two because they represent a quantity. For habitat I have created a dummy variable based on deviation coding, and introduced it as a variable in my model. Total sample size is 357, of which each sample is an observation at a particular time of day. A histogram of my response variable is not normally distributed and has a bit of a U-shape with lots of 0s and 1s, which means the animal was either completely engaged in that activity during the 10 min. observation or didn't show it at all. I have tried a series of transformations to normalize but have been unsuccessful (log, log(x+1), ln, sqrt, fourth root). What type of analyses have I tried? (1) Regression trees. Using categorical variables as categorical without changing into numerical. This was coded with package rpart and is the preferred analyses due to ease of interpretation. The response variable was untransformed and the distribution chosen Poisson. Result was a tree with immediately increasing error (cp) which picked 0 splits as the best tree. (2) Multiple regression Tried using package relaimpo to obtain a classification on the importance of explanatory variables. Used different transformations to analyze residuals and in all cases obtained a weird looking set of residuals with a portion normally distributed and another portion clustered to the side, giving the whole graph a clear trend (my guess is these are all the 1s and 0s in the data). I also tried non-linear regressions (glm) with package pscl (Poisson, negative binomial and zero inflated negative binomial. In all cases fit seemed adequate but variance explained was very small and coefficients estimated for my EVs very low. Any ideas??? I have lastly used Primer to analyze the response variable in response to each EV individually. That works well but limits my conclusions and doesn't allow me to account for variation in one of the EVs affecting others. I appreciate any help I can get, -- Andrew Halford Ph.D Associate Research Scientist Marine Laboratory University of Guam Ph: +1 671 734 2948
[R] goodness-of-fit test
Hi All, I have a dataset consisting of abundance counts of a fish and I want to test if my data are poisson in distribution or normal. My first question is whether it is more appropriate to model my data according to a poisson distribution (if my test says it conforms) or use transformed data to normalise the data distribution? I have been using the vcd package gf-goodfit(Y,type= poisson,method= MinChisq) but i get the following error message Warning message: In optimize(chi2, range(count)) : NA/Inf replaced by maximum positive value I then binned my count data to see if that might help V1 V2 1 5 34 2 10 30 3 15 10 4 20 8 5 25 7 6 30 0 7 35 3 8 40 2 9 45 3 10 50 1 11 55 0 12 60 1 but still received an error message Goodness-of-fit test for poisson distribution X^2 df P( X^2) Pearson 2573372 330 Warning message: In summary.goodfit(gf) : Chi-squared approximation may be incorrect Am I getting caught out because of zero counts or frequencies in my data? Andy -- Andrew Halford Ph.D Associate Research Scientist Marine Laboratory University of Guam Ph: +1 671 734 2948 [[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] penalized regression analysis
Hi All, I am using the package 'penalized' to perform a multiple regression on a dataset of 33 samples and 9 explanatory variables. The analysis appears to have performed as outlined and I have ended up with 4 explanatory variables and their respective regression coefficients. What I am struggling to understand is where do I get the variance explained information from and how do I determine the relative importance of the 4 variables selected? It does not appear to be a part of the penalized procedure. I submit the final call to 'penalized' with the estimated values of lambda1 and lambda2 fitfinal - penalized(CHAB~.,data=chabun,lambda1=356.0856,lambda2=3.458605,model = linear,steps=1,standardize = TRUE) # nonzero coefficients: 5 fitfinal Penalized linear regression object 10 regression coefficients of which 5 are non-zero Loglikelihood = -154.1055 L1 penalty = 4944.889 at lambda1 = 356.0856 L2 penalty = 234.7781 at lambda2 = 3.458605 coefficients (fitfinal) (Intercept)BC POC EXPFI 4.685739e+01 2.074521e-01 1.079459e-01 -1.373058e-05 -2.295339e+00 cheers Andy -- Andrew Halford Ph.D Associate Research Scientist Marine Laboratory University of Guam Ph: +1 671 734 2948 [[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] repeating an analysis
Hi All,
For those still interested the code submitted by Phil (see below) worked a
treat and produced a vector with the optimal 'nsplit' collated from 50 runs
of the rpart model. I then produced a histogram for the vector called answer
and chose my modal number for nsplits from which I had my appropriate tree
size.
many thanks to all who answered
Andy
On Wed, Oct 13, 2010 at 10:30 AM, Phil Spector spec...@stat.berkeley.eduwrote:
Andrew -
I think
answer = replicate(50,{fit1 - rpart(CHAB~.,data=chabun, method=anova,
control=rpart.control(minsplit=10,
cp=0.01, xval=10));
x = printcp(fit1);
x[which.min(x[,'xerror']),'nsplit']})
will put the numbers you want into answer, but there was no reproducible
example to test it on. Unfortunately, I don't know of any way to surpress
the printing from printcp().
- Phil Spector
Statistical Computing Facility
Department of Statistics
UC Berkeley
spec...@stat.berkeley.edu
On Wed, 13 Oct 2010, Andrew Halford wrote:
Hi All,
I have to say upfront that I am a complete neophyte when it comes to
programming. Nevertheless I enjoy the challenge of using R because of its
incredible statistical resources.
My problem is this .I am running a regression tree analysis using
rpart and I need to run the calculation repeatedly (say n=50 times) to
obtain a distribution of results from which I will pick the median one to
represent the most parsimonious tree size. Unfortunately rpart does not
contain this ability so it will have to be coded for.
Could anyone help me with this? I have provided the code (and relevant
output) for the analysis I am running. I need to run it n=50 times and
from
each output pick the appropriate tree size and post it to a datafile where
I
can then look at the frequency distribution of tree sizes.
Here is the code and output from a single run
fit1 - rpart(CHAB~.,data=chabun, method=anova,
control=rpart.control(minsplit=10, cp=0.01, xval=10))
printcp(fit1)
Regression tree:
rpart(formula = CHAB ~ ., data = chabun, method = anova, control =
rpart.control(minsplit = 10,
cp = 0.01, xval = 10))
Variables actually used in tree construction:
[1] EXP LAT POC RUG
Root node error: 35904/33 = 1088
n= 33
CP nsplit rel error xerrorxstd
1 0.539806 0 1.0 1.0337 0.41238
2 0.050516 1 0.46019 1.2149 0.38787
3 0.016788 2 0.40968 1.2719 0.41280
4 0.010221 3 0.39289 1.1852 0.38300
5 0.01 4 0.38267 1.1740 0.38333
Each time I re-run the model I will get a slightly different output. I
want
to extract the nsplit number corresponding to the lowest xerror for each
run
of the model (in this case it is for nsplit = 0) over 50 runs and then
look
at the distribution of nsplits after 50 runs.
Any help appreciated.
Andy
--
Andrew Halford
Associate Researcher
Marine Laboratory
University of Guam
Ph: +1 671 734 2948
[[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.
--
Andrew Halford Ph.D
Associate Research Scientist
Marine Laboratory
University of Guam
Ph: +1 671 734 2948
[[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] repeating an analysis
thanks Phil, I have your solution and another which I will attempt in the
next day or so and will post results to the list then.
cheers
andy
On Wed, Oct 13, 2010 at 10:30 AM, Phil Spector spec...@stat.berkeley.eduwrote:
Andrew -
I think
answer = replicate(50,{fit1 - rpart(CHAB~.,data=chabun, method=anova,
control=rpart.control(minsplit=10,
cp=0.01, xval=10));
x = printcp(fit1);
x[which.min(x[,'xerror']),'nsplit']})
will put the numbers you want into answer, but there was no reproducible
example to test it on. Unfortunately, I don't know of any way to surpress
the printing from printcp().
- Phil Spector
Statistical Computing Facility
Department of Statistics
UC Berkeley
spec...@stat.berkeley.edu
On Wed, 13 Oct 2010, Andrew Halford wrote:
Hi All,
I have to say upfront that I am a complete neophyte when it comes to
programming. Nevertheless I enjoy the challenge of using R because of its
incredible statistical resources.
My problem is this .I am running a regression tree analysis using
rpart and I need to run the calculation repeatedly (say n=50 times) to
obtain a distribution of results from which I will pick the median one to
represent the most parsimonious tree size. Unfortunately rpart does not
contain this ability so it will have to be coded for.
Could anyone help me with this? I have provided the code (and relevant
output) for the analysis I am running. I need to run it n=50 times and
from
each output pick the appropriate tree size and post it to a datafile where
I
can then look at the frequency distribution of tree sizes.
Here is the code and output from a single run
fit1 - rpart(CHAB~.,data=chabun, method=anova,
control=rpart.control(minsplit=10, cp=0.01, xval=10))
printcp(fit1)
Regression tree:
rpart(formula = CHAB ~ ., data = chabun, method = anova, control =
rpart.control(minsplit = 10,
cp = 0.01, xval = 10))
Variables actually used in tree construction:
[1] EXP LAT POC RUG
Root node error: 35904/33 = 1088
n= 33
CP nsplit rel error xerrorxstd
1 0.539806 0 1.0 1.0337 0.41238
2 0.050516 1 0.46019 1.2149 0.38787
3 0.016788 2 0.40968 1.2719 0.41280
4 0.010221 3 0.39289 1.1852 0.38300
5 0.01 4 0.38267 1.1740 0.38333
Each time I re-run the model I will get a slightly different output. I
want
to extract the nsplit number corresponding to the lowest xerror for each
run
of the model (in this case it is for nsplit = 0) over 50 runs and then
look
at the distribution of nsplits after 50 runs.
Any help appreciated.
Andy
--
Andrew Halford
Associate Researcher
Marine Laboratory
University of Guam
Ph: +1 671 734 2948
[[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.
--
Andrew Halford Ph.D
Associate Researcher Scientist
Marine Laboratory
University of Guam
Ph: +1 671 734 2948
[[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] repeating an analysis
Hi All, I have to say upfront that I am a complete neophyte when it comes to programming. Nevertheless I enjoy the challenge of using R because of its incredible statistical resources. My problem is this .I am running a regression tree analysis using rpart and I need to run the calculation repeatedly (say n=50 times) to obtain a distribution of results from which I will pick the median one to represent the most parsimonious tree size. Unfortunately rpart does not contain this ability so it will have to be coded for. Could anyone help me with this? I have provided the code (and relevant output) for the analysis I am running. I need to run it n=50 times and from each output pick the appropriate tree size and post it to a datafile where I can then look at the frequency distribution of tree sizes. Here is the code and output from a single run fit1 - rpart(CHAB~.,data=chabun, method=anova, control=rpart.control(minsplit=10, cp=0.01, xval=10)) printcp(fit1) Regression tree: rpart(formula = CHAB ~ ., data = chabun, method = anova, control = rpart.control(minsplit = 10, cp = 0.01, xval = 10)) Variables actually used in tree construction: [1] EXP LAT POC RUG Root node error: 35904/33 = 1088 n= 33 CP nsplit rel error xerrorxstd 1 0.539806 0 1.0 1.0337 0.41238 2 0.050516 1 0.46019 1.2149 0.38787 3 0.016788 2 0.40968 1.2719 0.41280 4 0.010221 3 0.39289 1.1852 0.38300 5 0.01 4 0.38267 1.1740 0.38333 Each time I re-run the model I will get a slightly different output. I want to extract the nsplit number corresponding to the lowest xerror for each run of the model (in this case it is for nsplit = 0) over 50 runs and then look at the distribution of nsplits after 50 runs. Any help appreciated. Andy -- Andrew Halford Associate Researcher Marine Laboratory University of Guam Ph: +1 671 734 2948 [[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.
