gmane.comp.lang.r.ecology

Re: Poisson regression

Highland Statistics Ltd — 2012-05-24T15:40:59







Check whether your Clupeidae variable contains non-integers.......

Alain

Re: Poisson regression

Gavin Simpson — 2012-05-24T11:16:01

I don't know if this is the problem or not, but you can't call the model
that way. if you don't name arguments then you must list them in the
order the function expects. You don;t say which `gam()` you use but
assuming it is `mgcv:::gam()` then the second argument is `data` and you
passed it a function `poisson`.

Does it work if you do:

model1<-gam(Clupeidae~s(depth)+s(temperature)+s(salinity),
            family = poisson)

?

Really though you should be passing it both a data and a family
argument. Assuming your data are in object named `mydata` then:

model1<-gam(Clupeidae~s(depth)+s(temperature)+s(salinity),
            data = mydata, family = poisson)

would be the correct way to work with the function.

G

Poisson regression

Mahnaz Rabbaniha — 2012-05-24T07:11:06

Dear all


to find relation between non-normal response with independents variable i
use this code:

model1<-gam(Clupeidae~s(depth)+s(temperature)+s(salinity),poisson)

the result shown is:

There were 50 or more warnings (use warnings() to see the first 50)
Warning messages:
1: In dpois(y, mu, log = TRUE) : non-integer x = 2.079542 Error in
cat(list(...), file, sep, fill, labels, append) :
  argument 2 (type 'list') cannot be handled by 'cat'


what is meaning it?

do i allow use it for continue?

thanks

Re: PCA as a predictive model

Marc Taylor — 2012-05-23T10:42:49

Thank you Bob - you and Jari seem to be of consensus here :-)

I will have to double check that what I am doing really gives the same
result as predict.prcomp. My problem is that I have set up my PCA in a
slightly different way than you have - or for that matter, different from
many of the ordination examples in R -

My matrix columns are the "samples" and the rows are the "variables" - just
the transpose of what is typical. So I have been calculating the PC
loadings based on the sample covariances and not the variable covariances.
If you don't mind my excursion from ecology, I'll explain that my data
consists of measured light spectra (350-800nm) with a value at each nm.
Thus my matrix consists of 451 rows and each column is a sample. I have
been applying my scaling to the samples and not the variables. So, the
length of my centers is the length of the number of samples. I thus run
into problems when I want to predict from a newdata that contains a
different number of samples. In the end, I think I am getting reasonable
predictions by doing this in the way that I described here:
http://stats.stackexchange.com/questions/28916/can-empirical-orthogonal-function-eof-analysis-be-used-as-a-predictive-model/28986#28986
It would, however, be comforting to reproduce this prediction in the
standard way... I'll keep experimenting.
Cheers, Marc

On Wed, May 23, 2012 at 11:15 AM, Bob O'Hara <bohara-ePdyZ8bJA/Uw22byECBumQ< at >public.gmane.org> wrote:

[[alternative HTML version deleted]]

Re: PCA as a predictive model

Bob O'Hara — 2012-05-23T09:15:31

I'll pretend to be Jari for a moment. :-)

PCA just scales and rotates the data in cunning ways, so with the new 
data you need to scale and rotate it in the same way. If you scale the 
values first then you've already changed the scaling.

What you need to do is either do PCA on the raw data or scale the new 
data using the mean and varianes of the old data.

library(MASS)

NVar=5; NObs=50
Sigma=matrix(c(
  10,0.2,   0, 0,0.4,
0.2,   5,0.1, 0,0.6,
    0,0.1,1.0, 0.2,0,
    0,   0,0.2, 5, 0,
0.4,0.6,  0, 0,1), nrow=5)

# simulate data
Data=mvrnorm(NObs, rnorm(NVar), Sigma=Sigma)
# Do PCA on scaled data
Data.Sc=scale(Data)
PC=princomp(Data.Sc)

# Simulate new data
NewData=mvrnorm(10, rnorm(NVar), Sigma=Sigma)
# Do PCA on new data. First do it wrong...
PC.wrong=predict(PC, newdata=scale(NewData))

# Now scale correctly

NewData.Sc=scale(NewData, center=attr(Data.Sc, "scaled:center"), 
scale=attr(Data.Sc, "scaled:scale")
PC.right=predict(PC, newdata=NewData.Sc)

HTH

Bob

Re: PCA as a predictive model

Marc Taylor — 2012-05-23T09:10:34

Excellent - thanks Jari.
Cheers,
Marc

On Wed, May 23, 2012 at 11:05 AM, Jari Oksanen <jari.oksanen-jSDM6y+8MbU< at >public.gmane.org> wrote:


[[alternative HTML version deleted]]

Re: PCA as a predictive model

Jari Oksanen — 2012-05-23T09:05:56

Marc,

I see no danger. When you predict with 'newdata', that 'newdata' should have no influence on the original PCA. You can use predict(..., newdata) for cross validation.  If you look at the code of stats:::predict.prcomp, you see that the prediction is only (i) scaling and centring your new data with the means and scale of the original data, and (ii) rotation of these scaled values to the original PCs. This is the last line of stats:::predict.prcomp:

scale(newdata, object$center, object$scale) %*% object$rotation

The PCA result 'object' saves the original 'center' and 'scale' (which can be FALSE or SD of variables) of your original variables, and applies those to your 'newdata' before rotating to PCs.

Cheers, Jari Oksanen
________________________________
From: Marc Taylor [marchtaylor-Re5JQEeQqe8AvxtiuMwx3w< at >public.gmane.org]
Sent: 23 May 2012 11:55
To: Jari Oksanen
Cc: r-sig-ecology-0bNBQ1PAWB4BXFe83j6qeQ< at >public.gmane.org
Subject: Re: [R-sig-eco] PCA as a predictive model

Hi Jari - one more question if you don't mind. Since the weights of the PCs are related to the the amount of variance that they explain in the original data - is it problematic to predict the PC scores with a second data set that has a different amount of variance (e.g. due to differing number of samples)? In both the 1st and 2nd data sets I have been using scaled values for the variables (mean=0 and sd=1 for each sample).
Cheers,
Marc


On Wed, May 23, 2012 at 9:59 AM, Marc Taylor <marchtaylor-Re5JQEeQqe8AvxtiuMwx3w< at >public.gmane.org<mailto:marchtaylor-Re5JQEeQqe8AvxtiuMwx3w< at >public.gmane.org>> wrote:
Hi Jari,

That's good to hear - I hadn't made the connection to cca/rda. This will help me find pertinent literature as well.

Many thanks,
Marc


On Wed, May 23, 2012 at 9:33 AM, Jari Oksanen <jari.oksanen-jSDM6y+8MbU< at >public.gmane.org<mailto:jari.oksanen-jSDM6y+8MbU< at >public.gmane.org>> wrote:
Marc,

Basic R stats functions like prcomp have a predict method that can be used to "predict" (calculate) scores with 'newdata'. This is standard, and has been in R for ever. Most textbooks of multivariate analysis should handle this issue.

In community ecological context, see functions predict.cca and predict.rda in vegan. These take argument 'newdata' which can be new community data -- depending on what you want to predict (argument 'type'). Function calibrate.cca documented in the same help page can be used to predict values of constraining variables in constrained ordination (CCA, RDA) from community composition, also with 'newdata' communities.

Cheers, Jari Oksanen

________________________________________
From: r-sig-ecology-bounces-0bNBQ1PAWB4BXFe83j6qeQ< at >public.gmane.org<mailto:r-sig-ecology-bounces-0bNBQ1PAWB4BXFe83j6qeQ< at >public.gmane.org> [r-sig-ecology-bounces-0bNBQ1PAWB4BXFe83j6qeQ< at >public.gmane.org<mailto:r-sig-ecology-bounces-0bNBQ1PAWB4BXFe83j6qeQ< at >public.gmane.org>] on behalf of Marc Taylor [marchtaylor-Re5JQEeQqe8AvxtiuMwx3w< at >public.gmane.org<mailto:marchtaylor-Re5JQEeQqe8AvxtiuMwx3w< at >public.gmane.org>]
Sent: 23 May 2012 10:19
To: r-sig-ecology-0bNBQ1PAWB4BXFe83j6qeQ< at >public.gmane.org<mailto:r-sig-ecology-0bNBQ1PAWB4BXFe83j6qeQ< at >public.gmane.org>
Subject: [R-sig-eco] PCA as a predictive model

Hello R-sig-ecology group,

I was wondering if anyone is aware of an example where PCA is used as a
predictive model? A community analysis example might be to predict the PC
values of a sample given its community composition. I had thrown this
question up on a statistics forum (
http://stats.stackexchange.com/questions/28916/can-empirical-orthogonal-function-eof-analysis-be-used-as-a-predictive-model)
but have gotten hardly any response. I imagined that there are some folks
here that would have some insight into this problem.

Many thanks,
Marc

       [[alternative HTML version deleted]]

_______________________________________________
R-sig-ecology mailing list
R-sig-ecology-0bNBQ1PAWB4BXFe83j6qeQ< at >public.gmane.org<mailto:R-sig-ecology-0bNBQ1PAWB4BXFe83j6qeQ< at >public.gmane.org>
https://stat.ethz.ch/mailman/listinfo/r-sig-ecology



[[alternative HTML version deleted]]

Re: PCA as a predictive model

Marc Taylor — 2012-05-23T08:55:04

Hi Jari - one more question if you don't mind. Since the weights of the PCs
are related to the the amount of variance that they explain in the original
data - is it problematic to predict the PC scores with a second data set
that has a different amount of variance (e.g. due to differing number of
samples)? In both the 1st and 2nd data sets I have been using scaled values
for the variables (mean=0 and sd=1 for each sample).
Cheers,
Marc


On Wed, May 23, 2012 at 9:59 AM, Marc Taylor <marchtaylor-Re5JQEeQqe8AvxtiuMwx3w< at >public.gmane.org> wrote:


[[alternative HTML version deleted]]

Re: PCA as a predictive model

Marc Taylor — 2012-05-23T07:59:19

Hi Jari,

That's good to hear - I hadn't made the connection to cca/rda. This will
help me find pertinent literature as well.

Many thanks,
Marc

On Wed, May 23, 2012 at 9:33 AM, Jari Oksanen <jari.oksanen-jSDM6y+8MbU< at >public.gmane.org> wrote:


[[alternative HTML version deleted]]

Re: PCA as a predictive model

Jari Oksanen — 2012-05-23T07:33:06

Marc,

Basic R stats functions like prcomp have a predict method that can be used to "predict" (calculate) scores with 'newdata'. This is standard, and has been in R for ever. Most textbooks of multivariate analysis should handle this issue. 

In community ecological context, see functions predict.cca and predict.rda in vegan. These take argument 'newdata' which can be new community data -- depending on what you want to predict (argument 'type'). Function calibrate.cca documented in the same help page can be used to predict values of constraining variables in constrained ordination (CCA, RDA) from community composition, also with 'newdata' communities.

Cheers, Jari Oksanen

________________________________________
From: r-sig-ecology-bounces-0bNBQ1PAWB4BXFe83j6qeQ< at >public.gmane.org [r-sig-ecology-bounces-0bNBQ1PAWB4BXFe83j6qeQ< at >public.gmane.org] on behalf of Marc Taylor [marchtaylor-Re5JQEeQqe8AvxtiuMwx3w< at >public.gmane.org]
Sent: 23 May 2012 10:19
To: r-sig-ecology-0bNBQ1PAWB4BXFe83j6qeQ< at >public.gmane.org
Subject: [R-sig-eco] PCA as a predictive model

Hello R-sig-ecology group,

I was wondering if anyone is aware of an example where PCA is used as a
predictive model? A community analysis example might be to predict the PC
values of a sample given its community composition. I had thrown this
question up on a statistics forum (
http://stats.stackexchange.com/questions/28916/can-empirical-orthogonal-function-eof-analysis-be-used-as-a-predictive-model)
but have gotten hardly any response. I imagined that there are some folks
here that would have some insight into this problem.

Many thanks,
Marc

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-ecology mailing list
R-sig-ecology-0bNBQ1PAWB4BXFe83j6qeQ< at >public.gmane.org
https://stat.ethz.ch/mailman/listinfo/r-sig-ecology

PCA as a predictive model

Marc Taylor — 2012-05-23T07:17:35

Hello R-sig-ecology group,

I was wondering if anyone is aware of an example where PCA is used as a
predictive model? A community analysis example might be to predict the PC
values of a sample given its community composition. I had thrown this
question up on a statistics forum (
http://stats.stackexchange.com/questions/28916/can-empirical-orthogonal-function-eof-analysis-be-used-as-a-predictive-model)
but have gotten hardly any response. I imagined that there are some folks
here that would have some insight into this problem.

Many thanks,
Marc

[[alternative HTML version deleted]]

Re: gam result

ONKELINX, Thierry — 2012-05-22T10:40:50

Given the failry complex model you want to fit on a relative small dataset, you should better seek some local (face-to-face) advise on your problem.

Best regards,

ir. Thierry Onkelinx
Instituut voor natuur- en bosonderzoek / Research Institute for Nature and Forest
team Biometrie & Kwaliteitszorg / team Biometrics & Quality Assurance
Kliniekstraat 25
1070 Anderlecht
Belgium
+ 32 2 525 02 51
+ 32 54 43 61 85
Thierry.Onkelinx-MuCAU19RsmE< at >public.gmane.org
www.inbo.be

To call in the statistician after the experiment is done may be no more than asking him to perform a post-mortem examination: he may be able to say what the experiment died of.
~ Sir Ronald Aylmer Fisher

The plural of anecdote is not data.
~ Roger Brinner

The combination of some data and an aching desire for an answer does not ensure that a reasonable answer can be extracted from a given body of data.
~ John Tukey


-----Oorspronkelijk bericht-----
Van: r-sig-ecology-bounces-0bNBQ1PAWB4BXFe83j6qeQ< at >public.gmane.org [mailto:r-sig-ecology-bounces-0bNBQ1PAWB4BXFe83j6qeQ< at >public.gmane.org] Namens Mahnaz Rabbaniha
Verzonden: maandag 21 mei 2012 18:33
Aan: r-sig-ecology
Onderwerp: [R-sig-eco] gam result

Dear all

for finding the relation between this variable i use gam ( in base of
non-normality)

[1] "temperature" "salinity"    "depth"       "Clupeidae"


model<-gam(Clupeidae~s(temperature)+s(salinity)+s(depth))


the result:

 Parametric coefficients:

                         Estimate        Std. Error       t value
    Pr(>|t|)

(Intercept)          0.40156        0.05467           7.345
2.12e-10 ***

Signif. codes:  0 '***'    0.001 '**'    0.01 '*'  0.05 '.' 0.1 ' ' 1



Approximate significance of smooth terms:

                               Edf        Ref.df         F           p-value

s(temperature)         4.321     5.326      1.470      0.207

s(salinity)               2.893      3.569      1.540       0.204

s(depth)                  1.000     1.000       0.015       0.903


R-sq.(adj) =    0.1       Deviance explained = 18.9%

GCV score = 0.28203       Scale est. = 0.25109        n = 84


this is my question;in base of R-sq.   or others parameters,is this
model robust ifor analyses?

also in others section this parameter is less than0.1


thanks all

_______________________________________________
R-sig-ecology mailing list
R-sig-ecology-0bNBQ1PAWB4BXFe83j6qeQ< at >public.gmane.org
https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
* * * * * * * * * * * * * D I S C L A I M E R * * * * * * * * * * * * *
Dit bericht en eventuele bijlagen geven enkel de visie van de schrijver weer en binden het INBO onder geen enkel beding, zolang dit bericht niet bevestigd is door een geldig ondertekend document.
The views expressed in this message and any annex are purely those of the writer and may not be regarded as stating an official position of INBO, as long as the message is not confirmed by a duly signed document.

gam result

Mahnaz Rabbaniha — 2012-05-21T16:33:00

Dear all

for finding the relation between this variable i use gam ( in base of
non-normality)

[1] "temperature" "salinity"    "depth"       "Clupeidae"


model<-gam(Clupeidae~s(temperature)+s(salinity)+s(depth))


the result:

 Parametric coefficients:

                         Estimate        Std. Error       t value
    Pr(>|t|)

(Intercept)          0.40156        0.05467           7.345
2.12e-10 ***

Signif. codes:  0 ‘***’    0.001 ‘**’    0.01 ‘*’  0.05 ‘.’ 0.1 ‘ ’ 1



Approximate significance of smooth terms:

                               Edf        Ref.df         F           p-value

s(temperature)         4.321     5.326      1.470      0.207

s(salinity)               2.893      3.569      1.540       0.204

s(depth)                  1.000     1.000       0.015       0.903


R-sq.(adj) =    0.1       Deviance explained = 18.9%

GCV score = 0.28203       Scale est. = 0.25109        n = 84


this is my question;in base of R-sq.   or others parameters,is this
model robust ifor analyses?

also in others section this parameter is less than0.1


thanks all

linear mixed-effect model

C Hess — 2012-05-19T13:26:47

Dear list,

I am working on a dendrochronological investigation with the use of a 
linear mixed-effect model.

I identified significant main and interactions effects in parameter 
estimation of fixed-effects

Has anyone a hint on how I can generate or receive results and 
statements like:

main/interaction effect of predictor A/AB explained X% of variation of 
growth as response variable.


thanks
CH

P.S.: thanks also for all replies concerning my previous questions. this 
is a great community.

Does weighting presence equal to background in species distribution modeling using glm in the stats package work correctly??

Marian K Talbert — 2012-05-18T14:26:28

In species distribution modeling where one uses a large sample of 
background points to capture background variation in 
presence\pseudo-absence or use\available models (0\1 response) it is 
frequently recommended that one weight the data so the sum of the absence 
weights is equal to the sum of presence weights so that the model isn?t 
swamped by an overwhelming and arbitrary number of background points.  My 
concern is I don't know that glm in the R stats package handles this 
correctly. 

When weights are set in glm they are ignored in the stepAIC and I'm 
wondering if this should be the case.  For example if I have a data set 
with 75 presence and 75 absence I fit glm.  If I repeat the 75 absences 10 
times say and then weigh the absences so total weight of absence equals 
total weight of presence fitting any individual glm to the two sets give 
almost the same answer (in coefficients, null and residual deviance but 
NOT in likelihood extracted using logLik).  Then when I use the step 
command in R to select the best from a series of models I get quite 
different answers for two reasons.  Step uses the likelihood which is not 
the same in the two models (and I'm not sure why given that they have the 
same deviance).  The second reason is that the step function in R uses AIC 
= - 2*log L + k * edf where edf is effective degrees of freedom.  It seems 
that if you're setting the weights then the the step function should 
adjust the degrees of freedom it uses accordingly.  I've tried adjusting 
the stepwise regression by multiplying the penalty term in step by 
sum(data$Weights)/length(data$Weights) where I generally set the penalty 
to 2 for AIC or  log(nrow(data)) for BIC.  Indeed when I use this approach 
I get much more similar predictions from the two models but because the 
likelihood is different I don't always get exactly the same predictions 
from the two models.  I was wondering if there is a statistical 
explanation for why the answers should be different and if not, why I 
can't get the correct answer.  I tried the same thing on STATA this 
afternoon and it was interesting to note that STATA did give the same 
likelihood for the two models and the same answer for stepwise selection. 

I'm aware that GLM in the standard stats package in R is intended to use 
weights only to specify an integer number of binomial trials under the 
assumption that the response is proportion of success but had read that 
the only consequence of using non-integer weights with a Bernoulli 
response was that starting values are difficult to choose and convergence 
to the maximum likelihood isn't guaranteed.  I was  willing to accept this 
possibility but it's seeming that the issues with using glm in this way 
might run more deep.  Another option would be to consider the 
presence\background as strata with different sampling probabilities and 
use the survey package svyglm but I'm not sure how to do the automatic 
model selection in that case. 

I've noticed that many papers and software programs use R code with non 
integer weights set in this or a similar way: BIOMOD uses the glm with 
weights and stepwise model selection without any consideration of the 
above issues, Phillips and Elith (Logistic Methods for Resource Selection 
Functions and Presence-Only Species Distribution Models), Guisan, Edwards 
and Hastie (Generalized linear and generalized additive models in studies 
of species distributions: setting the scene) though they use a custom 
stepwise selection based on the deviance rather than the likelihood. 

Anyway, if anyone has any insight on how to appropriately set weights in 
glm and use stepwise model selection, I'd really appreciate your 
recommendations. 

Marian Talbert

[[alternative HTML version deleted]]

adehabitatLT help in plot pathways

Conrado Galdino — 2012-05-17T19:19:07

Dear all,

I'm now using the package adehabitatLT to handle with animals movements.
Despite things are going ok  I cannot superimpose the various animals
pathways that plot.ltraj (or simple plot) function generates. My big
question is: Is it possible to do such procedure? Alternatively: Does some
of you can suggest some way to do this task?

Thanks in advance,
Cheers,

Re: gam

Bob O'Hara — 2012-05-17T17:44:30

It means that you don't have a lot of different values for one of your 
covariates. You can check which with

length(unique(depth))
length(unique(temperature))
length(unique(salinity))

Unless you have a lot of data, that the sensible thing to is probably to 
fit a straight line rather than a smooth curve. The other thing to do is 
to define the maximum flexibility of the curve, e.g.

model<-gam(Sillaginidae~s(temperature, k=6)+s(salinity)+s(depth))

I hope this helps.

Bob

gam

Mahnaz Rabbaniha — 2012-05-17T16:17:43

Dear all

i have done GAM for data ( they are non- normal and i find regression
between Sillaginidae with several hydrological factors),i done :

[1] "depth"        "temperature"  "salinity"     "Sillaginidae"

 > pairs(sc,panel=function(x,y){points(x,y);lines(lowess(x,y))})


but i take this message:

Error in smooth.construct.tp.smooth.spec(object, dk$data, dk$knots) :
  A term has fewer unique covariate combinations than specified
maximum degrees of freedom


please help me what can i do? or what is meaning this sentence?

thanks


Mahnaz

Re: Continuous (Non-Count) Skewed Data With Many Zeros

Mieke Zwart — 2012-05-17T10:22:52

Hi Alain,

Just wondering if the 2012 book you are talking about is the "Zero Inflated Models and Generalized Linear Mixed Models with R". You say to look at chapter 10 but this book has only 9 chapters according to the website. Is it the book you are talking about?

Cheers,

Mieke

 
-----Original Message-----
From: r-sig-ecology-bounces-0bNBQ1PAWB4BXFe83j6qeQ< at >public.gmane.org [mailto:r-sig-ecology-bounces-0bNBQ1PAWB4BXFe83j6qeQ< at >public.gmane.org] On Behalf Of Highland Statistics Ltd
Sent: 16 May 2012 14:09
To: r-sig-ecology-0bNBQ1PAWB4BXFe83j6qeQ< at >public.gmane.org
Subject: Re: [R-sig-eco] Continuous (Non-Count) Skewed Data With Many Zeros



Non-normality of your response variable is not a reason to apply a data transformation.
It all depends, and no sensible answer can be given. 15% of zeros can screw things up....but it is also possible that 80% of zeros comply with a regression or GLM. For a discussion with examples see Chapter 10 in our 2012 book.


Depends on the previous remark.....anything from linear regression to a zero inflated model for a continuous distributed response variable. 
There is just no simple answer possible. It all depends. But based on what you describe it will probably be something zero-inflated.

Alain

covariance structure issue with lme

Charles Determan Jr — 2012-05-16T16:19:01

Greetings,

I have previously tried the R mixed model mailing list but unfortunately I
remain
stuck with this question.  I certainly hope you will be able to shed some
light on this particular problem.  I am interested in ecology
investigations,
however this is just practice data to try and accomplish the modelling
parameters.
My sincere thanks for an assistance.

I have previously attempted to replicate the UN (unstructured),
CS (compound symmetry), and AR(1) (Autoregression)
covariance structures used in SAS PROC MIXED with the lme
function in this nlme package.  However, my efforts
have fallen short on replicating the *Variance Components (VC)*
structure.  I have read that it is also known as a diagonal structure.
 Below I have copied over all the models I have tried and their output
with no success.  Perhaps someone here will see my error or something
I have overlooked.  I have attached the data for this particular
model.  Thanks to all, I certainly cannot thank this help list enough.
 If you need any further information/clarification, please ask.

Cheers, Charles

library(nlme)

dat=read.table("C:/subset.csv"
,sep=",",header=TRUE, na.strings=".")
attach(dat)

dat34=dat[Group %in% c("3", "4"),]
attach(dat34)
liver34=within(dat34, {
      Group=factor(Group)
      Event_name=factor(Event_name)
      Died=factor(Died)
      ID=factor(ID)
})
attach(liver34)

##What is should be from SAS
#CV var
#Type 3 Tests of Fixed Effects
#Effect          NumDF DenDF F Value Pr > F
#Group                  1      22      0.25    0.6244
#Died                    1      22      6.55    0.0179
#Group*Died          1      22      4.43    0.0470

fit.1=lme(var~Group*Died,
      random=~1|ID,
      data=dat34)
anova(fit.1, type="marginal", adjustSigma=F)
#               numDF denDF   F-value p-value
#(Intercept)    1   101 227.58700  <.0001
#Group          1    22   0.18320  0.6728
#Died            1    22   3.63388  0.0698
#Group:Died   1    22   3.04103  0.0951

fit.2=lme(var~Group*Died,
      data=dat34,
      random=~1|ID/Died)
anova(fit.2, type="marginal", adjustSigma=F)
#               numDF denDF  F-value p-value
#(Intercept)    1   101 77.99004  <.0001
#Group          1    22  1.46275  0.2393
#Died             1    22  5.84535  0.0243
#Group:Died   1    22  3.04103  0.0951

fit.3=lme(var~Group*Died,
      random=list(ID=pdSymm(~Event_name)),
      data=dat34)
anova(fit.3, type="marginal", adjustSigma=F)
#               numDF denDF   F-value p-value
#(Intercept)    1   101 273.10918  <.0001
#Group          1    22   0.69692  0.4128
#Died             1    22   1.43316  0.2440
#Group:Died    1    22   5.74399  0.0255

fit.4=lme(var~Group*Died,
      random=list(ID=pdSymm(~Group)),
      data=dat34)
anova(fit.4, type="marginal", adjustSigma=F)
#               numDF denDF   F-value p-value
#(Intercept)    1   101 235.13889  <.0001
#Group          1    22   0.15878  0.6941
#Died             1    22   3.83253  0.0631
#Group:Died   1    22   3.01222  0.0966

fit.5=lme(var~Group*Died,
      random=list(ID=pdSymm(~Group)),
      data=dat34,
      weights=varIdent(form=~1|Event_name))
anova(fit.5, type="marginal", adjustSigma=F)
#               numDF denDF   F-value p-value
#(Intercept)    1   101 277.16705  <.0001
#Group          1    22   0.23901  0.6298
#Died            1    22   3.99283  0.0582
#Group:Died   1    22   3.23135  0.0860

fit.6=lme(var~Group*Died,
      random=list(ID=pdSymm(~Group)),
      data=dat34,
      weights=varIdent(form=~1|Event_name))
anova(fit.6, type="marginal", adjustSigma=F)
#               numDF denDF   F-value p-value
#(Intercept)    1   101 277.16705  <.0001
#Group          1    22   0.23901  0.6298
#Died             1    22   3.99283  0.0582
#Group:Died   1    22   3.23135  0.0860

fit.7=lme(var~(Group*Died),
      random=list(ID=pdCompSymm(~Died)),
      data=dat34)
anova(fit.7, type="marginal", adjustSigma=F)
#               numDF denDF  F-value p-value
#(Intercept)    1   101 85.83799  <.0001
#Group          1    22  1.60624  0.2183
#Died            1    22  4.71795  0.0409
#Group:Died   1    22  2.65379  0.1175

fit.8=lme(var~(Group*Died),
      data=dat34,
      random=~1|ID,
      corr=corSymm())
anova(fit.8, type="marginal", adjustSigma=F)
#               numDF denDF   F-value p-value
#(Intercept)    1   101 119.54403  <.0001
#Group          1    22   4.58972  0.0435
#Died            1    22   8.01715  0.0097
#Group:Died   1    22   5.27470  0.0315

fit.9=lme(var~(Group*Died),
      data=dat34,
      random=list(ID=pdDiag(~Group*Died)),
      corr=corSymm(, ~1|ID))
#  Error in lme.formula(var ~ (Group * Died), data = dat34, random =
list(ID = pdDiag(~Group *  :
#  nlminb problem, convergence error code = 1
#  message = iteration limit reached without convergence (9)

fit.10=lme(var~(Group*Died),
      data=dat34,
      random=list(ID=pdDiag(~Group*Died)),
      corr=NULL)
anova(fit.10, type="marginal", adjustSigma=F)
#               numDF denDF  F-value p-value
#(Intercept)    1   101 93.90211  <.0001
#Group          1    22  1.75311  0.1991
#Died            1    22  6.84379  0.0158
#Group:Died   1    22  3.11458  0.0915

fit.11=lme(var~Group*Died,
      data=dat34,
      random=list(ID=pdDiag(~Group*Died)))
anova(fit.11, type="marginal", adjustSigma=F)
#               numDF denDF  F-value p-value
#(Intercept)    1   101 93.90211  <.0001
#Group          1    22  1.75311  0.1991
#Died            1    22  6.84379  0.0158
#Group:Died   1    22  3.11458  0.0915

fit.12=lme(var~Group*Died,
      data=dat34,
      random=list(ID=pdDiag(~Event_name)))
anova(fit.12, type="marginal", adjustSigma=F)
#               numDF denDF  F-value p-value
#(Intercept)    1   101 87.33040  <.0001
#Group          1    22  1.09661  0.3064
#Died            1    22  5.46329  0.0289
#Group:Died   1    22  2.94589  0.1001

fit.13=lme(var~Group*Died,
      data=dat34,
      random=list(ID=pdDiag(~Group)))
anova(fit.13, type="marginal", adjustSigma=F)
#               numDF denDF  F-value p-value
#(Intercept)    1   101 77.99004  <.0001
#Group          1    22  1.46275  0.2393
#Died            1    22  5.84535  0.0243
#Group:Died   1    22  3.04103  0.0951

fit.14=lme(var~Group*Died,
      data=dat34,
      random=list(ID=pdDiag(~Died)))
anova(fit.14, type="marginal", adjustSigma=F)
#               numDF denDF  F-value p-value
#(Intercept)    1   101 85.83800  <.0001
#Group          1    22  1.60624  0.2183
#Died            1    22  4.71795  0.0409
#Group:Died   1    22  2.65379  0.1175

fit.15=lme(var~Group*Died,
      data=dat34,
      random=~1|ID,
      corr=corCompSymm())
anova(fit.15, type="marginal", adjustSigma=F)
#same as fit.13

fit.16=lme(var~Group*Died,
      data=dat34,
      random=~1|ID/Event_name)
anova(fit.16, type="marginal", adjustSigma=F)
#same as fit.13
#######################
_______________________________________________
R-sig-ecology mailing list
R-sig-ecology-0bNBQ1PAWB4BXFe83j6qeQ< at >public.gmane.org
https://stat.ethz.ch/mailman/listinfo/r-sig-ecology

Re: Continuous (Non-Count) Skewed Data With Many Zeros

Rich Shepard — 2012-05-16T15:01:12


Brian,

   Gotcha'. I'm learning of options that I need to carefully learn and
evaluate. Will install and learn the quantreg package.

Thanks,

Rich