Re: specifying normal distribution for terminal node in aster
Ruth Shaw
Professor, Dept of Ecology, Evolution and Behavior
St. Paul MN 55108
Hi David and Damian,
The short answer to what I think you are looking for (specifying normal distribution for terminal node in aster) is this:
famlist <- list(fam.bernoulli(),fam.normal.location(standard deviation))
Below is a my longer answer that is formulated around some example script. Please let me know if this works for you. If not perhaps if I can see your script I might be able to get at it more quickly.
### Set up aster model
# first I omitted the zeros for the dead trees in my dataframe, “aster_data”. HT12 was the column in my data frame for height
ZHT12<- subset(aster_data,HT12 > 0,)
# next I visually checked the distribution
hist(ZHT12$HT12)
# then I used the script below to get the standard distribution for height
(ZHT12.param <- fitdistr(ZHT12$HT12, "normal"))
# In this case the standard distribution was 518.86003.
# I plugged that into the following aster script using the optional argument “famlist”
# note the first node was Bernoulli for survival
# the second position indicates to aster that the second node is normal using “fam.normal.location”
famlist <- list(fam.bernoulli(),fam.normal.location(518.86003))
# Following this example through a bit more below when the aster graph is specified aster will use Bernoulli which is assigned with the number 1 in fam
# argument for pred 0 through 5 which is survival data. The number 1 used in the argument fam corresponds with the position in the
# order listed above in famlist for the fam argument. The above specified normal distribution assigned using the number 2 in the fam argument will be
# used for pred 6 which is the height data. The number 2 comes from the position in the order listed above in famlist.
pred <- c(0, 1, 2, 3, 4, 5, 6)
fam <- c(1, 1, 1, 1, 1, 1, 2)
### in your call be sure to include (famlist = famlist) for example see the call below
aster.formula(formula = resp ~ varb + BLK, pred = pred, fam = fam,
varvar = varb, idvar = id, root = root, data = chamae2, famlist = famlist)
Best Wishes,
Marcus
Marcus Warwell, Ph.D.
GeneticistForest Service
Rocky Mountain Research Station, FWE program
1221 S. Main Street
Moscow, ID 83843
www.fs.fed.us
Caring for the land and serving people
From: Ruth Shaw [mailto:shaw...@umn.edu]
Sent: Tuesday, April 03, 2018 8:44 AM
To: David Lowry <davidbry...@gmail.com>
Cc: Damian Popovic (popovi...@gmail.com) <popovi...@gmail.com>; Warwell, Marcus V -FS <mwar...@fs.fed.us>
Subject: Re: Aster
Hi David (and Damian),
It's good to hear from you, and of course I'm very glad to know you're using aster. I have not myself used the capability of specifying the terminal node as normal, but my former student, Marcus Warwell, has. I'm taking the liberty of introducing you to each other. Marcus, can you help Damian past this sticking point?
Best to you all!
Ruth
Ruth G. Shaw
Professor, Dept of Ecology, Evolution and BehaviorResident Fellow, Minnesota Center for Philosophy of Science
612 624 7206
Mail to:1479 Gortner Ave., 140 Gortner Laboratory
University of Minnesota
St. Paul MN 55108
On Tue, Apr 3, 2018 at 8:37 AM, David Lowry <davidbry...@gmail.com> wrote:
Dear Ruth,
How are you? I hope all is well and that Spring comes sometime soon.
My graduate student, Damian, is trying to use Aster to analyze his field reciprocal transplant experiment. He has two variables that he is focused on: Survival to the end of the season and Final Dry Weight. We both thought that it would make the most sense to treat dry weight as normally distributed. However, we cannot figure out the proper code to enter it into Aster that way (see below). Could you possibly tell us how to do that properly? I think once we get past this step we should be able to analyze the data.
Thanks a ton!
David
---------- Forwarded message ----------
From: Damian Popovic <popovi...@gmail.com>
Date: Tue, Apr 3, 2018 at 9:27 AM
Subject: Re: Aster
To: David Lowry <davidbry...@gmail.com>
Morning, David!
So that was my initial inkling. When I previously tried the above, I received the following error:
fam <- c(fam.bernoulli(), fam.normal.location(sd))
Error: is.numeric(sd) is not TRUE
After which, I decided to plug in the actual standard deviation of my mass data like so:
sd.mass <- sd(Mass)
fam <- c(fam.bernoulli(), fam.normal.location(sd.mass))
Once I plugged that into my first model, I received the following error:
aout1 <- aster(resp ~ varb + Garden, pred, fam, varb, id, root, data=redata)
Error: length(fam) == ncol(x) is not TRUE
Upon viewing fam, I noticed that the standard deviation was listed as a third column or component to the object. Thus, my best guess is that this method makes the length of fam incongruent with pred - both of which are supposed to have the same number of components (here: 2). When simply using fam.normal.location without specifying standard deviation, I only have the two listed components (bernoulli & normal.location) yet the standard deviation is still somehow included in the object (without bloating the list with an extra component).
Cheers!
Damian
On Tue, Apr 3, 2018 at 6:45 AM, David Lowry <davidbry...@gmail.com> wrote:
I'll see you this morning! What happens if you try the following:
fam <- c(fam.bernoulli(), fam.normal.location(sd))with the (sd) added?
D
On Mon, Apr 2, 2018 at 10:01 PM, Damian Popovic <popovi...@gmail.com> wrote:
Evening, David!
So I've been wracking my brain to solve the "family" conundrum we had earlier. The first pdf lists the families (or response models) known to the package. They are as follows:
fam.bernoulli()
fam.poisson()
fam.truncated.poisson(truncation)
fam.negative.binomial(size)
fam.truncated.negative.binomial(size, truncation)
fam.normal.location(sd)
fam.default()
famfun(fam, deriv, theta)
The numeric identifiers (1-4) used in your code and all of Ruth's examples are a limited list of families hard coded into fam.default(). She has this to say about it:
"For all but fam.default, a list of class "astfam" giving name and values of any hyperparameters.
For fam.default, a list each element of which is of class "astfam". The list of families which
were hard coded in earlier versions of the package."
and
"fam specifies the one-parameter exponential families for the nodes: fam[j] gives the index of the family for node j, the index being for the list of family specifications that is supplied to model fitting functions, the default being the list returned by fam.default()."
Since gaussian functions aren't represented in fam.default(), I decided to write out the code in longhand. Thus I ended up with fam <- c(fam.bernoulli(), fam.normal.location). This worked until I attempted to fit my first model - the code for the model being aout1 <- aster(resp ~ varb + Garden, pred, fam, varb, id, root, data=redata). I received the following error: "Error in stopifnot(all(fam == as.integer(fam))) : (list) object cannot be coerced to type 'integer'". I suppose the problem is that upon coding fam in longhand, R programs the object as a list rather than integers as I assume it would if I used the fam.default() identifiers. I hope this isn't too confusing; regardless, I'm stumped. Any ideas? I'll be in lab starting ~10:15am tomorrow if you want to chat in person for further clarification.
Cheers!
Damian
On Mon, Apr 2, 2018 at 3:25 PM, David Lowry <davidbry...@gmail.com> wrote:
I can also e-mail Ruth Shaw. I know her decently well.
D
--David B. Lowry
Assistant Professor
Plant Biology Department
Michigan State University
Plant Biology Laboratories Room 268
517-432-4882
http://davidbryantlowry.wordpress.com/
--
Damian Popovic
Graduate Student
Department of Plant Biology (PLB)
Michigan State University
Plant Biology Laboratories, Room 268
--David B. Lowry
Assistant Professor
Plant Biology Department
Michigan State University
Plant Biology Laboratories Room 268
517-432-4882
http://davidbryantlowry.wordpress.com/
--
Damian Popovic
Graduate Student
Department of Plant Biology (PLB)
Michigan State University
Plant Biology Laboratories, Room 268
--David B. Lowry
Assistant Professor
Plant Biology Department
Michigan State University
Plant Biology Laboratories Room 268
517-432-4882
http://davidbryantlowry.wordpress.com/
This electronic message contains information generated by the USDA solely for the intended recipients. Any unauthorized interception of this message or the use or disclosure of the information it contains may violate the law and subject the violator to civil or criminal penalties. If you believe you have received this message in error, please notify the sender and delete the email immediately.
Ruth Shaw
Professor, Dept of Ecology, Evolution and Behavior
St. Paul MN 55108
Hi Damain,
I do not have a good answer for either of your questions. I adopted the use of “fitdistr” to obtain the sd following the direction of Charlie Geyer and didn’t think to check for differences or look at alternative means to calculate sd. I do not know what the difference is between sd calculated using “fitdistr” in MASS and the traditional sd(). The Mass manual doesn’t present a formula. I guess you probably already looked there too. Sorry I can’t help more, but please do let me know if you find a good answer.
Best of luck.
Marcus
From: Damian Popovic [mailto:popovicdamian@gmail.com]
Sent: Tuesday, April 03, 2018 11:44 AM
To: Ruth Shaw <shaw...@umn.edu>
Cc: Warwell, Marcus V -FS <mwar...@fs.fed.us>; David Lowry <davidbry...@gmail.com>; aster-analysis-user-group@googlegroups.com
Subject: Re: specifying normal distribution for terminal node in aster
This is fantastic! Again, thank you for your time & aid Marcus! A couple quick questions before I continue hammering out this code. Is it necessary to calculate the normal distribution using the "fitdistr" function from the MASS package or is it entirely kosher to plug in the calculated standard deviation of the nonzero masses directly into the famlist code? It appears that the standard deviation computed using the traditional sd() function comes out slightly different to that which is estimated by the fitdistr function. What causes this disparity and which should I use?
Cheers!
Damian