Forum | VSN International
http://www.vsni.co.uk/forum/
The forum for discussion of VSN International software and general statistical issuesenCopyright 2000-2010 VSN International Ltd.support@vsni.co.uksupport@vsni.co.ukhttp://blogs.law.harvard.edu/tech/rss60Sat, 25 Mar 2017 03:54:53 GMTSat, 25 Mar 2017 03:54:53 GMThttp://www.vsni.co.uk/common/images/asreml145.gifForum | VSN International
http://www.vsni.co.uk/forum/
ASReml - glmm with asreml_R
http://www.vsni.co.uk/forum/viewtopic.php?p=5042#5042
Dear all,
<br />
We are currently trying to fit GLMMs using asreml with spatial autocorrelation. As a preliminary step, we have been comparing the results from asreml and lmer, without the spatial component, using simulated data. We used the function sim.glmm published by Johnson et al. to simulate binomial data and compared two simple models fit with glmer and asreml.
<br />
In our comparisons, the reported effects were different for GLMER and ASReml. This is probably because we did not fit the same model. As you can see in the following example, we have been using asreml with the proportion data, and specified the number of trial with the weight option (lmer uses the matrix of success,failures). Specifying weights did not affect the likelihood but changed the residual variance (the overdispersion effect). In this simulated data, n is constant, but in our real dataset it changes, and specifying weights may be important. Is the weight option the correct way to specify this information in asreml?
<br />
We would be extremely grateful if you could provide us with some help and pieces of advice.
<br />
<br />
Example
<br />
<br />
<div class="code">Code:<div class="inside_code">library(RCurl)
<br />
options(RCurlOptions=list(cainfo=system.file("CurlSSL","cacert.pem",package="RCurl")))
<br />
eval(expr=parse(text=getURL("https://raw.githubusercontent.com/pcdjohnson/sim.glmm/master/sim.glmm.R")))
<br />
<br />
# simulate data
<br />
latsq <-
<br />
rbind(
<br />
c("C", "E1", "E2", "E3", "E4", "E5"),
<br />
c("E5", "C", "E1", "E2", "E3", "E4"),
<br />
c("E4", "E5", "C", "E1", "E2", "E3"),
<br />
c("E3", "E4", "E5", "C", "E1", "E2"),
<br />
c("E2", "E3", "E4", "E5", "C", "E1"),
<br />
c("E1", "E2", "E3", "E4", "E5", "C"))
<br />
colnames(latsq) <- paste("hut", 1:nrow(latsq), sep = "")
<br />
rownames(latsq) <- paste("week", 1:ncol(latsq), sep = "")
<br />
latsq
<br />
<br />
mosdata <-
<br />
expand.grid(
<br />
hut = factor(1:ncol(latsq)),
<br />
week = factor(1:nrow(latsq)),
<br />
night = factor(1:6))
<br />
mosdata <- mosdata[order(mosdata$hut, mosdata$week, mosdata$night),]
<br />
mosdata$observation <- factor(formatC(1:nrow(mosdata), flag="0", width=3))
<br />
mosdata$net <- factor(diag(latsq[mosdata$week, mosdata$hut]))
<br />
mosdata$n <- 25
<br />
set.seed(321123) # set seed for random number generator to give repeatable results
<br />
mosdata <-
<br />
sim.glmm(design.data = mosdata,
<br />
fixed.eff = list(
<br />
intercept = qlogis(0.7),
<br />
net = log(c(C = 1, E1 = 1.7, E2 = 1.7, E3 = 1.7, E4 = 1.7, E5 = 1.7))),
<br />
rand.V = c(hut = 0.5, week = 0.5, observation = 1),
<br />
distribution = "binomial")
<br />
mosdata
<br />
<br />
# test with lme4
<br />
library(lme4)
<br />
fit <-
<br />
glmer(
<br />
cbind(response, n - response) ~
<br />
net + (1 | hut) + (1 | week) + (1 | observation),
<br />
family = binomial, data = mosdata)
<br />
summary(fit)$AICtab
<br />
as.data.frame(VarCorr(fit))
<br />
<br />
# test with asreml
<br />
library(asreml)
<br />
mosdata$prop<-mosdata$response/mosdata$n
<br />
m1 <- asreml(prop ~ net, random = ~ hut + week , data = mosdata,
<br />
family=asreml.binomial(link='logit',dispersion=NA),na.method.X='omit',weight=n,trace=F)
<br />
m1$loglik
<br />
summary(m1)$varcomp
<br />
<br />
m2 <- asreml(prop ~ net, random = ~ hut + week , data = mosdata,
<br />
family=asreml.binomial(link='logit',dispersion=NA),na.method.X='omit',trace=F)
<br />
m2$loglik
<br />
summary(m2)$varcomp
<br />
</div></div>ASRemlhttp://www.vsni.co.uk/forum/posting.php?mode=reply&t=1701gaylThu, 23 Mar 2017 14:28:20 GMThttp://www.vsni.co.uk/forum/viewtopic.php?p=5042#5042ASReml - RE: Predicting from binomial mixed effects model
http://www.vsni.co.uk/forum/viewtopic.php?p=5041#5041
<div class="quote">kiwi wrote:<div class="inside_quote">If I get predictions as so:
<br />
<div class="code">Code:<div class="inside_code">
<br />
<br />
useLogCTs <- seq(min(xx$logCT), max(xx$logCT), len = 5)
<br />
preds <- predict(xx.asr2, classify = "Species:Stage:Temperature:logCT",
<br />
levels = list(logCT = useLogCTs),
<br />
type = "response",
<br />
weights = "Total")
<br />
<br />
pred.df <- as.data.frame(preds$predictions$pvals)
<br />
<br />
> summary(pred.df$transformed)
<br />
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
<br />
0.000017 1.000000 1.000000 0.827300 1.000000 1.000000 10
<br />
</div></div>
<br />
<br />
Why am I getting only extreme predictions? Looking at the
<br />
coefficients in xx.asr2, all appear to be about 5 or 6 orders of
<br />
magnitude too high. Is that to be expected?
<br />
<br />
<br />
PS: I tried to show how I got xx.asr2 but it would vanish.
<br />
<br />
TIA</div></div>
<br />
Though I had hoped to get feedback on the syntax used producing my xx.asr2 object, the software seems to be unable to cope with that amount of text, though it's scarcely in kilobytes. It would be sufficient to learn if anyone had succeeded in getting useful predictions from a binomial model. I tried scaling down the coefficients by a number of orders of magnitude, but apparently, the coefficients stored in the asreml object aren't the ones used in the prediction process. One thing that is slightly reassuring is that the predictions change from 0 to 1 at believable points so something is working.ASRemlhttp://www.vsni.co.uk/forum/posting.php?mode=reply&t=1700kiwiWed, 22 Mar 2017 20:41:57 GMThttp://www.vsni.co.uk/forum/viewtopic.php?p=5041#5041ASReml - RE: New licensing process for ASReml4?
http://www.vsni.co.uk/forum/viewtopic.php?p=5040#5040
<div class="quote">kiwi wrote:<div class="inside_quote">I can use ASReml-R version 3 without any problem, but when I try to load asreml4, there's a pop-up screen complaining about a missing or invalid license file.
<br />
<br />
It seems to expect it to be in the /home/ directory instead of where it's indicated in an environment variable. So I made a copy there but the complaint persists.
<br />
<br />
Where is the information describing the change?
<br />
<br />
Thanks</div></div>
<br />
Just a follow up since that post had quite a few views, but no answers.
<br />
I eventually found out from VSNi support that the "reason" was that ASReml4 for R has not been released.
<br />
<br />
Nobody should get the impression that there's something wrong with their asreml.lic file.
<br />
<br />
KiwiASRemlhttp://www.vsni.co.uk/forum/posting.php?mode=reply&t=1698kiwiMon, 20 Mar 2017 20:22:14 GMThttp://www.vsni.co.uk/forum/viewtopic.php?p=5040#5040ASReml - Predicting from binomial mixed effects model
http://www.vsni.co.uk/forum/viewtopic.php?p=5039#5039
If I get predictions as so:
<br />
<div class="code">Code:<div class="inside_code">
<br />
<br />
useLogCTs <- seq(min(xx$logCT), max(xx$logCT), len = 5)
<br />
preds <- predict(xx.asr2, classify = "Species:Stage:Temperature:logCT",
<br />
levels = list(logCT = useLogCTs),
<br />
type = "response",
<br />
weights = "Total")
<br />
<br />
pred.df <- as.data.frame(preds$predictions$pvals)
<br />
<br />
> summary(pred.df$transformed)
<br />
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
<br />
0.000017 1.000000 1.000000 0.827300 1.000000 1.000000 10
<br />
</div></div>
<br />
<br />
Why am I getting only extreme predictions? Looking at the
<br />
coefficients in xx.asr2, all appear to be about 5 or 6 orders of
<br />
magnitude too high. Is that to be expected?
<br />
<br />
<br />
PS: I tried to show how I got xx.asr2 but it would vanish.
<br />
<br />
TIAASRemlhttp://www.vsni.co.uk/forum/posting.php?mode=reply&t=1700kiwiSun, 19 Mar 2017 20:45:09 GMThttp://www.vsni.co.uk/forum/viewtopic.php?p=5039#5039ASReml - RE: Correlations in bivariate analyses
http://www.vsni.co.uk/forum/viewtopic.php?p=5038#5038
Thank you for your help Arthur. That worked for me.ASRemlhttp://www.vsni.co.uk/forum/posting.php?mode=reply&t=1699god4rettaTue, 07 Mar 2017 20:42:36 GMThttp://www.vsni.co.uk/forum/viewtopic.php?p=5038#5038ASReml - RE: Correlations in bivariate analyses
http://www.vsni.co.uk/forum/viewtopic.php?p=5037#5037
Dear god4retta,
<br />
<br />
Change the model line to
<br />
<br />
PCV WWT - Trait Tr.damage Tr.brtype !r !{ Tr.id at(Tr,2).dam !} at(Tr,2).ide(dam) Tr.litter
<br />
<br />
and the G structure lines for Tr.id to
<br />
<br />
Tr.id 2
<br />
3 0 US !GP
<br />
.1 .01 .1 .001 .001 .05
<br />
id 0 AINV
<br />
<br />
Where !{ !} make sure that the enclosed model terms are not reorder so we can
<br />
define a 3 cross N variance structure for the 3N effects (N direct PCV effects, N direct WWT effects, N maternal WWT effects)
<br />
<br />
and we define a variance structure for Tr.id which is 3 cross N and so also includes the at(Tr,2).dam effects.
<br />
<br />
The only issue then is that you probably need to provide some plausible starting values.
<br />
The values I have given are plausible if the variance of both traits is 1.
<br />
If the variances were say 100 and 25 (SD 10 and 5), you should change the start values to
<br />
10 .5 2.5 .05 .025 1.25
<br />
<br />
calculated as
<br />
10*10 * 0.1 is 10
<br />
10*5 *.01 is 0.5
<br />
5*5*0.1 is 2.5
<br />
10*5*.001 is 0.05
<br />
5*5*.001 is 0.025
<br />
5*5*.05 is 1.25
<br />
<br />
<br />
In ASReml4 functional suntax, the model would be written as
<br />
<br />
PCV WWT - Trait Tr.damage Tr.brtype !r !{ Tr.id at(Tr,2).dam us(3).nrm(id) !} at(Tr,2).ide(dam) us(Tr).litter
<br />
residual units.us(Tr)
<br />
<br />
and you should not need to worry about start values.ASRemlhttp://www.vsni.co.uk/forum/posting.php?mode=reply&t=1699ArthurFri, 03 Mar 2017 21:42:34 GMThttp://www.vsni.co.uk/forum/viewtopic.php?p=5037#5037ASReml - Correlations in bivariate analyses
http://www.vsni.co.uk/forum/viewtopic.php?p=5036#5036
Hello forum, please can someone help me on how to get genetic correlation between direct effects on one trait and maternal genetic effect on another trait, in a bivariate analysis.
<br />
<br />
For example, in a bivariate analysis of PCV and weaning weight,
<br />
<br />
PCV WWT - Trait Tr.damage Tr.brtype !r Tr.id at(Tr,2).dam at(Tr,2).ide(dam) Tr.litter
<br />
1 2 2
<br />
0
<br />
Trait 0 US
<br />
3*0
<br />
Tr.id 2
<br />
Trait 0 US !GP
<br />
3*0
<br />
id
<br />
Tr.litter 2
<br />
Trait 0 US !GP
<br />
3*0
<br />
Litter
<br />
<br />
How do I get the correlation between direct effect on PCV and maternal effect on weaning weight?ASRemlhttp://www.vsni.co.uk/forum/posting.php?mode=reply&t=1699god4rettaThu, 02 Mar 2017 15:25:20 GMThttp://www.vsni.co.uk/forum/viewtopic.php?p=5036#5036ASReml - New licensing process for ASReml4?
http://www.vsni.co.uk/forum/viewtopic.php?p=5035#5035
I can use ASReml-R version 3 without any problem, but when I try to load asreml4, there's a pop-up screen complaining about a missing or invalid license file.
<br />
<br />
It seems to expect it to be in the /home/ directory instead of where it's indicated in an environment variable. So I made a copy there but the complaint persists.
<br />
<br />
Where is the information describing the change?
<br />
<br />
ThanksASRemlhttp://www.vsni.co.uk/forum/posting.php?mode=reply&t=1698kiwiThu, 02 Mar 2017 01:14:23 GMThttp://www.vsni.co.uk/forum/viewtopic.php?p=5035#5035ASReml - RE: Heteregeneous AR1 x Ar1 model for MET
http://www.vsni.co.uk/forum/viewtopic.php?p=5033#5033
Dear Subash,
<br />
<br />
I typical initial model for a variety trial MET which is plausible across 50 trials
<br />
with common genotypes (on a 64bit PC with say 4 or more Gb) is
<br />
using functional syntax
<br />
<br />
yield ~ mu Trial mv !r xfa1(Trial).Genotype
<br />
residual at(Trial).ar1(row).ar1(column)ASRemlhttp://www.vsni.co.uk/forum/posting.php?mode=reply&t=1696ArthurMon, 20 Feb 2017 00:08:20 GMThttp://www.vsni.co.uk/forum/viewtopic.php?p=5033#5033ASReml - Heteregeneous AR1 x Ar1 model for MET
http://www.vsni.co.uk/forum/viewtopic.php?p=5032#5032
Dear All,
<br />
<br />
I am trying to fit a heteregeneous AR1xAR1 structure to a data from MET where there are 50+ Fields or Environments. For a single field one can use
<br />
<br />
row row AR1 0.1
<br />
col col ARl 0.1
<br />
<br />
Wondering what is the correct syntax to specify this structure for Multi Environments and are there any concerns in terms of memory requirements.
<br />
<br />
<br />
Thanks.
<br />
Best regards,
<br />
-SubashASRemlhttp://www.vsni.co.uk/forum/posting.php?mode=reply&t=1696skrishFri, 17 Feb 2017 18:41:25 GMThttp://www.vsni.co.uk/forum/viewtopic.php?p=5032#5032