I would like to use wormplots to visually inspect calculated model of birthweight percentiles. xvar is set do D$T that has weeks of gestation (18..42).
When I use wp() function from GAMLSS module, some of wormplots are just empty spaces. No error is shown in console of RStudio. Is that a bug? ( I don't know how to attach these wormplot charts to this post).
Could you please turn this into a self-contained reprex (short for minimal reproducible example)? It will help us help you if we can be sure we're all working with/looking at the same stuff.
If you've never heard of a reprex before, you might want to start by reading the tidyverse.org help page. The reprex dos and don'ts are also useful.
For pointers specific to the community site, check out the reprex FAQ, linked to below.
1 Like
@mara suggested using a reprex. Using a reprex not only makes it easier for us to help you it also makes it easy to show what you code produces, including plots.
Here is an example of a reprex that includes a plot:
suppressPackageStartupMessages(library(tidyverse))
suppressPackageStartupMessages(library(ggplot2))
# table some data
tbl <- tribble(
~x, ~y,
1,2,
2,3,
5,7,
3,1,
)
# shows the table content
tbl
#> # A tibble: 4 x 2
#> x y
#> <dbl> <dbl>
#> 1 1. 2.
#> 2 2. 3.
#> 3 5. 7.
#> 4 3. 1.
plt <- ggplot(tbl, aes(x,y)) + geom_line()
# shows the plot
plt
Created on 2018-03-16 by the reprex package (v0.2.0).
And more info about reprex's
A prose description isn't sufficient, you also need to make a simple reprex that:
- Builds the input data you are using.
- The function you are trying to write, even if it doesn't work.
- Usage of the function you are trying to write, even if it doesn't work.
- Builds the output data you want the function to produce.
You can learn more about reprex's here:
https://www.jessemaegan.com/post/so-you-ve-been-asked-to-make-a-reprex/
Right now the is an issue with the version of reprex that is in CRAN so you should download it directly from github.
Until CRAN catches up with the latest version install reprex with
devtools::install_github("tidyverse/reprex")
The reason we ask for a reprex is that it is the easiest and quickest way for someone to understand the issue you are running into and answer it.
Nearly everyone here who is answering questions is doing it on their own time and really appreciate anything you can do to minimize that time.
2 Likes
I did the reprex but for new users only one image is allowed so I attach only one wrong wormplot as an example.
wp(G_D,xvar=D$T1,n.inter=9,ylim.worm=2) #NOT OK - 7 instead of 9
#> number of missing points from plot= 0 out of 266
#> number of missing points from plot= 0 out of 305
#> number of missing points from plot= 0 out of 261
#> number of missing points from plot= 0 out of 352
#> number of missing points from plot= 0 out of 403
#> number of missing points from plot= 0 out of 216
#> number of missing points from plot= 0 out of 70
and here is the code creating GAMLSS model...
library(readxl)
D <- read_excel("C:/ncbk.xls")
library(gamlss)
#> Ĺadowanie wymaganego pakietu: splines
#> Ĺadowanie wymaganego pakietu: gamlss.data
#> Ĺadowanie wymaganego pakietu: gamlss.dist
#> Ĺadowanie wymaganego pakietu: MASS
#> Ĺadowanie wymaganego pakietu: nlme
#> Ĺadowanie wymaganego pakietu: parallel
#> ********** GAMLSS Version 5.0-6 **********
#> For more on GAMLSS look at http://www.gamlss.org/
#> Type gamlssNews() to see new features/changes/bug fixes.
G_D <- lms(y=D$M1,x=D$T1,k=3,cent=c(5,10,25,50,75,90,95),families="BCTo")
#> *** Initial fit***
#> Warning in is.na(data): 'is.na()' zastosowane do nie-listy lub nie-wektora
#> typu 'NULL'
#> GAMLSS-RS iteration 1: Global Deviance = 27233.87
#> GAMLSS-RS iteration 2: Global Deviance = 27233.87
#> *** Fitting BCTo ***
#> Warning in is.na(data): 'is.na()' zastosowane do nie-listy lub nie-wektora
#> typu 'NULL'
#> GAMLSS-RS iteration 1: Global Deviance = 27069.67
#> GAMLSS-RS iteration 2: Global Deviance = 27070.51
#> GAMLSS-RS iteration 3: Global Deviance = 27071.76
#> GAMLSS-RS iteration 4: Global Deviance = 27072.68
#> GAMLSS-RS iteration 5: Global Deviance = 27073.42
#> GAMLSS-RS iteration 6: Global Deviance = 27073.89
#> GAMLSS-RS iteration 7: Global Deviance = 27074.12
#> GAMLSS-RS iteration 8: Global Deviance = 27074.24
#> GAMLSS-RS iteration 9: Global Deviance = 27074.28
#> GAMLSS-RS iteration 10: Global Deviance = 27074.31
#> GAMLSS-RS iteration 11: Global Deviance = 27074.34
#> GAMLSS-RS iteration 12: Global Deviance = 27074.37
#> GAMLSS-RS iteration 13: Global Deviance = 27074.39
#> GAMLSS-RS iteration 14: Global Deviance = 27074.41
#> GAMLSS-RS iteration 15: Global Deviance = 27074.43
#> GAMLSS-RS iteration 16: Global Deviance = 27074.44
#> GAMLSS-RS iteration 17: Global Deviance = 27074.46
#> GAMLSS-RS iteration 18: Global Deviance = 27074.47
#> GAMLSS-RS iteration 19: Global Deviance = 27074.48
#> GAMLSS-RS iteration 20: Global Deviance = 27074.49
#> Warning in RS(): Algorithm RS has not yet converged
#> % of cases below 4.704867 centile is 5.018687
#> % of cases below 10.33432 centile is 10.03737
#> % of cases below 25.20193 centile is 25.14682
#> % of cases below 49.24801 centile is 50.08009
#> % of cases below 75.64799 centile is 75.17352
#> % of cases below 89.61366 centile is 90.17619
#> % of cases below 95.00755 centile is 94.98131
any ideas what might be wrong with missing wormplots or is it a bug...?
There is no bug you have seven intervals each wormplot corresponds to one of these intervals of your x var. I think the ordering is bottom to top left to right. The plot above your wormplot is called a shingle.
Many thanks, I did not notice it before.
Do you have any idea why my explanatory variable is divided to seven intervals, instead of nine? (n.inter=9)?
regards,
MarekK
can you attach the data making the error the built in example data produce nine intervals
library(gamlss)
#> Loading required package: splines
#> Loading required package: gamlss.data
#>
#> Attaching package: 'gamlss.data'
#> The following object is masked from 'package:datasets':
#>
#> sleep
#> Loading required package: gamlss.dist
#> Loading required package: MASS
#> Loading required package: nlme
#> Loading required package: parallel
#> ********** GAMLSS Version 5.1-3 **********
#> For more on GAMLSS look at http://www.gamlss.org/
#> Type gamlssNews() to see new features/changes/bug fixes.
data(abdom)
m1 <- lms(y,x , data=abdom, cent=c(5,10,25,50,75,90,95),families="BCTo")
#> *** Initial fit***
#> GAMLSS-RS iteration 1: Global Deviance = 4935.853
#> GAMLSS-RS iteration 2: Global Deviance = 4935.853
#> *** Fitting BCTo ***
#> GAMLSS-RS iteration 1: Global Deviance = 4762.119
#> GAMLSS-RS iteration 2: Global Deviance = 4759.821
#> GAMLSS-RS iteration 3: Global Deviance = 4759.8
#> GAMLSS-RS iteration 4: Global Deviance = 4759.814
#> GAMLSS-RS iteration 5: Global Deviance = 4759.823
#> GAMLSS-RS iteration 6: Global Deviance = 4759.825
#> GAMLSS-RS iteration 7: Global Deviance = 4759.828
#> GAMLSS-RS iteration 8: Global Deviance = 4759.828
#> % of cases below 5.711605 centile is 5.081967
#> % of cases below 11.14911 centile is 10
#> % of cases below 24.40083 centile is 25.08197
#> % of cases below 49.40677 centile is 50
#> % of cases below 76.86124 centile is 74.91803
#> % of cases below 89.78305 centile is 90
#> % of cases below 94.56518 centile is 94.91803
wp(m1,xvar=abdom$x,n.inter=9)
#> number of missing points from plot= 0 out of 68
#> number of missing points from plot= 0 out of 71
#> number of missing points from plot= 0 out of 67
#> number of missing points from plot= 0 out of 67
#> number of missing points from plot= 0 out of 66
#> number of missing points from plot= 0 out of 71
#> number of missing points from plot= 0 out of 65
#> number of missing points from plot= 0 out of 69
#> number of missing points from plot= 0 out of 66
Created on 2019-03-26 by the reprex package (v0.2.1)