Re: [R] axis labels not showing
It's very simple my boy!! Do you already to play with mar? So..Try to change the values this object. For example, par(..., mar=c(*5*,2,2,2)). Bye! -- View this message in context: http://r.789695.n4.nabble.com/axis-labels-not-showing-tp4541268p4541354.html Sent from the R help mailing list archive at Nabble.com. __ 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] Help ordinal mixed model!
Good afternoon, gentlemen! After several days studying and researching on
categorical data (various forums with answers from the owner of the library
- all incipient) how to interpret the output the function MCMCglmm, come to
enlist the help of you, if someone has already worked with MCMCglmm function
in the case of variables ordinal dependent. I've read and reread all the
pdf's of the package, the coursenotes Jarrod, finally, I'm exhausted. To
clarify the database, the treatment (called fases) consist of three levels
(1-pre, 2-propolis and 3-vincris) and the ordinal variable response has
three categories (1-normal, 2-agudo, 3 - cronico). See table!
du -
transform(read.table('http://dl.dropbox.com/u/33619290/Dados/Dtest.txt',h=T),FASES=factor(FASES),ALT.RENAIS=ordered(ALT.RENAIS))
summary(du)
library('MCMCglmm')
du - subset(du, ALT.RENAIS != 'NA')
tabela - table(du[,c(2,4)])
tabela
colnames(tabela) - c('Normal','Aguda','Crônica')
rownames(tabela) - c('Pre','Propolis','Vincr')
tabela
#the mixed model:
set.seed(1)
mod1 - MCMCglmm(ALT.RENAIS ~-1+FASES, random= ~ ANIMAIS,
family='ordinal',pl=TRUE,data=du)
summary(mod1)
Then the pain starts, since the documentation is insufficient in this case.
According to him Jarrod (forums), the a posteriori means of the coefficients
of the covariates are the probit scale. According to my study, these
coefficients are the scores of standard normal distribution. More scores
should not correspond to cutpoints? In this case, we would have j (response
variable categories) -1 cutpoints, ie, two cutpoints. The output shows me
only one cutpoint. How can then calculate the probabilities with only one
cutpoint? According to the documentation (Vignettes, page 22), if P (y = k)
= F (yk | l (vlatente), 1) - F (yk-1 | l, 1), this '1' would probably be the
category '1' of the dependent variable? Anyway gentlemen, how can I extract
the probabilities for the stages for each category of the dependent
variable? I thank everyone's attention.
#Absurd results!
latentv - mean(mod1$Liab)
cutpoint - mean(mod1$CP)
pnorm(-(latentv), 0, sqrt(2))
pnorm(cutpoint - (latentv),0, sqrt(2)) - pnorm((latentv),0, sqrt(2))
1- pnorm(cutpoint - (latentv),0, sqrt(2))
#this would have a logical outcome to some extent, more would just be
referring to category '1' of the dependent variable. And the other?
bre -
c(mean(mod1$Liab),mean(mod1$Sol[,1]),mean(mod1$Sol[,2]),mean(mod1$Sol[,3]))
pnorm(bre[2])-pnorm(bre[1])
pnorm(bre[3])-pnorm(bre[2])
pnorm(bre[4])-pnorm(bre[3])# negative probability?
--
View this message in context:
http://r.789695.n4.nabble.com/Help-ordinal-mixed-model-tp4501943p4501943.html
Sent from the R help mailing list archive at Nabble.com.
__
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] To add cut-off points in surface response with lattice
Good morning gentlemen!
I'm not a fan of the lattice due to a large number of procedures what should
be done to reach a simple goal, but have confess that in some cases the
graphics are way better than the graphics. Some days I have been searching
without success as is to add a cut-off point on a graph of response surface.
It is interesting that the researcher to look at the graph spotting cut-off
point. Here is a minimal code reproducible.
require(plotrix)
jet.colors - colorRampPalette( c(blue, green) )
x - seq(-1.95, 1.95, length=30)
y - seq(-1.95, 1.95, length=35)
da - expand.grid(x=x, y=y)
da$z - with(da, x*y^2)
require(lattice)
panel.3d.contour -
function(x, y, z, rot.mat, distance,
nlevels = 20, zlim.scaled, ...)
{
add.line - trellis.par.get(add.line)
panel.3dwire(x, y, z, rot.mat, distance,
zlim.scaled = zlim.scaled, ...)
clines -
contourLines(x, y, matrix(z, nrow = length(x), byrow = TRUE),
nlevels = nlevels)
for (ll in clines) {
m - ltransform3dto3d(rbind(ll$x, ll$y, zlim.scaled[1]),
rot.mat, distance)
panel.lines(m[1,], m[2,], col = add.line$col,
lty = add.line$lty, lwd = add.line$lwd)
}
}
wireframe(z~x+y, da, drape=TRUE,
scales=list(arrows=FALSE),col.regions=jet.colors(100),panel.3d.wireframe=panel.3d.contour)
Tank' s!
--
View this message in context:
http://r.789695.n4.nabble.com/To-add-cut-off-points-in-surface-response-with-lattice-tp3590414p3590414.html
Sent from the R help mailing list archive at Nabble.com.
__
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] The nls2 function automatically prints the object!
Good morning gentlemen! When I use the function nls2, and store it in an object, that object is automatically printed, without the summary or to draw the object. For example. model - nls2 (...) Number of iterations to convergence: ... Achieved convergence tolerance: ... Nonlinear regression model model: ... ~ ... Date: NULL The B k ... ... ... residual sum-of-squares: ... Number of iterations to convergence: ... Achieved convergence tolerance: ... I would not print automatically, I would like to remain stored in the object model, and only if I authorize impressed, calling the model or a summary. Does anyone know how to get around this? -- View this message in context: http://r.789695.n4.nabble.com/The-nls2-function-automatically-prints-the-object-tp3057350p3057350.html Sent from the R help mailing list archive at Nabble.com. __ 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] How using the weights argument in nls2?
Good morning gentlemen! How using a weighted model in nls2? Values with the nls are logical since values with nls2 are not. I believe that this discrepancy is due to I did not include the weights argument in nls2. Here's an example: MOISTURE - c(28.41640, 28.47340, 29.05821, 28.52201, 30.92055, 31.07901, 31.35840, 31.69617, 32.07168, 31.87296, 31.35525, 32.66118, 33.23385, 32.72256, 32.57929, 32.12674, 52.35225, 52.77275, 64.90770, 64.90770, 85.23800, 84.43300, 68.96560, 68.41395, 70.82880, 71.18400, 96.13240, 96.07920, 95.35160, 94.71660, 87.59190, 88.63250, 89.78760, 90.17820, 88.46160, 87.10860, 94.86660, 94.51830, 75.79000, 76.98780, 144.70950, 143.88950, 111.58620, 112.71510, 120.85300, 121.43100, 116.34840, 114.87420, 195.35040, 191.36040, 265.35220, 267.25450, 227.13700, 228.78000, 238.37120, 242.70700, 299.54890, 291.04110, 220.09920, 219.82650, 236.79150, 243.70710, 208.79880, 208.12260, 417.21420, 429.59480, 360.91080, 371.66400, 357.72520, 360.53640, 383.82600, 383.82600, 434.02700, 432.57500, 440.56260, 438.32340, 468.69600, 469.82140, 497.93680, 497.17010) YEARS - rep(c(86,109, 132, 158, 184, 220, 254, 276, 310, 337),c(8,8,8,8,8,8,8,8,8,8)) VARIANCE - c(2.0879048 , 2.0879048,2.0879048,2.0879048, 2.0879048, 2.0879048,2.0879048,2.0879048,0.3442724,0.3442724, 0.3442724,0.3442724,0.3442724,0.3442724, 0.3442724, 0.3442724, 151.9481710, 151.9481710, 151.9481710, 151.9481710, 151.9481710, 151.9481710, 151.9481710, 151.9481710, 115.3208995, 115.3208995, 115.3208995, 115.3208995, 115.3208995, 115.3208995, 115.3208995, 115.3208995, 51.9965027, 51.9965027, 51.9965027, 51.9965027, 51.9965027, 51.9965027, 51.9965027, 51.9965027, 180.0496045, 180.0496045, 180.0496045, 180.0496045, 180.0496045, 180.0496045, 180.0496045, 180.0496045, 791.3223240, 791.3223240, 791.3223240, 791.3223240, 791.3223240, 791.3223240, 791.3223240, 791.3223240, 1280.0179973, 1280.0179973, 1280.0179973, 1280.0179973, 1280.0179973, 1280.0179973, 1280.0179973, 1280.0179973, 728.9582154, 728.9582154, 728.9582154, 728.9582154, 728.9582154, 728.9582154, 728.9582154, 728.9582154, 752.4591144, 752.4591144, 752.4591144, 752.4591144, 752.4591144, 752.4591144, 752.4591144, 752.4591144) test - data.frame(YEARS,MOISTURE,VARIANCE) mod.nls - nls(MOISTURE ~ A/(1+B*exp(-k*YEARS)), data = test, weights = 1/VARIANCE, start = list(A=1500, B=200, k=0.03345), control=list(maxiter = 500), trace=TRUE) summary(mod.nls) Following the example of pdf! st1 - expand.grid(A = seq(0, 2000, len = 10), B = seq(0, 500, len = 10), k = seq(-1, 10, len = 10)) st1 mod.nls2 -nls2(MOISTURE ~ A/(1+B*exp(-k*YEARS)), data = test, start = st1, algorithm=brute-force) mod.nls2 I appreciate everyone's attention. -- View this message in context: http://r.789695.n4.nabble.com/How-using-the-weights-argument-in-nls2-tp2524328p2524328.html Sent from the R help mailing list archive at Nabble.com. __ 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] Res: How to use the function glht of multcomp package to test interaction?
Hi Richard, First thank you for your attention. Actually the way it approached the examples of statements do not like a lot, because the calculations are done separately for each factor of interest to the interaction. Why will not it pleases me so much? Tukey's tests as for example using the mean square error (MSE) for the calculation of minimum significant differences (MSD). When the analysis is done separately, the contrasts are made with the new model, rather than the original model which contained the interactions. Thus the contrasts are different, since they use different MSE. What is your opinion? M.Sc Ivan Bezerra Allaman Zootecnista Doutorando em Produção Animal/Aquicultura - UFLA email e msn - ivanala...@yahoo.com.br Tel: (35)3826-6608/9925-6428 De: Richard M. Heiberger [via R] ml-node+2236373-99253941-109...@n4.nabble.com Enviadas: Domingo, 30 de Maio de 2010 13:09:37 Assunto: Re: How to use the function glht of multcomp package to test interaction? Usually you will want to look at simple effects of the factors when there is interaction. Please look at the WoodEnergy demos in the HH package. These examples use glht. install.packages(HH) ## if you don't currently have HH library(HH) demo(MMC.WoodEnergy-aov) demo(MMC.WoodEnergy) Rich [[alternative HTML version deleted]] __ [hidden email] 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. View message @ http://r.789695.n4.nabble.com/How-to-use-the-function-glht-of-multcomp-package-to-test-interaction-tp2236315p2236373.html To unsubscribe from How to use the function glht of multcomp package to test interaction?, click here. -- View this message in context: http://r.789695.n4.nabble.com/How-to-use-the-function-glht-of-multcomp-package-to-test-interaction-tp2236315p2236706.html Sent from the R help mailing list archive at Nabble.com. [[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] How to use the function glht of multcomp package to test interaction?
It's been a few weeks I'm racking my brains on how to use the function glht the package multcomp to test interactions. Unfortunately, the creator of the package forgot to put a sample in pdf package how to do it. I have looked in several places, but found nothing. If someone for the love of God can help me I'll be extremely grateful. The model is glm. -- View this message in context: http://r.789695.n4.nabble.com/How-to-use-the-function-glht-of-multcomp-package-to-test-interaction-tp2236315p2236315.html Sent from the R help mailing list archive at Nabble.com. __ 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] Split-plot design in GLM with only fixed factors.
Good evening gentlemen! I have a test in split-plot with randomized block design where my answer is a binomial variable. I wonder if there is any way I can calculate the probability of my factors considering the design errors in the case are two. I looked at various threads here and elsewhere, and unfortunately no to answer objective my problem that is very simple. My interest isn't to estimate variance components, so I see no reason to use functions like LM as found on most topics. If anyone knows a way to calculate the probability of my factors as we do in an LM model, ie, announcing the types of errors so that the probability factor in interest is not prejudiced, I'd be grateful. -- View this message in context: http://r.789695.n4.nabble.com/Split-plot-design-in-GLM-with-only-fixed-factors-tp2228126p2228126.html Sent from the R help mailing list archive at Nabble.com. __ 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] Res: Using the zero-inflated binomial in experimental designs
Hi Ben! Following his recommendations I did the following: 1st step: I compared the best model for binomial and binomial inflates. 1.1 Best model for Binomial. dg$resp.mumi - cbind(dg$MUMI,dg$NT - dg$MUMI) dg names(dg) mod.mumi.binomial - glm(resp.mumi ~ factor(PARTO)*REG, family=binomial, data = dg) summary(mod.mumi.binomial) mod.mumi.binomial1 - glm(resp.mumi ~ factor(PARTO) + REG, family=binomial, data = dg) summary(mod.mumi.binomial1) mod.mumi.binomial2 - glm(resp.mumi ~ REG, family=binomial, data = dg) summary(mod.mumi.binomial2) pchisq(-2*(logLik(mod.mumi.binomial1)-logLik(mod.mumi.binomial)),lower.tail=FALSE, df= 15) [1] 0.1354171 pchisq(-2*(logLik(mod.mumi.binomial2)-logLik(mod.mumi.binomial1)),lower.tail=FALSE, df= 5) [1] 0.06030012 The 5% significance level, we can choose the most parsimonious model, ie the model mod.mumi.binomial2. 1.2 Best model for binomial inflates library(VGAM) mod.mumi.binomialinflacionada - vglm(resp.mumi ~ factor(PARTO)*REG,zibinomial, data = dg) summary(mod.mumi.binomialinflacionada) mod.mumi.binomialinflacionada1 - vglm(resp.mumi ~ factor(PARTO)+REG,zibinomial, data = dg) summary(mod.mumi.binomialinflacionada1) mod.mumi.binomialinflacionada2 - vglm(resp.mumi ~ REG,zibinomial, data = dg) summary(mod.mumi.binomialinflacionada2) pchisq(-2*logLik(mod.mumi.binomialinflacionada1)-logLik(mod.mumi.binomialinflacionada)),lower.tail=FALSE, df= 15) [1] 0.1477837 pchisq(-2*logLik(mod.mumi.binomialinflacionada2)-logLik(mod.mumi.binomialinflacionada1)),lower.tail=FALSE, df= 5) [1] 0.0989934 The 5% significance level, we can choose the most parsimonious model, ie the model mod.mumi.binomialinflacionada2. 2st step: Compare the best model of the binomial model with the best of inflated binomial. pchisq(-2*(logLik(mod.mumi.binomial2)-logLik(mod.mumi.binomialinflacionada2)),lower.tail=FALSE, df= 1) [1] 0.1929690 There wasn't difference between the models. Must i choose the most parsimonious model? Thanks for your attention and sorry for the inconvenience. -- View this message in context: http://r.789695.n4.nabble.com/Using-the-zero-inflated-binomial-in-experimental-designs-tp2221254p2223819.html Sent from the R help mailing list archive at Nabble.com. __ 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] Using the zero-inflated binomial in experimental designs
I'm trying to use the inflated binomial distribution of zeros (since 75% of the values are zeros) in a randomized block experiment with four quantitative treatments (0, 0.5, 1, 1.5), but I'm finding it difficult, since the examples available in VGAM packages like for example, leave us unsure of how it should be the data.frame for such analysis. Unfortunately the function glm does not have an option to place a family of this kind I'm about, because if I had, it would be easy, made that my goal is simple, just wanting to compare the treatments. For that you have an idea, here is an example of my database. BLOCK NIVNT MUMI Inicial 0 18 0 Inicial 0 15 0 Inicial 0.5 9 0 Inicial 0.5 19 1 Inicial 1 13 1 Inicial 1 11 0 Inicial 1.5 12 2 Inicial 1.5 10 1 Meio 0 13 0 Meio 0 10 2 Meio 0.5 17 0 Meio 0.5 14 1 Meio 1 13 0 Meio 1 9 0 Meio 1.5 110 Meio 1.5 12 1 where: NIV are the treatments; NT is the total number of piglets born; Mumi is the number of mummified piglets NT. Mumi The variable is of interest. If someone can tell me some stuff on how I can do these tests in R, similar to what I would do using the function glm, I'd be grateful. I thank everyone's attention. -- View this message in context: http://r.789695.n4.nabble.com/Using-the-zero-inflated-binomial-in-experimental-designs-tp2221254p2221254.html Sent from the R help mailing list archive at Nabble.com. __ 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] Res: Using the zero-inflated binomial in experimental designs
Hi Ben! First I thank you for your attention. Unfortunately, the ANOVA does not work with vglm. In another email, Rafael warned me that actually a lot of zeros does not necessarily imply a distribution of zeros binomail inflated. So how could I test if my variable is or not a binomial zero inflated? Thanks. M.Sc Ivan Bezerra Allaman Zootecnista Doutorando em Produção Animal/Aquicultura - UFLA email e msn - ivanala...@yahoo.com.br Tel: (35)3826-6608/9925-6428 De: Ben Bolker [via R] ml-node+2221578-2039137178-109...@n4.nabble.com Enviadas: Terça-feira, 18 de Maio de 2010 13:34:01 Assunto: Re: Using the zero-inflated binomial in experimental designs Ivan Allaman ivanalaman at yahoo.com.br writes: I'm trying to use the inflated binomial distribution of zeros (since 75% of the values are zeros) in a randomized block experiment with four quantitative treatments (0, 0.5, 1, 1.5), but I'm finding it difficult, since the examples available in VGAM packages like for example, leave us unsure of how it should be the data.frame for such analysis. Unfortunately the function glm does not have an option to place a family of this kind I'm about, because if I had, it would be easy, made that my goal is simple, just wanting to compare the treatments. For that you have an idea, here is an example of my database. BLOCK NIVNT MUMI Inicial 0 18 0 [snip] where: NIV are the treatments; NT is the total number of piglets born; Mumi is the number of mummified piglets NT. Mumi The variable is of interest. If someone can tell me some stuff on how I can do these tests in R, similar to what I would do using the function glm, I'd be grateful. I thank everyone's attention. something like comparing the likelihoods of m1 - vglm(cbind(MUMI,NT-MUMI)~NIV*BLOCK,zibinomial,data=mydata) m2 - vglm(cbind(MUMI,NT-MUMI)~NIV+BLOCK,zibinomial,data=mydata) m3 - vglm(cbind(MUMI,NT-MUMI)~BLOCK,zibinomial,data=mydata) I don't know whether the anova() method works for VGLM objects or not. By the way, 75% zeroes doesn't necessarily imply zero-inflation -- perhaps it just means a low incidence? __ [hidden email] 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. View message @ http://r.789695.n4.nabble.com/Using-the-zero-inflated-binomial-in-experimental-designs-tp2221254p2221578.html To unsubscribe from Using the zero-inflated binomial in experimental designs, click here. -- View this message in context: http://r.789695.n4.nabble.com/Using-the-zero-inflated-binomial-in-experimental-designs-tp2221254p2221635.html Sent from the R help mailing list archive at Nabble.com. [[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.
