GAMs V1.0 Generalized Additive Models (through the mgcv R package)

[#1]

HI everyone,

I have been a little frustrated with the fact that JMP still does not have a GAM modeling platform (even if it exists in SAS for many years) so I decided to create an add-in that is based on the mgcv package in R. So to run this add-in you will need to have R v. 3.3.3 running on your computer. It basically creates a JMP dialog window with many of the options available in the R package and retrieves the essential results from R and adds some graphs and options. If you need additional features, please let me know. Here are some screen shots of the interface and results window.

Best regards, Yves

[Raaed]

Dear Sir

Please, update this add-ins to run in R 3.6

[yvesprairie]

Dear Raewdy,

Unfortunately, this has nothing to do with the script, it is a limitation of the JMP R integration functions. Only JMP can fix that. Presumably, they will fix it in the next version of JMP... That being said, you can have multiple versions of R on your computer, including R 3.3.3.

Regards, Yves

[yvesprairie]

Hi Raewdy,

Here is the script for it. Regards, Yves

Names Default To Here( 1 );
Clear Globals();
//
//	Check host and R connection
//
macsys = Host is( "Mac" );
winsys = Host is( "Windows" );

connected = R Is Connected();

R Term();
If( connected == 0, 

	If( macsys == 1,
		Set Environment Variable( "R_HOME", "/Library/Frameworks/R.framework/Versions/Current/Resources" )
	);
	If( winsys == 1,
	
		ppath = Pick Directory( "Locate the directory of your R program version 3.3.3" );
		ppath=substr(ppath,2,length(ppath)-1);
	 
		part = Parse( Eval Insert( "
	Set Environment Variable( \!"R_HOME\!",\!"^ppath^\!" );" ) );
		Eval( part );
	);
	
);
R Init();
connected = R Is Connected();
R Term();
If( connected == 0,
	ww2 = New Window( "Warning", 
	
		Spacer Box(),
		Text Box( "    CANNOT ESTABLISH CONNECTION WITH R    " ),
		Spacer Box(),
		H Center Box(
			Button Box( "OK",
				ww2 << close window;
				
			) 
		
		)
	);
	Stop();
);
	
//
//   CHECK R VERSION
//				


R Init();
connected = R Is Connected();

If( connected == 1,
	Rver = Char( R Get Version() );
	If( Rver == "[3 3 3]",
		,
		ww1 = New Window( "Warning",
			<<Modal, 
		
			Spacer Box(),
			Text Box( "    THIS ADD-IN REQUIRES R version 3.3.3    " ),
			Spacer Box(),
			H Center Box(
				Button Box( "OK", ww1 << close window )

	
		
			)
		);
		Set Environment Variable( "R_HOME", "C:/Program Files/R/R-3.3.3/" );
		Stop();
	);
);

dt = Current Data Table();
If( Is Empty( dt ),
	Try( dt = Open(), Throw( "No data table found" ) )
);
allvars = dt << get column names( string );
families = ({"gaussian", "binomial", "Gamma", "poisson"});
meths = ({"REML", "GCV Cp", "GACV Cp", "ML"});
aa = Associative Array( {1, 2, 3, 4}, {"REML", "GCV.Cp", "GACV.Cp", "ML"} );
aa2 = Associative Array( {1, 2, 3, 4}, {"gaussian", "binomial", "Gamma", "poisson"} );
aa3 = Associative Array(
	{"gaussian", "binomial", "Gamma", "poisson"},
	{lnkgauss, lnkbinomial, lnkgamma, lnkpoisson}
);
gamma = 1;
lnkgauss = ({"identity", "log", "inverse"});
lnkbinomial = ({"logit", "probit", "cauchit", "log", "cloglog"});
lnkgamma = ({"identity", "log", "inverse"});
lnkpoisson = ({"identity", "log", "sqrt"});
tt = lnkgauss;
sellnk = "identity";
selfam = "gaussian";
pog = 0;
cw = 0;

nw = New Window( "GAM Modelling", 
	
	H List Box(
		V List Box(
			Panel Box( "Select Columns", MDSColList = Col List Box( All, width( 140 ), nlines( 10 ) ) ),
			H List Box(
				Outline Box( "Family",
					rb1 = Combo Box(
						families,
						selfam = rb1 << getselected();
						tt = aa3 << getvalue( selfam );
						cb << set items( tt );
			
					)
				), 
				
				ob = Outline Box( "Link",
					cb = Combo Box(
						tt,
						sellnk = cb << Getselected();
			
					)
				)
			), 
						
			H List Box(
				Outline Box( "Method", rb2 = Radio Box( meths ) ),
				Outline Box( "Smoothness",
					H List Box( tb = Text Box( "Gamma: " ), sb = Number Edit Box( gamma, 3 ) ), 

				)
			), 

		)
						
	,
		V List Box(
			Panel Box( "Cast Selected Columns into Roles",
				Lineup Box( N Col( 2 ), Spacing( 3 ),
					Button Box( "Y", GAM_YVar << Append( MDSColList << GetSelected ) ),
					GAM_YVar = Col List Box( width( 140 ), nLines( 1 ), MinItems( 1 ) ),
					Mbb = Button Box( "Linear X variables", GAM_L_XVar << Append( MDSColList << GetSelected ) ),   //for SAS - matrix var to specify diff subjects;
					GAM_L_XVar = Col List Box( width( 140 ), nLines( 5 ) ),
					Button Box( "Smoothed X variables", GAM_S_XVar << Append( MDSColList << GetSelected ) ),
					GAM_S_XVar = Col List Box( width( 140 ), nLines( 5 ), MinItems( 1 ) ), 
								
								
					Button Box( "Remove",
						GAM_YVar << RemoveSelected;
						GAM_L_XVar << RemoveSelected;
						GAM_S_XVar << RemoveSelected;
					
					), 

				), 
			
			),
			tb5 = Text Box( "Options:" ),
			cbp = Check Box( {"Show data points in smoother plots"}, pog = cbp << Get() ),
			kdo = Check Box( {"Keep dialog open"} ),
			kdo << set( 1, 1 );
			cw = kdo << Get();
		
		
		, 					
//	), //end HListBox
					
			Text Box( " " ),
			H List Box(
				Button Box( "Run",
					Yvar = GAM_YVar << Get Items;
//			
					S_X_vars = GAM_S_XVar << Get Items;
					L_X_vars = GAM_L_XVar << Get Items;
					fam = rb1 << Get;
					method = rb2 << Get;
					npterms = N Items( L_X_vars ) + 1;
					Nxvar = N Items( S_X_vars );
					gamma = sb << get;
					If( pog == 1,
						Rpog = "TRUE",
						Rpog = "FALSE"
					);
			//
					//	Main portion
					//
			
			
					xvars = S_X_vars;
					xcode = "";
					R Init();
					pterms = List( "intercept" );
					TYV = Dt << Get as matrix( Column( Yvar ) );
		//			LY = Min( YV ) * .95;
		//			HY = Max( YV ) * 1.05;
					ndt = N Row( dt );
					SV = "";
					LV= "";
					

					R Send( dt );
					R Submit( "Rnames<-names(dt)" );
					RnamesinJMP = R Get( Rnames );
					aan = Associative Array( Allvars, RnamesinJMP );
					sterms = List( "" );
			
					RYvar = aan << Get value( Yvar[1] );

			
					For( i = 1, i < npterms, i++, 

						LXvar = aan << Get value( L_X_vars[i] );
						Insert Into( pterms, Char( L_X_vars[i] ) );
						xcode ||= Eval Insert( "^LXvar^+" );
				
					);
					For( i = 1, i <= Nxvar, i++, 

						RXvar = aan << Get value( xvars[i] );
						Insert Into( sterms, Char( xvars[i] ) );
						xcode ||= Eval Insert( "s(^RXvar^)+" );
		
					);
					TSV = dt << Get as matrix( S_X_vars );
					TLV = dt << Get as matrix( L_X_vars );
					TAV = TYV;
					If( npterms > 1, TAV ||= TLV );
					If( Nxvar >= 1,
						TAV ||= TSV,
						Print( "Please select variables" )
					);
					nm = Loc Nonmissing( TAV );
					
					nnm=n row(nm);
					SV=TSV[nm];
					LV=TLV[nm];
					YV=TYV[nm];
					LY = Min( YV ) * .95;
					HY = Max( YV ) * 1.05;
					
					if(nnm<npterms+Nxvar, 
					ww3 = New Window( "Warning", 
	
		Spacer Box(),
		Text Box( "   FEWER OBSERVATIONS THAN PARAMETERS    " ),
		Spacer Box(),
		H Center Box(
			Button Box( "OK",
				ww3 << close window;
				
			) 
		
		)
	);
					stop();
						
					);
			
					sterms = sterms[2 :: nxvar + 1];
	
					xcode = Substr( xcode, 1, Length( xcode ) - 1 );
		
					Rcode = "library(mgcv)";
					Rcode ||= Eval Insert(
						"\[
b <- gam(^RYvar^~^xcode^,data=dt,method="^aa<<getvalue(method)^", family="^aa2<<getvalue(fam)^(link=^sellnk^)",gamma=^gamma^)

sumb<-summary(b)
fv<-b$fitted.values
rsq<-as.list(sumb$r.sq)
dev<-as.list(sumb$dev.expl)	
ptable<-as.list(sumb$p.table)
spcrit<-as.list(sumb$sp.criterion)
stable<-as.list(sumb$s.table)
nobs<-as.list(sumb$n)
parest<-as.list(sumb$p.coeff)
pt<-(sumb$p.t)
plot.gam(b,pages=1,residuals=^Rpog^,cex=2,all.terms=TRUE,pers=TRUE, shade=TRUE, pch=20)
]\"
					);

		
			R Init();
					R Submit( Rcode );
			
	
					SmPl = R Get Graphics( "png" );
	
					ptable = R Get( ptable );
					stable = R Get( stable );
					spcrit = R Get( spcrit );
					rsq = R Get( rsq );
					dev = R Get( dev );
					nobs = R Get( nobs );
					fitvalues = J( N Row( TYV ), 1, . );
					tfv = R Get( fv );
					fitvalues = tfv;
					LYH = Min( fitvalues ) * .95;
					HYH = Max( fitvalues ) * 1.05;
					Resids = TYV[nm] - transpose(fitvalues);
					LR = Min( Resids ) * 1.05;
					HR = Max( Resids ) * 1.05;
					If( Abs( LR ) >= HR,
						HR = Abs( LR ),
						LR = -HR
					);
	
					R Term();
	
	
					genfit = rsq || dev || nobs;
					pest = J( npterms, 1, . );
					pstderr = J( npterms, 1, . );
					pt = J( npterms, 1, . );
					pprob = J( npterms, 1, . );

					For( i = 1, i <= npterms, i++,
						pest[i] = ptable[i];
						pstderr[i] = ptable[i + 1 * npterms];
						pt[i] = ptable[i + 2 * npterms];
						pprob[i] = If( ptable[i + 3 * npterms] <= 0.0001,
							0.0001,
							ptable[i + 3 * npterms]
						);
	
					);
	

	
					sedf = J( nxvar, 1, . );
					srdf = J( nxvar, 1, . );
					sf = J( nxvar, 1, . );
					sprob = J( nxvar, 1, . );
	
	
					For( i = 1, i <= nxvar, i++,
						sedf[i] = stable[i];
						srdf[i] = stable[i + nxvar * 1];
						sf[i] = stable[i + 2 * nxvar];
						sprob[i] = If( stable[i + 3 * nxvar] <= 0.0001,
							0.0001,
							stable[i + 3 * nxvar]
						);
	
					);
	
			
					New Window( "GAM Model Fit",
						Outline Box( Eval Insert( "Family: ^selfam^" ) ),
						Outline Box( Eval Insert( "Link function: ^sellnk^" ) ),
						Outline Box( "Whole Model",
							{"Save predicted", New Column( "Predicted values", values( fitvalues ) ),
							"Save residuals", New Column( "Residuals", values( YV - fitvalues ) )},
							H List Box(
								V List Box(
									Outline Box( "Obs vs Pred.",
										Graph Box(
											X Scale( LYH, HYH ),
											Y Scale( LY, HY ),
											Yname( "Obs Y" ),
											Xname( "Pred Y" ),
											Frame Size( 250, 250 ),
											Line( {LY, LY}, {HY, HY} ),
											Marker( Marker State( 0 ), fitvalues, YV )
										)
									),
									Outline Box( "Residuals vs Pred.",
										Graph Box(
											Y Scale( LR, HR ),
											X Scale( LYH, HYH ),
											Yname( "Residuals" ),
											Xname( "Pred Y" ),
											Frame Size( 250, 120 ),
											Line( {LYH, 0}, {HYH, 0} ),
											Marker( Marker State( 0 ), fitvalues, Resids )
										)
									)
								),
								Outline Box( "Smoother and Linear leverage plots", Picture Box( SmPl ) )
							)
						), 
	
						Outline Box( "Summary of Fit", 
		
							Table Box(
								String Col Box( , {"adj. R-square", "Deviance", "Number of Obs."} ),
								Number Col Box( , genfit )
							)
						), 
		
						Outline Box( "Parametric terms estimates", 
		
							Table Box(
								String Col Box( "Term", pterms ),
								Number Col Box( "Estimate", pest ),
								Number Col Box( "Std Error", pstderr ),
								Number Col Box( "Prob.", pprob ), 

							)
						),
						Outline Box( "Smoothed terms estimates", 
		
							Table Box(
								String Col Box( "Term", sterms ),
								Number Col Box( "Est. df", sedf ),
								Number Col Box( "Ref df", srdf ),
								Number Col Box( "F", sf ),
								Number Col Box( "Prob.", sprob ), 

							)
						)	
			
		
					);
					If( cw == 0,
						nw << close window
					);
				)
			)
		)
	)
	
);

[Raaed]

Dear sir

Thank you very much, I will make some edit to run in R-3.6.1, and sending here

Regards, Raewdy

[inn]

Hi Yves @yvesprairie , 

Thank you for posting this script, I run my GAMs through this. 

Is it possible to add to the script a method of computing AIC?

Kind regards, 

Ingrid Rognes

[yvesprairie]

Hi Ingrid,

Sorry I hadn't seen your post. If mccv outputs it then yes it would be easy to add it. I'll look into it.

Thanks for your interest, Yves

[ptadger]

Hello Yves Praire!

Really interesting Addin, I wonder if is compatible with JMP Pro 16.1?

I give me an error when I tried to install it.

[yvesprairie]

Hi Ptadger,

I don't use it very often but I just tried it and it works well with 16. Not sure what could be the problem. Have you tried just the script as opposed to the add-in?

Yves

[ptadger]

Hi  Yves ! 

Thank you for such a prompt response!!

Indeed, it works with the script (maybe the reason is that I didn't have R 3.3 install before, now I have it). I wonder if you can recommend me something for my ordinary predictors that have a good fit as a SASH or Zero-inflated Poisson). I have been reading that using a penalized spline can be used in ordinary using a specific amount of knots (related with the number of ordinary levels).

[Nickki_T]

@yvesprairie This add-in doesn't work for JMP17, correct? 

[yvesprairie]

@Nickki_T The script will work with JMP17 , it is the R version that might be problematic. Also try the script instead of the add-in, might be easier.

[Nickki_T]

Thanks, @yvesprairie. I got the GAM modelling box open with the script, Thanks. However, when click run, I got "Send data for "dt" failed. Why is that? 

[diriba]

I'm trying to use Generalized Additive model in JMP18 student edition. I received the following error:
Add_in "ca.uqam prairie.gam - GAMs V1.0" is not valid for this version of JMP while trying to install the add-in. For which version of JMP is it valid? or How can I get it for my JMP version?
Thanks

This issue is resolved using the code. But I wonder how can I compare between treatments using the GAM platform like in graph builder similar to this: