Forum | VSN International

GenStat - RE: splines?

hedderleyd — Fri, 24 May 2013 03:31:35 GMT

Carole
I've found Verbyla, Cullis, Kenward and Welham (Applied Statistics, 1999 'The analysis of designed experiments and longitudinal data by using smoothing splines') a good guide, with several examples and a few recipes (!). One catch with splines is that you need to include a slope term as well as a spline term, to catch any overall changes in level with time
Hope the reference helps
Duncan

The contents of this e-mail are confidential and may be subject to legal privilege.
If you are not the intended recipient you must not use, disseminate, distribute or
reproduce all or any part of this e-mail or attachments. If you have received this
e-mail in error, please notify the sender and delete all material pertaining to this
e-mail. Any opinion or views expressed in this e-mail are those of the individual
sender and may not represent those of The New Zealand Institute for Plant and
Food Research Limited.

Post generated using Mail2Forum (http://www.mail2forum.com)

GenStat - splines?

Anonymous — Fri, 24 May 2013 02:55:59 GMT

I have data from a RCBD (8 tmts by 3 reps). Measurements were taken at 6 time points. The response over time is not linear and doesn't suit fitting any standard curves. There is one rep which drops at one time point for all treatments (these plants were apparently more advanced in development and the growth stage may have caused this effect).

What is the best way to test for different treatment responses over time?

I've looked at fitting splines in REML but am not confident with choosing what to put where. I'm thinking this might be the best approach??

I've also used the pol function in a regression model (and ANOVA) but I think I need something less constraining.

Any help would be greatly appreciated.
Thank you
Carole

Post generated using Mail2Forum (http://www.mail2forum.com)

GenStat - RE: Death codes

slangton — Wed, 22 May 2013 09:22:03 GMT

Here's a variant that uses factors for the codes so that A00=100, A01=101, etc. It also uses the Biometris procedure FPOINTER to set up scalars with the appropriate names and values, thus allowing the text form to be used in the RESTRICT statement.

Steve

text codes;val=!t(C00,C01,D22,C02,C03,C04,D01,D07)

vari [val=(100...799)] vcode
txcon [tcode] !t(100(A,B,C,D,E,F,G)),vcode;dec=0;sk=0,1
fpoint tcode;scpo "biometris procedure"
scal scpo[];#vcode

fact [lab=tcode;lev=vcode] fcodes;val=codes

print codes,fcodes
rest fcodes;fcodes.ge.C03.and.fcodes.le.D10
print fcodes
rest fcodes

*******************************************************************
Steve Langton
Defra Agricultural Change & Environment Observatory
Foss House, Kings Pool Tel: 01904 456410
1-2 Peasholme Green GTN: 5137
York YO1 7PX
E-mail: steve.d.langton@defra.gsi.gov.uk
*******************************************************************

Department for Environment, Food and Rural Affairs (Defra)

This email and any attachments is intended for the named recipient only. If you have received it in error you have no authority to use, disclose,
store or copy any of its contents and you should destroy it and inform the sender.
Whilst this email and associated attachments will have been checked for known viruses whilst within Defra systems we can accept no responsibility once it has left our systems.
Communications on Defra's computer systems may be monitored and/or recorded to secure the effective operation of the system and for other lawful purposes.

Post generated using Mail2Forum (http://www.mail2forum.com)

GenStat - RE: Death codes

david@vsn.co.nz — Wed, 22 May 2013 00:44:46 GMT

Hi Brian,

Quote:

Setting up analyses of cause-specific mortality in GenStat used to be simple
when ICD codes for causes were numeric, even for wide ranges of codes,
because GenStat has good ways of specifying ranges of values, which can be
used with 'in' to restrict to selected groups of causes.

However, now that we're in the era of ICD10 codes, which start with a letter,
things are not so straightforward. Suppose I want to analyse for the group of
all neoplasms, i.e. all causes in the range C00-D49. Can anyone suggest a way
of restricting deaths explicitly to this range easily, please? Assume I've created
a factor that contains all the ICD10 codes, with labels in 'Ann' format as above.

Well this is much longer than Peter's code, but this covers the problem in generality.

FACTOR [labels=!t(C00,C01,D22,C02,C03,C04,D01,D07); values=1...8] ICD10

"Add appropriate numeric levels to this factor"
TXCONST [tIDC10] ICD10 "Convert to text"
TEXT Letters; !T(A,B,C,D,E,F,G,H,I,J,K,L,M,N,O,P,Q,R,S,T,U,V,W,X,Y,Z)
FACTOR [LABELS=Letters] FirstLetter
"Read letters and numbers from text"
READ [C=tIDC10;PRINT=*;LAYOUT=fixed] FirstLetter,Number; FIELD=1,2; SKIP=0; FREP=labels
GROUPS 100*FirstLetter + Number; ICD10 "Combine back again with numeric levels"

PRINT 2(ICD10); FREP=Labels,Levels; DEC=0

ICD10 ICD10
C00 300
C01 301
D07 422
D22 302
C02 303
C03 304
C04 401
D01 407

"Restrict to codes between C02...D10 A=100's, B = 200's, C=200's etc."
RESTRICT ICD10; ICD10 >= 302 .and. ICD10 <= 410
PRINT 2(ICD10); FREP=Labels,Levels; DEC=0

ICD10 ICD10
D22 302
C02 303
C03 304
C04 401
D01 407

Regards, David.
______________________________________________
Dr David Baird Statistical Consultant and GenStat Developer
VSN (NZ) Limited (David@VSN.CO.NZ)
8 Mariposa Crescent, Aidanfield, Christchurch 8025, New Zealand
Ph +64 3 3350588 Cell +64 21 1160803

Post generated using Mail2Forum (http://www.mail2forum.com)

GenStat - RE: Death codes

AndrewMc — Wed, 22 May 2013 00:36:05 GMT

I have a variant on Peter's solution.
Instead of first using EDIT to replace the letters with numbers, I use READ directly to read the first letter of the code into a text (Letter) and the following two digits into a variate (Number). Then, you can restrict separately on the Letter and Number values.

Cheers, Andrew Mc

delete [r=y] ICD10,Letter,Number,tICD10
factor [labels=!t(C00,C01,D22,C02,C03,C04,D01); values=1...5] ICD10

ftext ICD10; tICD10
text Letter
read [CHAN=tICD10; LAYOUT=fixed] Letter,Number; FIELD=1,2

restrict ICD10
restrict ICD10; cond=(Letter.in.'C').and.(Number.in.!(1...50))
print ICD10

From: GENSTAT-Request [mailto:GENSTAT@jiscmail.ac.uk] On Behalf Of Peter Lane
Sent: Wednesday, 22 May 2013 11:40 a.m.
To: GENSTAT@JISCMAIL.AC.UK
Subject: Re: Death codes

The best way I can find requires three steps: factor to text, edit text to replace letters with numbers, then text to variate. But it seems to work.

factor [labels=!t(C00,C01,D22,C02,C03,C04,D01); values=1...5] ICD10

ftext ICD10; tICD10
edit [chan=!t('G/C/1/ G-/D/2/ :')] tICD10
read [chan=tICD10] vICD10

restrict ICD10; vICD10.in.!(101...103)
print ICD10
Peter Lane
Statistical Consultant, Fenner School of the Environment & Society, Australian National University

On 22 May 2013 00:39, Brian Miller <brian.miller@iom-world.org> wrote:
Setting up analyses of cause-specific mortality in GenStat used to be simple when ICD codes for causes were numeric, even for wide ranges of codes, because GenStat has good ways of specifying ranges of values, which can be used with 'in' to restrict to selected groups of causes.

However, now that we're in the era of ICD10 codes, which start with a letter, things are not so straightforward. �Suppose I want to analyse for the group of all neoplasms, i.e. all causes in the range C00-D49. �Can anyone suggest a way of restricting deaths explicitly to this range easily, please? Assume I've created a factor that contains all the ICD10 codes, with labels in 'Ann' format as above.

Thanks
Brian

The contents of this e-mail are confidential and may be subject to legal privilege.
If you are not the intended recipient you must not use, disseminate, distribute or
reproduce all or any part of this e-mail or attachments. If you have received this
e-mail in error, please notify the sender and delete all material pertaining to this
e-mail. Any opinion or views expressed in this e-mail are those of the individual
sender and may not represent those of The New Zealand Institute for Plant and
Food Research Limited.

Post generated using Mail2Forum (http://www.mail2forum.com)

GenStat - RE: Death codes

Anonymous — Tue, 21 May 2013 23:39:54 GMT

The best way I can find requires three steps: factor to text, edit text to replace letters with numbers, then text to variate. But it seems to work.

factor [labels=!t(C00,C01,D22,C02,C03,C04,D01); values=1...5] ICD10

ftext ICD10; tICD10
edit [chan=!t('G/C/1/ G-/D/2/ :')] tICD10
read [chan=tICD10] vICD10

restrict ICD10; vICD10.in.!(101...103)
print ICD10

Peter Lane

Statistical Consultant, Fenner School of the Environment & Society, Australian National University

On 22 May 2013 00:39, Brian Miller <brian.miller@iom-world.org (brian.miller@iom-world.org)> wrote:

Quote:

Setting up analyses of cause-specific mortality in GenStat used to be simple when ICD codes for causes were numeric, even for wide ranges of codes, because GenStat has good ways of specifying ranges of values, which can be used with 'in' to restrict to selected groups of causes.

However, now that we're in the era of ICD10 codes, which start with a letter, things are not so straightforward. Suppose I want to analyse for the group of all neoplasms, i.e. all causes in the range C00-D49. Can anyone suggest a way of restricting deaths explicitly to this range easily, please? Assume I've created a factor that contains all the ICD10 codes, with labels in 'Ann' format as above.

Thanks
Brian

Post generated using Mail2Forum (http://www.mail2forum.com)

GenStat - Death codes

Anonymous — Tue, 21 May 2013 14:39:52 GMT

Setting up analyses of cause-specific mortality in GenStat used to be simple when ICD codes for causes were numeric, even for wide ranges of codes, because GenStat has good ways of specifying ranges of values, which can be used with 'in' to restrict to selected groups of causes.

However, now that we're in the era of ICD10 codes, which start with a letter, things are not so straightforward. Suppose I want to analyse for the group of all neoplasms, i.e. all causes in the range C00-D49. Can anyone suggest a way of restricting deaths explicitly to this range easily, please? Assume I've created a factor that contains all the ICD10 codes, with labels in 'Ann' format as above.

Thanks
Brian

Post generated using Mail2Forum (http://www.mail2forum.com)

GenStat - RE: Plotting histograms with common x-axis settings

AndrewMc — Tue, 14 May 2013 02:40:17 GMT

Hi All

Following David's suggestion to use TRELLIS, I realised that if the TRELLIS procedure does much of what I want, then perhaps I could use and modify code from the procedure to get the graphs that I need. In the example below, I have adapted about a dozen or so lines of code from TRELLIs for calculating histogram limits. The values LHLIM and UHLIM are then passed as options to the DHISTOGRAM command which enables the same limits to be used for more than one graph.

Cheers, Andrew Mc

"Plot 3 histograms on a page all with common x-axis scales"
VARIATE Y1; VALUE=!(13,31,0,54,70,14,19,82,12,55,71,99,9,17,\
95,50,58,27,33,47,52,60,11,30,19,46,44,77,31,10)
VARIATE Y2; VALUE=!(22,9,32,12,61,30)
VARIATE Y3; VALUE=!(42,35,58,71,90,37)
APPEND [NEW=Y; GROUPS=Variate] Y1,Y2,Y3 "stack the data"

"Setup pens, axes and graph frames"
PEN [RESET=y] 2...5; TAREA=100
FRAME [RESET=y] 1
FRAME [RESET=y] 2,3; YUP=1; XLO=0.70,0.85; XUP=0.85,1; \
YLO=0.25; XMLO=0.05; XMUP=0
XAXIS [RESET=y] 2,3; TI=' '"; LPOS=*"
YAXIS [RESET=y] 1,2,3; LPOS=*

"Borrowed code from TRELLIS for calculating histogram limits"
DELETE [REDEF=y] CW,HLIMITS,LHLIM,UHLIM,NG,NGM1,NGROUPS
scalar NGROUPS; 5 "<- user-specified number of groups"
CALC CW=(MAX(Y)-MIN(Y))/NGROUPS
CALC HLIMITS = MIN(Y) + !(1...#NGROUPS)*CW
SCALAR LHLIM; VALUE=MIN(Y)
SCALAR UHLIM; VALUE=MAX(Y)
CALC NG,NGM1=NVAL(HLIMITS)-0,1
IF (HLIMITS$[1]-LHLIM).LE.(HLIMITS$[2]-HLIMITS$[1])/2
CALC LHLIM = 1.5*HLIMITS$[1] - 0.5*HLIMITS$[2]
ENDIF
IF (UHLIM-HLIMITS$[NG]).LE.(HLIMITS$[NG]-HLIMITS$[NGM1])/2
CALC UHLIM = 1.5*HLIMITS$[NG] - 0.5*HLIMITS$[NGM1]
ENDIF

"Plot histograms"
DHIST [WIN=1; KEYW=0; TI='Y1'; SCREEN=clear; ORIENT=h; \
NGROUP=NGROUPS; low=LHLIM; upp=UHLIM] Y1; PEN=2
DHIST [WIN=2; KEYW=0; TI='Y2'; SCREEN=keep; ORIENT=h; \
NGROUP=NGROUPS; low=LHLIM; upp=UHLIM] Y2; PEN=3
DHIST [WIN=3; KEYW=0; TI='Y3'; SCREEN=keep; ORIENT=h; \
NGROUP=NGROUPS; low=LHLIM; upp=UHLIM] Y3; PEN=4

-----Original Message-----
From: GENSTAT-Request [mailto:GENSTAT@jiscmail.ac.uk] On Behalf Of Andrew McLachlan
Sent: Tuesday, 14 May 2013 9:13 a.m.
To: GENSTAT@JISCMAIL.AC.UK
Subject: Re: Plotting histograms with common x-axis settings

Thanks David. I had not thought of this approach. Using TRELLIS certainly does provide histograms with a common x-axis. A limitation is that it is more work to (manually) customise the individual plots.

Cheers, Andrew Mc

-----Original Message-----
From: GENSTAT-Request [mailto:GENSTAT@jiscmail.ac.uk] On Behalf Of David Baird
Sent: Monday, 13 May 2013 4:35 p.m.
To: GENSTAT@JISCMAIL.AC.UK
Subject: Re: Plotting histograms with common x-axis settings

Hi Andrew;

If you stack the data, you could use TRELLIS that does this for you:

VARIATE Y1; VALUE=!(13,31,0,54,70,14,19,82,12,55,71,99,9,17,\
95,50,58,27,33,47,52,60,11,30,19,46,44,77,31,10)
VARIATE Y2; VALUE=!(22,9,32,12,61,30)

APPEND [NEW=Y; GROUPS=Variate] Y1,Y2
TRELLIS [GROUPS=Variate] Y; METHOD=histogram

Regards, David.
______________________________________________
Dr David Baird�� Statistical Consultant and GenStat Developer VSN (NZ) Limited (David@VSN.CO.NZ)
8 Mariposa Crescent, Aidanfield, Christchurch 8025, New Zealand Ph +64 3 3350588 Cell +64 21 1160803

Quote:

-----Original Message-----
From: GENSTAT-Request [mailto:GENSTAT@jiscmail.ac.uk] On Behalf Of
Andrew McLachlan
Sent: 13 May 2013 4:14 p.m.
To: GENSTAT@JISCMAIL.AC.UK
Subject: Plotting histograms with common x-axis settings

Hi All

I am trying to plot several histograms with common x-axis settings
(axis

limits

Quote:

and tickmarks). This is relatively easy to do when plotting one-way
tables

Quote:

values - simply use the same factor to define each of the one-way tables.
However, when plotting a histogram for a variate, DHISTOGRAM creates
the data bins, sets the axis limits and the placement of the
tickmarks. Since

I am

Quote:

plotting histograms separately for several variates, I am having

difficulty in

Quote:

matching the x-axis settings from one graph to the other. While I can

easily set

Quote:

the same number of bins for each histogram, does anyone have any
suggestions on how I can get x-axis tickmarks to match from one
histogram

Quote:

the next. I have included some example code to plot two histograms.

VARIATE Y1; VALUE=!(13,31,0,54,70,14,19,82,12,55,71,99,9,17,\
95,50,58,27,33,47,52,60,11,30,19,46,44,77,31,10)
VARIATE Y2; VALUE=!(22,9,32,12,61,30)
PEN [RESET=y] 2...5; TAREA=100
FFRAME [ROW=2; COL=1; RSKIP=0.02; YMLO=0.1] YAXIS [RESET=y] 1,2;
LPOS=* XAXIS [RESET=y] 2,1; TITLE='Base x-axis','Need to match scale
below'
DHISTOGRAM [WIN=2; KEYW=0; SCREEN=clear; NGROUP=5] Y1; PEN=2
DHISTOGRAM [WIN=1; KEYW=0; SCREEN=keep; NGROUP=5] Y2; PEN=3

Cheers, Andrew Mc

Andrew McLachlan
Biometrician
Plant & Food Research

T: +64 6 355 6110
F: +64 6 351 7050
andrew.mclachlan@plantandfood.co.nz
www.plantandfood.co.nz
The New Zealand Institute for Plant & Food Research Limited

Postal Address: Plant & Food Research Palmerston North Private Bag
11600, Manawatu Mail Centre, Palmerston North, 4442, New Zealand
Physical Address: Plant & Food Research Palmerston North Batchelar
Road, Palmerston North, 4474, New Zealand

The contents of this e-mail are confidential and may be subject to
legal privilege.
If you are not the intended recipient you must not use, disseminate,

distribute

Quote:

or
reproduce all or any part of this e-mail or attachments. If you have

received

Quote:

this
e-mail in error, please notify the sender and delete all material

pertaining to

Quote:

this
e-mail. Any opinion or views expressed in this e-mail are those of
the individual sender and may not represent those of The New Zealand
Institute for Plant and Food Research Limited.

The contents of this e-mail are confidential and may be subject to legal privilege.
If you are not the intended recipient you must not use, disseminate, distribute or reproduce all or any part of this e-mail or attachments. If you have received this e-mail in error, please notify the sender and delete all material pertaining to this e-mail. Any opinion or views expressed in this e-mail are those of the individual sender and may not represent those of The New Zealand Institute for Plant and Food Research Limited.

The contents of this e-mail are confidential and may be subject to legal privilege.
If you are not the intended recipient you must not use, disseminate, distribute or
reproduce all or any part of this e-mail or attachments. If you have received this
e-mail in error, please notify the sender and delete all material pertaining to this
e-mail. Any opinion or views expressed in this e-mail are those of the individual
sender and may not represent those of The New Zealand Institute for Plant and
Food Research Limited.

Post generated using Mail2Forum (http://www.mail2forum.com)

GenStat - RE: Plotting histograms with common x-axis settings

AndrewMc — Mon, 13 May 2013 21:12:57 GMT

Thanks David. I had not thought of this approach. Using TRELLIS certainly does provide histograms with a common x-axis. A limitation is that it is more work to (manually) customise the individual plots.

Cheers, Andrew Mc

-----Original Message-----
From: GENSTAT-Request [mailto:GENSTAT@jiscmail.ac.uk] On Behalf Of David Baird
Sent: Monday, 13 May 2013 4:35 p.m.
To: GENSTAT@JISCMAIL.AC.UK
Subject: Re: Plotting histograms with common x-axis settings

Hi Andrew;

If you stack the data, you could use TRELLIS that does this for you:

VARIATE Y1; VALUE=!(13,31,0,54,70,14,19,82,12,55,71,99,9,17,\
95,50,58,27,33,47,52,60,11,30,19,46,44,77,31,10)
VARIATE Y2; VALUE=!(22,9,32,12,61,30)

APPEND [NEW=Y; GROUPS=Variate] Y1,Y2
TRELLIS [GROUPS=Variate] Y; METHOD=histogram

Regards, David.
______________________________________________
Dr David Baird�� Statistical Consultant and GenStat Developer VSN (NZ) Limited (David@VSN.CO.NZ)
8 Mariposa Crescent, Aidanfield, Christchurch 8025, New Zealand Ph +64 3 3350588 Cell +64 21 1160803

Quote:

limits

Quote:

and tickmarks). This is relatively easy to do when plotting one-way
tables

Quote:

I am

Quote:

plotting histograms separately for several variates, I am having

difficulty in

Quote:

matching the x-axis settings from one graph to the other. While I can

easily set

Quote:

the same number of bins for each histogram, does anyone have any
suggestions on how I can get x-axis tickmarks to match from one
histogram

Quote:

distribute

Quote:

or
reproduce all or any part of this e-mail or attachments. If you have

received

Quote:

this
e-mail in error, please notify the sender and delete all material

pertaining to

Quote:

this
e-mail. Any opinion or views expressed in this e-mail are those of
the individual sender and may not represent those of The New Zealand
Institute for Plant and Food Research Limited.

GenStat - RE: Plotting histograms with common x-axis settings

david@vsn.co.nz — Mon, 13 May 2013 04:35:38 GMT

Hi Andrew;

If you stack the data, you could use TRELLIS that does this for you:

VARIATE Y1; VALUE=!(13,31,0,54,70,14,19,82,12,55,71,99,9,17,\
95,50,58,27,33,47,52,60,11,30,19,46,44,77,31,10)
VARIATE Y2; VALUE=!(22,9,32,12,61,30)

APPEND [NEW=Y; GROUPS=Variate] Y1,Y2
TRELLIS [GROUPS=Variate] Y; METHOD=histogram

Regards, David.
______________________________________________
Dr David Baird�� Statistical Consultant and GenStat Developer
VSN (NZ) Limited (David@VSN.CO.NZ)
8 Mariposa Crescent, Aidanfield, Christchurch 8025, New Zealand
Ph +64 3 3350588 Cell +64 21 1160803

Quote:

limits

Quote:

and tickmarks). This is relatively easy to do when plotting one-way tables

Quote:

values - simply use the same factor to define each of the one-way tables.
However, when plotting a histogram for a variate, DHISTOGRAM creates the
data bins, sets the axis limits and the placement of the tickmarks. Since

I am

Quote:

plotting histograms separately for several variates, I am having

difficulty in

Quote:

matching the x-axis settings from one graph to the other. While I can

easily set

Quote:

the same number of bins for each histogram, does anyone have any
suggestions on how I can get x-axis tickmarks to match from one histogram

Quote:

the next. I have included some example code to plot two histograms.

VARIATE Y1; VALUE=!(13,31,0,54,70,14,19,82,12,55,71,99,9,17,\
95,50,58,27,33,47,52,60,11,30,19,46,44,77,31,10)
VARIATE Y2; VALUE=!(22,9,32,12,61,30)
PEN [RESET=y] 2...5; TAREA=100
FFRAME [ROW=2; COL=1; RSKIP=0.02; YMLO=0.1]
YAXIS [RESET=y] 1,2; LPOS=*
XAXIS [RESET=y] 2,1; TITLE='Base x-axis','Need to match scale below'
DHISTOGRAM [WIN=2; KEYW=0; SCREEN=clear; NGROUP=5] Y1; PEN=2
DHISTOGRAM [WIN=1; KEYW=0; SCREEN=keep; NGROUP=5] Y2; PEN=3

Cheers, Andrew Mc

Andrew McLachlan
Biometrician
Plant & Food Research

T: +64 6 355 6110
F: +64 6 351 7050
andrew.mclachlan@plantandfood.co.nz
www.plantandfood.co.nz
The New Zealand Institute for Plant & Food Research Limited

Postal Address: Plant & Food Research Palmerston North
Private Bag 11600, Manawatu Mail Centre, Palmerston North, 4442, New
Zealand
Physical Address: Plant & Food Research Palmerston North
Batchelar Road, Palmerston North, 4474, New Zealand

The contents of this e-mail are confidential and may be subject to legal
privilege.
If you are not the intended recipient you must not use, disseminate,

distribute

Quote:

or
reproduce all or any part of this e-mail or attachments. If you have

received

Quote:

this
e-mail in error, please notify the sender and delete all material

pertaining to

Quote:

this
e-mail. Any opinion or views expressed in this e-mail are those of the
individual
sender and may not represent those of The New Zealand Institute for Plant
and
Food Research Limited.

Post generated using Mail2Forum (http://www.mail2forum.com)