title: Piecewise nonlinear regression by MLE in R software

text: I have been working with nonlinear models whose parameters have been estimated by Maximum Likelihood. Also, I use the "Vuong" test for model selection in non-nested models, this test requires to have every single value of log- likelihood of the sample. Lately, I found the GNM function really useful for estimating models by MLE model. The problem is that in this function is difficult (at least for me) to perform piecewise regression, it means, only one variable or term of the model changing its performance along the time. My data:                 1.29023348  0.55990170 -1.81628062  0.04587681  0.06915447 -0.64716236      -0.39194582 -0.17585932 -0.11706221 -0.05392234  0.94070437  0.84266829       2.12744151  1.72599850 -0.77428606  0.31521075  1.79509484  0.85176806      -0.31293698  0.39605241  1.63390712  0.82201408  2.32559314  1.62486421       2.45086803  1.71495391  1.52439007  4.18112171  1.01513944  4.10700194       4.81114802  2.77305277  4.38637816  3.37523675  4.41486264  4.76289968       4.22082055  4.93576180  3.39264254  4.17968804  3.40033996  4.36463976       4.34403336  4.24994731  5.28361166  6.02246253  5.88749010  5.60443358       6.52085529  5.02354714  4.52130778  5.66370450  7.28031355  5.28144249       4.77022080  4.70795500  4.52891649  7.88707463  7.38768136  5.18242267       My models are:               1) Y = Alpha1*time + Alpha2*sin(2*pi*time/Alpha3) + Alpha4*(-1)^time          2) Y = Alpha0(time>=30) + Alpha1*time + Alpha2*sin(2*pi*time/Alpha3) + Alpha4*(-1)^time          3) Y = Alpha0(time>=30) + Alpha1*time*(time<30) + Alpha2*time*(time>=30) +             Alpha3*sin(2*pi*time/Alpha4) + Alpha5*(-1)^time      `Alpha0`,`Alpha1`,`Alpha2`,`Alpha3`,`Alpha4`,`Alpha5` are the parameters to estimate by MLE. So far, my code is:               time <- seq(1,60)     CCPR <- function(x) {       list(predictors <- list(alpha1=1,alpha2=1,alpha3=1,alpha4=1),            variables <- list(substitute(x)) ,            term <- function(predictors,variables) {              paste(predictors[1], "*", variables[1], "+",                    predictors[2], "*sin(2*pi*", variables[1], "/",                    predictors[3], ")+", predictors[4], "*(-1)^", variables[1]                   )             } )     }     class(CCPR) <- "nonlin"     modelo.general <- gnm(ye~CCPR(time))       This is the code for the first model, and I do not know how to add the breakpoint in the second and third model. I might have to include a dummy variable, if so, how should I do it in the GNM function?


title: Using R-Studio for data mining to compare different data sets against one another, asynchronously

text: For a project, I've collected some basic information on weather, and some basic information on births. Both data sets are focused on a single state, Massachusetts. The weather set contains very basic information: the weather event (hail, winter storm, etc.), and the month it occurred in. Not very granular. Similarly, the birth records are merely total births in a given month. The birth data set contains records from 2007-2011. The winter weather set contains records from 1996-2013. So I'll only be using some of the winter weather records. The goal of this project is somewhat tongue-in-cheek, but serious in nature. I hope to see if there is a correlation between serious winter weather events (in a given month), and a spike in births nine months later; i.e., do bad storms keep people indoors, and does that lead to more births. I would expect to see an increase from whatever the baseline is, when there is a similar spike in serious winter weather events. I'm not looking for anybody to do this project for me, but hints and suggestions are very welcome. I was planning to use associative rules mining (and the "arules" package) to accomplish this goal, but I haven't had any success. I guess my data isn't organized as "transaction data", so I'm not sure what type of data mining technique to use here (I've only had limited experience with arules). **Any suggestions?** Using rbind, I've already combined the two data sets into two distinct files, which contain all the relevant information. **combined.weather1** <\- includes three variables (YEAR, MONTH, EVENT_TYPE), and 2750 observations. It contains info from 1996-2013. The YEAR variable is integer, while the MONTH and EVENT_TYPE variables are factors. **birth.data** <\- includes three variable (YEAR, MONTH, BIRTH_TOTAL), and sixty observations. It contains data from 2007-2011. The YEAR variable is integer, while the MONTH and BIRTH_TOTAL variables are factors. To further complicate things, I'd like to rule out some of the less serious winter weather events, such as "Hail", "Frost/Freeze", and "Winter Weather". Worth noting, there are no missing or NA values in any of these records. The data is clean and complete (and took forever to prepare!). So I need to select a data mining technique, first and foremost. From there, I would very much appreciate if anybody can offer help on how to set it up to compare these two simple data sets, while factoring in a nine month difference between certain winter weather events and (possible) spikes in births. EDIT: adding examples, per request. dput(combined.weather1[1:10],)               structure(list(YEAR = c(1996L, 1996L, 1996L, 1996L, 1996L, 1996L,      1996L, 1996L, 1996L, 1996L), MONTH_NAME = structure(c(5L, 5L,      5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L), .Label = c("April    ", "August   ",      "December ", "February ", "January  ", "July     ", "June     ",      "March    ", "May      ", "November ", "September", "October  ",      "April", "August", "December", "February", "January", "July",      "June", "March", "May", "November", "October"), class = "factor"),          EVENT_TYPE = structure(c(2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,          2L, 2L), .Label = c("Hail", "Heavy Snow", "Winter Storm",          "Winter Weather", "Ice Storm", "Frost/Freeze", "WINTER WEATHER",          "Blizzard"), class = "factor")), .Names = c("YEAR", "MONTH_NAME",      "EVENT_TYPE"), row.names = c(NA, 10L), class = "data.frame")      dput(birth.data[1:10,])               structure(list(YEAR = c(2007L, 2007L, 2007L, 2007L, 2007L, 2007L,      2007L, 2007L, 2007L, 2007L), MONTH = structure(c(5L, 4L, 8L,      1L, 9L, 7L, 6L, 2L, 12L, 11L), .Label = c("April", "August",      "December", "February", "January", "July", "June", "March", "May",      "November", "October", "September"), class = "factor"), BIRTH_TOTAL = structure(c(29L,      9L, 49L, 27L, 57L, 52L, 58L, 59L, 36L, 41L), .Label = c("5,414",      "5,459", "5,515", "5,627", "5,667", "5,730", "5,754", "5,801",      "5,833", "5,853", "5,925", "5,976", "5,979", "5,981", "6,040",      "6,051", "6,070", "6,085", "6,135", "6,158", "6,162", "6,194",      "6,198", "6,208", "6,212", "6,221", "6,227", "6,242", "6,250",      "6,260", "6,261", "6,318", "6,341", "6,342", "6,380", "6,385",      "6,396", "6,438", "6,444", "6,459", "6,466", "6,469", "6,489",      "6,506", "6,509", "6,510", "6,520", "6,531", "6,570", "6,583",      "6,616", "6,735", "6,781", "6,803", "6,820", "6,834", "6,858",      "6,933", "7,291"), class = "factor")), .Names = c("YEAR", "MONTH",      "BIRTH_TOTAL"), row.names = c(NA, 10L), class = "data.frame")      EDIT 2: I have completed the suggestions below. Each data frame has been further refined (eliminated spaces, converting to numeric, adding DATE variables, etc.). I'm now attempting to use the ccf() function to find a correlation between extreme winter weather events, and (possible) spikes in births nine months later. I have included below my stripped down basic code. Right now, I've written two attempts at ccf(): the first one uses the newly created birth.data$DATE, and compares to combined.weather$DATE; the second is the same, except that it adds a second variable, combined.weather$EVENT_TYPE. So I've made some basic progress. I just figured out to create a different weather data frame, using rbind again, with just the years I want:               combined.weather.birth <- rbind(winter2007,winter2008,winter2009,winter2010,winter2011)      My other concern - and this is basic programming but I'm not sure how to do this - is **how to filter out the less extreme winter weather events, so that they don't unduly influence my correlation model** i.e. a simple frost/freeze warning is probably insufficient to keep people indoors in the winter, and thus lead to a spike in births. Right now frost/freeze and hail are weighted the same as winter storm, and I need to get address that. **Any tips**? NOTE: please note, I've gotten rid of combined.weather1, and am now working exclusively with combined.weather and combined.weather.birth, to keep things simple. Here is some of the new code:               birth.data$DATE = ymd(paste0(birth.data$YEAR, '-', birth.data$MONTH, "-", 01)) #adds DATE variable     combined.weather$DATE = ymd(paste0(combined.weather$YEAR, '-', combined.weather$MONTH, "-", 01)) #adds DATE variable     #Create special weather comparison for birth records:     combined.weather.birth <- rbind(winter2007,winter2008,winter2009,winter2010,winter2011)     combined.weather.birth$DATE = ymd(paste0(combined.weather.birth$YEAR, '-', combined.weather.birth$MONTH, "-", 01)) #adds DATE variable         correlation1 <- ccf(birth.data$DATE, combined.weather.birth$DATE, lag.max = NULL, type="correlation", plot=TRUE)         correlation2 <- ccf(birth.data$DATE, combined.weather.birth$DATE+combined.weather$EVENT_TYPE, lag.max = NULL, type="correlation", plot=TRUE)


title: How to split a numeric variable into a binary low-high variable

text: I have measured frequency of a certain behavior on 15 individuals. I would like to create two groups based on the amount of this behaviour that was observed (i.e., a group exhibiting high levels of the behaviour and a group exhibiting low levels of the behaviour). I want to see whether this new binary variable predicts a dependent variable that I have measured.


title: What are the null and alternate hypotheses? Will you use a left, right, or two-tailed test?

text: > Ben has a coin which he claims is weighted in a way so that when he flips > it, heads appears more often than 50% of the time. He tries to prove it to > you by flipping the coin 100 times, which results in 60 heads. > > Test the claim that the population proportion of heads that result from > flipping this coin is actually more than 50%. Use a 1% level of > significance. > > Identify the sampling distribution you will use, what is the value of the > test statistics, then find P-value. I know the level of significance is 0.01; I'm just stuck on how to go about finding the null and alternate hypotheses. $\rm H_0$: would it be μ = 60, since that's the average of how many heads he got, and$\rm H_1: μ > 60$? Thank you.


title: Bias from increased information in FLAME clustering

text: I hope this is an appropriate forum for this question...if not, any pointers on a place to ask would be great. If my questions is not clear, please just let me know and I'll try to add information/explanation where I can. Say I have a set of data points that I run through a FLAME clustering algorithm, and I get my set of cluster supporting objects, and the fuzzy memberships among them. Next, let's say I get an additional set of data, for argument's sake, let's say it doubles the amount of data. If I want to add that to my clustering, without having to reprocess the old data (i.e. just using the cluster supporting objects, and adding these new observations to the cluster), will I be introducing some kind of bias? To get the most "accurate" picture, would it be better to reprocess ALL the data, re-establish CSO's, re-compute distances, etc.? To add to the complexity, let's say I'm going to be periodically adding similar quantities of data points. Does the answer to the above change? And I guess a final question - is there a clustering algorithm (or family of algorithms) that is not subject to extreme bias in the beginning, such that added data points can be added relatively pain-free? I'm sort of thinking that FLAME clustering, with its nice fuzziness properties is a reasonable way of going about this, but any more suggestions would be helpful. Thanks!


title: Clustering quality: normalized mean square error or absolute error?

text: I have divided my data (a matrix of proximities expressed by cosines between 94 objects) into clusters with Ward hierarchical method, and I am very happy with the results from a visual point of view. However, I'd like to check whether the quality of the clusters is good. I've been told that for Ward method, statistics that can be good to check how good are the clusters are:   * cophenetic correlation coefficient,   * normalized mean square error,   * normalized mean absolute error. However, I don't know how to read them. I have a very low correlation coeff. (0.56) which sounds bad, since good correlations are averagely above .80. Then I have a 351.14 NMSE and a 17.24 NMAE. Is that bad? Is that good? How do I tell?


title: Examples of visualization featuring comparisons between countries?

text: Other then gapminder.org, could you please direct me to any good examples (interactive or static) of visualizations that compare between countries? Thanks.


title: Comparing an LME to a LM

text: I am finishing up some analysis on a large set of data. I would like to take the linear model used in the first part of the work and re-fit it using an LME. The LME would be very similar with the exception that one of the variables used in the model would be used in the random effects. This data comes from many observations (>1000) in a small group of subjects (~10) and I know that modeling the effect of subject is better done as a random effect (this is a variable that I want to shift). The R code would look like:               my_modelB <- lm(formula = A ~ B + C + D)         lme_model <- lme(fixed=A ~ B + C, random=~1|D, data=my_data, method='REML')      Everything runs fine and the results are vastly similar. It would be nice if I could use something like RLRsim or an AIC/BIC to compare these two models and decide which is the most appropriate. My colleagues don't want to report the LME because there isn't an easily accessible way of choosing which is "better", even though I think the LME is the more appropriate model. Any suggestions?


title: Cross-correlation in matlab: symmetrical plot and sign convention

text: I have two cosine signals in MATLAB both having a frequency of 100 Hz and the second signal has been delayed in time by 0.003 seconds. Both signals have been sampled at 1 Khz. When I use the cross correlation function and plot it to verify the delay between the signals, I expect to see a single point in time at which the signals are aligned but the plot looks symmetrical about the y axis. Why is that? Also I am not sure about the sign convention here. Does a negative x axis value at the cross correlation peak denote time lag or does a positive x axis value denote time lag?               st=1/1000;     t=[0:999]*st;          s1=cos(2*pi*100*t);     t1=t-(3*st);     s2=cos(2*pi*100*t1);          [c,lags]=xcorr(s1,s2);     plot(lags,c/max(c))


title: Constrained and Weighted Regression

text: I'd like to modify the answer to this question to allow weighted observations. I think all I need to do is weight the inputs X and Y.               X = w * X     Y = w * Y      the other parts of the procedure should follow. Please correct me if I am wrong.


title: GLM, Regression analysis, or other

text: I need to resolve a question with the use of SPSS, but don't really know what model to use. The question talks about the environment (scale variable (from not at all concerned - very concerned)), I need to compare it to different scale variables like gender (male/female), age group, education level (no formal education -...) and tell which groups are more/less concerned about the environment (without taking into account effects that are caused by other variables). How should I go about answering this question? I'm pretty sure that I would need to use the GLM model.


title: Correction for interaction model as it is significant in MR

text: I have used 1 DV and 5 IVs. I performed MR to test a model, entered age and gender in first block, IVs in second block, and interaction terms in third block. All the models are coming significant. How can I correct for the significance of interaction model as it is showing significant effect on DV now. And by the way, will it still be called hierarchical multiple regression, since I am not controlling for any moderator variable (except age and gender which are not actually moderators in my study)? That means can I add interaction terms in the same block or they have to be added in another block?


title: Question on Precision and Recall

text: Is it possible to have both Precision and Recall to be equal to 1. If such a case may exist, does it mean the system is over-fitting the dataset?


title: Cross-correlation with Spearman's method

text: I should investigate the lead-lag relationship among 2 sets of data, annual incidence of a disease and annual average data of an environmental factor for 13 consecutive years. Both sets are not normally distributed and therefore, correlation should be assessed by Spearman's correlation. But "cross- correlation analysis" in SPSS 16, checks correlation in lags by Pearson's method. SPSS does not have any option to choose the method for cross- correlation. what can I do to obtain Spearman's rho for lags?


title: Calculating return from a factor model

text: I've been trying to figure this one for weeks so your advise is very appreciated. I have log returns for the past 5 yrs on the following indexes: NYSE, Communication, Growth and Large Cap. I performed PCA on those and kept 2 ‘unobservable factors’ with loadings as follows:               NYSE Index  0.2599870   0.6090115     Communications  0.7952713   -0.5660120     Growth          0.3918862   0.2710737     Large Cap   0.3825905   0.4850304      Orthogonal regression of NYSE returns on 2 factor scores gives me 2 betas equal to factor loadings (0.609 and 0.2599). Therefore betas for NYSE are:                           NYSE       Index Communications Growth  Large Cap     NYSE Index  0.316      0.337                0.309    0.359      respectively To calculate delta R on NYSE in terms of the factors:               Return on NYSE= 0.316* R NYSE + 0.337* R Comm+0.309* R Growth +0.359 R Large Cap.      However this doesn’t actually work… What am I missing? How can I calculate Return on NYSE Index using PCA's common factors?


title: Sample Size Determination for Rare Event

text: Suppose a survey contained a dichotomous variable of chief interest, such as smoker/non-smoker, along with ten other dichotomous variables. Then suppose that the incidence of smokers is a rare event with an unknown population proportion. The goal is to infer differences between smokers and non-smokers from the answers given to the other ten variables e.g. 70% of smokers are male. How would one determine the appropriate sample size (say given 95% confidence and +-5% interval)? I'm assuming $N=Z^2 \frac{P(1-P)}{D^2}$ would not be appropriate.


title: How to determine percentiles/quantiles in SPSS?

text: I am really new to statistics and now I have some data and I want to make a frequency table with the yearly salary (a variable in my data). In the end, I need to be able to say what is the maximum salary of the first 25% people and the minimum salary of the top 10%. How should I do it in SPSS?


title: Constructing 95% confidence interval based on profile likelihood

text: In my elementary statistics course, I learnt how to construct 95% confidence interval such as population mean, mu based on **asymptotic normality** for "large" sample size. Apart from **resampling methods** , there is another approach based on **"profile likelihood"**. Could someone elucidate this approach? Under what situations, the constructed 95% CI based on asymptotic normality and profile likelihood are comparable? I could not find any references on this topic, any suggested references, please? Why isn't it more widely used?


title: How would you calculate average time to make a breath?

text: Suppose I wanted to know a one minute average of the time it takes to make a breath. So basically I would record the number of times I made a breath during one minute (one inhale and one exhale) and divide by 60 seconds to get the 1 minute average time to make a single breath? Would this be a distribution with a standard deviation?


title: Market mix modelling with R

text: Can somebody give me references (book/online resource) on using R for Marketing Mix Modelling?


title: Missing values on clinical scoring system in SPSS

text: I have a SPSS database with 747 patients with sickness symptoms and lab results as variables. With those symptoms and results I have generate two scoring systems. I sum the scores of all the variables to come to the final score. However, on one variable I have 546 missing values. For the 201 other patients I do have this value. How can I best deal with these missing values? Can I use sampling to generate a score for these missing values? I need a score on this variable to generate the final score. Does anyone have any suggestions? Big thanks!


title: Systems for symbolically solving expectations?

text: Do there exist any systems for symbolically solving expectations? This is sort of a follow-up to my question List of Tricks for Solving Messy Expectations? Basically, I'm looking for ways to solve a messy expectation after I've exhausted all obvious routes. **EDIT: BACKGROUND** I'm trying to solve the following for $\alpha$ (constrained to be greater than 0) as a function of $\sigma_X^2$, $\sigma_Y^2$, and $p$ $E\left[\ln(F) F^\alpha X^2 (X + Y)^2 + \ln(F) F^{2\alpha}(X + Y)^4\right] = 0$ where: $Y \sim N(0,\sigma_Y^2)$ $ X = \left\\{ \begin{array}{cc} N(0,\sigma_X^2) & p \\\ 0 & (1-p) \end{array}\right.$ where it's assumed $\sigma_X^2 \gg \sigma_Y^2$ and that $p$ is very small. (i.e., $X$ is a jump process with most of its mass at 0) $F = F_{|Z|}(|X+Y|)$ where $Z \sim N(0,\sigma_Y^2)$ (i.e., $F$ is the CDF for the absolute value of a normal) **EDIT: EVEN MORE BACKGROUND** The equation I'm trying to solve above is the first order condition for the min-MSE problem: $\min_{\alpha > 0} \left(X^2 - \widehat{X^2}\right)^2 $ where $ \widehat{X^2} = F_{|Z|}\left(|R|\right)^{\alpha}R^2$ and $R = X + Y\,$ is the only observed variable. Basically, I'm trying to estimate the square of the jump, $X^2$ (given that I can only observe the aggregated process $R$) by smoothing down $R^2$. If $|R|$ is large, the smoothing function $F_{|Z|}\left(|R|\right)^{\alpha}$ should be close to 1 and the estimate of $X^2$ would be close to $R^2$. If $|R|$ is small, the estimate of $X^2$ would be close to 0.


title: What methods exist for finding optimal splits to discretize continuous data with respect to a target variable

text: I'm doing some research into methods for discretizing a continuous variable coupled with a binary target variable to find the optimal split points to maxamise a measure of impurity (gini/entropy). First off I'm having some trouble coming up with google-able terms. "Optimal Split" seems to relate to choosing which variable to split on in a decision tree as opposed to how to split that variable up. Is there a defined title on this problem? I imagine to convert a continuous variable into a binary variable I could potentially try setting a split point at each distinct value in the data and take the one that returns the maximum entropy or gini. But to extend it beyond binary splits the search space would grow to make that method pretty expensive. Are there any popular methods for solving this?


title: Defining precision in a mechanical movement

text: I have a device which rotates on a stepper motor through a belt and pulley system. I would like to know both the positional accuracy and precision of this movement. It is of my opinion that this should be determined by taking a large number of measurements to random angles to give a bell curve then using the offset of the mean from 0 error to define the positional accuracy and defining the precision of the movement as 3 sigma. The rest of the office are of the opinion that accuracy is the range of errors when rotating to random points and precision is the range of errors when rotating to a determined point, away to a random point, then back to the determined point. Is this just me applying my past experience in manufacturing to a problem that doesnt require it or is there merit in my method and everyone else is just specifying range, which is not what they really want?


title: Back propagation in Convolutional neural networks

text: I am trying to understand how CNN works. I want to use them in object recognition task. I thouhgt that CNN is unsupervised networks. My main question is how can I implement the back propagation without using image targets (pre-defined labels), in order to calculate the error function which will affect the weights of the network.


title: What does "X"stand for in "XARMAX"?

text: What does "X" stand for in "XARMAX"? I searched online but couldn't find it. Is there "XARMA"? Thanks!


title: What kind of statistical test could help show this?

text: I have a log of customer call data, and would like to know what kind of statistical test I could use to test for statistical significance. Customers call with a unique call reason, and fit into either one customer type, young customers or elderly customers.               Call Reason      Young customers      Elderly Customers     Bill Enquiry     5                    28     Product Problem  15                   20     Product Upgrade  25                   12      The numbers represent the number of calls for each call reason from each customer type. I want to know how I can test if there is statistically significant difference (with 95% confidence, say) in the type of customers that call for each call reason. For example, looking at the data, it appears that many more elderly customers call with billing enquiries than young customers. How might I rigorously test this?


title: How to show that $\mathrm{mgf}$ $M(s)$ and $\mathrm{pgf}$ $P(s)$ are related?

text: Let $X$ be an integer-valued $rv$ with $\mathrm{pgf}$ $P(s)$ (probability generating functions) and suppose that $\mathrm{mgf}$ $M(s)$ (moment generating functions) exist for $s∈(-s_0,s_0),s_0>0$. How can we proof that $M(s)$ and $P(s)$ are related??


title: What is the name for a regression undertaken with robust variance estimators for clustered data?

text: Just wondering what you call a regression undertaken with robust variance estimators for clustered data. Is it ok to say you have used OLS regression with robust variance estimators for clustered data, or is this fundamentally not the case (i.e. once you are using robust variance estimators for clustered data you are no longer doing OLS)?. Reference would be great.


title: When doing a t-test for the significance of a regression coefficient, why is the degrees of freedom n - predictors - 1?

text: I read at http://www.unesco.org/webworld/idams/advguide/Chapt5_2.htm that this was the degrees of freedom I should use, but I don't understand why. My understanding was that t-tests generally had n-1 degrees of freedom.


title: What do statisticians do that can't be automated?

text: Will software eventually make statisticians obsolete? What is done that can't be programmed into a computer?


title: Calculation of seasonal (annual) component of time-series: use of cross-validation?

text: I've been working for almost a year on electricity load forecasting in collaboration with some climate scientists, using temperature data obtained from models. Instead of using directly temperature values, we use the 'anomaly' from the annual average (called 'climatology'). It seems, in climate sciences, that is good practice to use 'cross-validation' to calculate this anomaly: the anomalies of year X are calculated considering the climatology of all the other years except X, a kind of 'K-FOLD' method with K equals to the number of years considered in the whole calculation. I'm a computer scientist and we are writing for an engineering journal and I was wondering if this whole procedure to calculate anomalies with this kind of 'K-FOLD' it is really considered necessary to make our entire work more robust. I've never read of similar procedure in other forecasting or engineering papers and for this reason I'd like to avoid unnecessary work (and code writing) at this stage of the work. Thank you for your comments,


title: How to obtain optimally scaled data from homals in R?

text: I am using the `homals`-package in R to optimally quantify my data that contains categorical (nominal and ordinal) as well as numerical variables. My goal is to quantify the categorical variables using optimal scaling. To achieve this, I use the rank constraint ` rank = 1`, since this "reduces homogeneity analysis to principal component analysis (with optimal scaling of the variables)". In order to preserve all the information contained in the data I set `ndim` to the number of variables in the data. The categorical variables contain NAs. How do I obtain the optimally quantified data (rather than the dimensions/principal components) from the results of `homals()`? Since I am not looking for the principal components, I am not looking for `objscores`. I am almost able to extract the desired quantifications, as via `low.rank` I know exactly how each category is transformed. However, these transformations do not include the NAs in the data and I would like to know too how these missing observations were quantified. From a variety of documentation (Package 'homals' (the package documentation), Homogeneity Analysis in R: The Package homals, and HOMALS) I have not been able to figure out how to do this. Is anybody familiar with the `homals`-package that knows how to do this or does anybody know how to compute the quantified data from the results?


title: Common standards for data import file formats for statistics software?

text: I'm looking for common standards (file formats) for how most statistics software packages import data. The reason is that I'm going to make an export function for a simulation software package. The data to be exported consists of time-series of (1 or more) values for all members of a population. The population size can be in the 100,000's. The length of the time-series is typically in the 1,000's. The number of values per member per timestep can range from 1 to 20. How is such data typically imported by statistics packages?


title: KNN imputation R packages

text: I am looking for a KNN imputation package. I have been looking at imputation package (http://cran.r-project.org/web/packages/imputation/imputation.pdf) but for some reason the KNN impute function (even when following the example from the description) only seems to impute zero values (as per below). I have been looking around but cannot find something yet, and hence was wondering if anyone has other suggestions for good KNN imputation packages? W In the code per below - the NA values are replaced by zero's - not by the Knn mean value               require(imputation)     x = matrix(rnorm(100),10,10)     x.missing = x > 1     x[x.missing] = NA     kNNImpute(x, 3)     x


title: SAS fails to fit a mixed model

text:
Here are some simulated data:                   library(mvtnorm)         I <- 3 # positions (fixed factor)         J <- 4 # tubes (random factor)         K <- 4 # repeats          n <- I*J*K         set.seed(123)         tube <- rep(1:J, each=I)         position <- rep(LETTERS[1:I], times=J)         Mu_i <- 3*(1:I)         Mu_ij <- c(t(rmvnorm(J, mean=Mu_i)) )           tube <- rep(tube, each=K)         position <- rep(position, each=K)         Mu_ij <- rep(Mu_ij, each=K)         dat <- data.frame(tube, position, Mu_ij)         sigmaw <- 2         dat$y <- rnorm(n, dat$Mu_ij, sigmaw)         dat$tube <- factor(dat$tube)          > str(dat)     'data.frame':   48 obs. of  4 variables:      $ tube    : Factor w/ 4 levels "1","2","3","4": 1 1 1 1 1 1 1 1 1 1 ...
          $ position: Factor w/ 3 levels "A","B","C": 1 1 1 1 2 2 2 2 3 3 ...      $ Mu_ij   : num  2.44 2.44 2.44 2.44 6.13 ...
          $ y       : num  3.24 2.66 1.33 6.01 7.12 ...     > head(dat)       tube position    Mu_ij        y     1    1        A 2.439524 3.241067     2    1        A 2.439524 2.660890     3    1        A 2.439524 1.327842     4    1        A 2.439524 6.013351     5    1        B 6.129288 7.124989     6    1        B 6.129288 2.196053      I fit a mixed model with R, it works well:               > library(lme4)     > lmer(y ~ position + (0+position|tube), data=dat)     Linear mixed model fit by REML      Formula: y ~ position + (0 + position | tube)         Data: dat         AIC   BIC logLik deviance REMLdev      212.6 231.3  -96.3    194.8   192.6     Random effects:      Groups   Name      Variance Std.Dev. Corr                tube     positionA 0.30123  0.54885                               positionB 0.68317  0.82654  -0.695                       positionC 1.66666  1.29099  -0.408  0.940       Residual           3.14003  1.77201                     Number of obs: 48, groups: tube, 4          Fixed effects:                 Estimate Std. Error t value     (Intercept)   3.3533     0.5211   6.435     positionB     3.1098     0.8923   3.485     positionC     5.6138     1.0144   5.534          Correlation of Fixed Effects:               (Intr) postnB     positionB -0.753            positionC -0.651  0.744      But the same model does not work well with SAS:               PROC MIXED DATA=dat ;     CLASS POSITION TUBE ;     MODEL y = POSITION ;     RANDOM POSITION / subject=TUBE type=UN G GCORR ;     RUN; QUIT;      gives                Estimated G matrix is not positive definite.                              Estimated G Matrix                         Row    Effect      position    tube        Col1        Col2        Col3                           1    position    A           1        0.08895     -0.5823     -0.1545                      2    position    B           1        -0.5823      0.1455      1.2431                      3    position    C           1        -0.1545      1.2431      1.4835      Is it possible to remedy this failure ?


title: Mediator versus moderators - not sure which is more appropriate

text: I am undertaking a study which has the following aims: The main aim of the current study is to investigate the effects of food messages about the healthiness of foods on actual food intake. Additionally, a secondary aim is to discover the extent to which impulsivity/sensation seeking, taste preference and BMI may moderate/mediate the effects of food-related messages on food intake. I am not sure which is more appropriate in this case to mediate or moderate. My DV is food consumption and IV's are food messages, impulsivity/sensation seeking, taste preference and BMI. It is an experimental study which involves three treatment groups which will receive either a healthy message, unhealthy message or no message about the cookie. They will be asked to taste and rate and consume 'cookies'. The taste rating will also ask questions about taste preference. I will then provide a questionnaire which gathers info about impulsivity/sensation seeking and also BMI. I am leaning towards mediation but am not sure. I believe that those in the healthy message group will consume more cookies based on prior research but I also deem that the other variables will have a influence as well on cookie consumption. Basically those who have high impulsivity, taste preference for the cookie and higher BMI either individually or in combination will consume more cookies. So I am not sure as I've indicated whether it would be more suitable to have moderator's or mediators. I am somehow leaning towards mediation but would like some clarity around this. Also if it's either mediation or moderation which IV's are which? Also what would a hypothesis look like? I am guessing that if I went with moderation for example then a hypothesis may look something like: the relation between health message and cookie consumption could be stronger for those who score higher on the impSS as compared to those who don’t score as high But I may be wrong. Any help would be appreciated. Thanks. MO


title: Include two correlated factors in mixed model?

text: In my experiment in which participants make a yes-no decision, I use two measures of empathy because I predict these to be both of influence. The two questionnaire scales are weakly correlated but do not measure exactly the same. Can I include these two scales in my mixed model ANOVA or will they sort of average each other out and should I, if I want to use both measures of empathy run two separate models?               MIXED decision BY treatment target WITH empathic_concern perspective_taking      /CRITERIA=CIN(95) MXITER(100) MXSTEP(10) SCORING(1) SINGULAR(0.000000000001) HCONVERGE(0, ABSOLUTE) LCONVERGE(0, ABSOLUTE) PCONVERGE(0.000001, ABSOLUTE)        /FIXED=treatment target empathic_concern perspective_taking     treatment * target perspective_taking     treatment * empathic_concern     treatment * perspective_taking     empathic_concern * perspective_taking     treatment * perspective_taking * empathic_concern       | SSTYPE(3)        /METHOD=ML        /PRINT=CORB DESCRIPTIVES  SOLUTION        /RANDOM=INTERCEPT | SUBJECT(id) COVTYPE(VC)        /EMMEANS=TABLES(OVERALL).


title: Do data and residuals of a VAR model have to be of normal distribution?

text: Does vector autoregression (VAR) model require data to be of normal distribution? What are the pitfalls if the residuals are not of normal distribution?


title: Multiple comparison issue in multivariate logistic regression

text: I am assessing a bunch of risk factors related to HIV infection by logistic regression. First I will conduct a couple of chi square tests or t-tests to see if each factor is associated with HIV infection or not. Then I will pick those factors with p value less than 0.1 out of the bivariate analysis above and include them in the multivariable logistic regression model (adjusted model). In this case, is multiple comparison a concern for me? Why and why not?


title: What is the difference between a hazard ratio and the e^coef of a Cox equation?

text: As I understand it, the definition of the hazard ratio is the ratio of two hazard rates. Often the exp(coef) from a Cox model is also used as an estimate of the hazard ratio. These methods give two different, although similar, results. Are there circumstances in which it's appropriate to use one method vs. the other?


title: How to analyse periodic time series volatility?

text: I am wondering what is the best way to extract knowledge about periodic time series data? In my case, I am trying to analyse historical hourly electricity prices to get information about its volatility, when the peaks occur, what is its distribution among the hours of the day or something similar to be able to make good decisions on when to throttle devices to save money and stay energy efficient. My original idea was to simply aggregate the time series by calculating mean values for each hour over a period of several months and then sort them to pick the top M% hours and say that they are the most expensive on the average and should be avoided. This is a very crude approach, so I'm wondering whether I'm making some obvious mistake, disregarding some information? Is there some better approach? It seems to me that the information on how high the current prices are in comparison to the average should also be used to sometimes exclude, say, 3 hours from a day and on other times only 1 hour.


title: Nonlinear regression with genetic programming for data from full-factorial design

text: I am trying to do nonlinear regression on a dataset obtained from full factorial design (2 and 3 levels) with the help of genetic programming. Is there something fundamentally wrong with trying to construct a nonlinear regression model with full-factorial design? My background is not in statistics, and I did read somewhere that full factorial design (esp when few levels are used) are supposed to be used at most to identify contribution of interaction factors, and key factors, and not for full blown nonlinear regression. Would appreciate your answers...


title: Best process for drawing numbers

text: This is a follow up to my previous question. Assume I want to randomly draw two numbers from the numbers 1 to 8, one even number and one odd number. I would create two urns, one with the even, one with the odd numbers, and draw one number from each. Thus, all numbers have the same probability to be drawn. This was the accepted answer to my previous question. Now assume I have the numbers from 1 to 9. I want to draw one number that is smaller or equal to 5, i.e. {1,2,3,4,5}, and one number that is larger or equal to 5, i.e. {5,6,7,8,9}. How would I draw them so that all numbers have the same probability to be drawn? My idea would be to create two urns, one with the numbers {1, 1, 2, 2, 3, 3, 4, 4, 5}, the other with the numbers {5, 6, 6, 7, 7, 8, 8, 9, 9}. If I flip a coin to randomly chose one urn, then randomly draw one number, all numbers have the same probability to be drawn. But if I draw a 5 and then switch to the other urn, there is a possibility that I draw a second 5 -- which is not what I want to allow to happen. How would you do it? _If you know a more meaningful title for this question, please edit it._


title: Occam's razor obsolete?

text: I saw Vapnik's books about statistical learning... I read the first few chapters. Anyway what surprised me the most was that he thought that the Occam's razor was obsolete. I thought it was related to the situation in which assuming higher dimension improves the fit significantly. Am I understood right? Is it correct that Occam's razor can't be no longer right like Vapnik said? Is there any argument that Occam's razor shouldn't be considered as the default? The exact sentences come from the Preface to the Second Edition of The Nature of Statistical Learning, which are: > The years since the first edition of the book have also changed the general > philosophy in our understanding the of nature of the induction problem. > After many successful experiments with SVM, researchers became more > determined in criticism of the classical philosophy of generalization based > on the principle of Occam's razor." I wish someone could elaborate on the criticism on Occam's razor.


title: Is the symmetry of u function for  the robust M estimation mandatory?

text: I have build a $\rho$ function which has the following definition: \begin{equation} \rho(x)= \left\\{ \begin{array}{ll} 4- \frac{8}{x^2} \text{if } x \lt-2\\\ \frac{x^2}{2} \text{if } x \in [-2,3] \\\ 9-\frac{81}{2*x^2} \text{if } x \gt 3 \end{array} \right. \end{equation} The function has the continuity and derivability properties, but is asymetric. Can be solved using Iteratively reweighted least squares (IRLS) method? The $u$ function is: $ u(x)=\begin{equation} \left\\{ \begin{array}{ll} \frac{16}{x^4} \text{ if } x \lt-2\\\ 1 \text{ if } x \in [-2,3] \\\ \frac{81}{x^4} \text{ if } x \gt 3 \end{array} \right. \end{equation} $ For the IRLS method, the symmetry of $u$ function is mandatory?


title: Comparing coefficients in logistic regression

text: I have some problems I need help with. I am running a binary logistic regression.               DV: Brand choice (0/1)     IV1: Attitude towards product (p<0.05)     IV2: Price sensitivity (p<0.05)      (Both IV1 and IV2 are measured on the same scale) I have found that “Attitude towards product” is more powerful ($β=1.308$) than “Price sensitivity” ($β=0.956$) in predicting brand choice. However, according to my thesis supervisor, in order to claim this fact, I need to run an additional test. That is, I need to find out if this difference is significant. So, what I want to know is if: $β_{Attitude}> β_{Price}$ (or if $β_{Attitude}≠ β_{Price}$) and if this is statistically significant. I have been searching a lot, but I cannot find how I am supposed to test this. Does anybody have any ideas?? I have read that it is possible to run a t-test (with brand choice as grouping variable) and look at the t-values and if the t-value is higher for “Attitude towards product” is higher, then it is the stronger one. However, I have 3 groups (conditions) so I’m not sure that I can run a t-test. Can I use an ANOVA instead and look at the F-values? I also came across some information that I could include an interaction-term (Attitude*Price) and that the p-value for each interaction term gives me a significance test for the difference in those coefficients. Is this a valid method??


title: SVM model selection for datasets with sharp corners

text: I'm working with an artificially generated dataset that is separated by many sharp corners. As an example, imagine an H-shape in a 3D (or higher- dimensional) space. Points within the H are positive, points outside are negative. In my application, many new points will be generated at runtime. I would like to classify them as being inside the H-shape or not. I currently do an exact check for these points using standard collision detection methods, but since this is proving expensive I would like to do much faster classification at the expense of some accuracy. For instance, 90-95% correct classifications at an order of magnitude or more faster speed than exact checking is ideal for me. I'm using an SVM with the radial basis function (RBF) kernel since my dataset is not linearly separable in the original space. I've trained it using grid- search on the parameters gamma and nu, with 3-fold cross validation, varying both the range and granularity of the search. However, the best I can do is about 80-85% accuracy, and a running time of 6-10s for 100,000 new datapoints, which is about the same time as it takes to do the exact check. Where am I going wrong? I particularly suspect that the RBF is poor at handling sharp corners (e.g. the H-shape that I described, whose boundary separates positive and negative points). Should I be using a different kernel that is more suited to these parameters? Should I use some combination of kernels? Should I do some form of region decomposition, and train separate kernels for each region? Should I be using a different ML approach (e.g. neural networks) instead?


title: How to test whether button pressing behaviour is random or is in response to a stimulus?

text: I have a set of univariate data, where the response variable in a trial is a sequence of time intervals between button presses (my experimental subjs have to press a button intermittently to react to some task conditions). Eventually I have to run a within-subject regression to see if my hypothesized predictors are indeed significant. However, the first thing I want to do is to show that the subjects didn't just press the button randomly. In other words, I need to make sure that the person didn't just press a button at a regular pace, (averaged response times + random variability). One idea I have is to check if my data resemble a Poisson process, in which arrival times occur randomly at a certain rate. Based on what I've looked up on using Poisson process so far, there are two approaches to fit the poisson process to my time interval data: 1) count how many times the person pressed a button in each trial of fixed time period. Then I can follow the Poisson distribution fitting process to count data.   2) the wait times between events in a poisson process is said to follow an exponential distribution. So I can try to fit "exponential" distribution to the time intervals. In that case, should I try to fit a bunch of different distributions such as Exponential, Gamma, and Weibull? In either case, I would then use a chi-square goodness of fit test and/or a QQ-plot to assess whether the distribution fit can be rejected. Any advice would be greatly appreciated. Is fitting to Poisson Process to my data even the right way to go? I would be happy to hear about any other suggestions. Thanks in advance.


title: How do I fix this error in mice multiple imputation?

text: I am trying to do multiple imputation using the mice package in R, and the imputation keeps stopping with the following error:               Error in mice.impute.logreg(c(1L, 1L, 2L, 1L, 1L, 1L, 2L, 1L, 1L, 2L,  :        dims [product 145] do not match the length of object [146]     In addition: There were 50 or more warnings (use warnings() to see the first 50)     > warnings()     Warning messages:     1: In rchisq(1, sum(ry) - ncol(x)) : NAs produced     2: In runif(length(d), 0, a1/10^10) : NAs produced     3: In runif(length(d), 0, a1/10^10) : NAs produced     4: In runif(length(d), 0, a1/10^10) : NAs produced     5: In runif(length(d), 0, a1/10^10) : NAs produced     6: In runif(length(d), 0, a1/10^10) : NAs produced     7: In runif(length(d), 0, a1/10^10) : NAs produced     8: In rchisq(1, sum(ry) - ncol(x)) : NAs produced     9: In runif(length(d), 0, a1/10^10) : NAs produced     10: In runif(length(d), 0, a1/10^10) : NAs produced     11: In runif(length(d), 0, a1/10^10) : NAs produced     12: In runif(length(d), 0, a1/10^10) : NAs produced     13: In runif(length(d), 0, a1/10^10) : NAs produced     14: In runif(length(d), 0, a1/10^10) : NAs produced     15: In runif(length(d), 0, a1/10^10) : NAs produced     16: In runif(length(d), 0, a1/10^10) : NAs produced     17: In runif(length(d), 0, a1/10^10) : NAs produced     18: In runif(length(d), 0, a1/10^10) : NAs produced     19: In rchisq(1, sum(ry) - ncol(x)) : NAs produced     20: In runif(length(d), 0, a1/10^10) : NAs produced     21: In runif(length(d), 0, a1/10^10) : NAs produced     22: In runif(length(d), 0, a1/10^10) : NAs produced     23: In runif(length(d), 0, a1/10^10) : NAs produced     24: In runif(length(d), 0, a1/10^10) : NAs produced     25: In runif(length(d), 0, a1/10^10) : NAs produced     26: In runif(length(d), 0, a1/10^10) : NAs produced     27: In runif(length(d), 0, a1/10^10) : NAs produced     28: In runif(length(d), 0, a1/10^10) : NAs produced     29: In runif(length(d), 0, a1/10^10) : NAs produced     30: In rchisq(1, sum(ry) - ncol(x)) : NAs produced     31: In runif(length(d), 0, a1/10^10) : NAs produced     32: In runif(length(d), 0, a1/10^10) : NAs produced     33: In runif(length(d), 0, a1/10^10) : NAs produced     34: In runif(length(d), 0, a1/10^10) : NAs produced     35: In runif(length(d), 0, a1/10^10) : NAs produced     36: In rchisq(1, sum(ry) - ncol(x)) : NAs produced     37: In runif(length(d), 0, a1/10^10) : NAs produced     38: In runif(length(d), 0, a1/10^10) : NAs produced     39: In runif(length(d), 0, a1/10^10) : NAs produced     40: In runif(length(d), 0, a1/10^10) : NAs produced     41: In rchisq(1, sum(ry) - ncol(x)) : NAs produced     42: In runif(length(d), 0, a1/10^10) : NAs produced     43: In runif(length(d), 0, a1/10^10) : NAs produced     44: In runif(length(d), 0, a1/10^10) : NAs produced     45: In runif(length(d), 0, a1/10^10) : NAs produced     46: In rchisq(1, sum(ry) - ncol(x)) : NAs produced     47: In runif(length(d), 0, a1/10^10) : NAs produced     48: In runif(length(d), 0, a1/10^10) : NAs produced     49: In runif(length(d), 0, a1/10^10) : NAs produced     50: In runif(length(d), 0, a1/10^10) : NAs produced      If the seed is set to a particular value, the error occurs at the same point, otherwise the error will occur at different iterations, imputations and variables in different runs. OK. So here is the data structure. This was imported from Stata using read.dta. The list of variables is truncated and does not show the ones that are showing the problem. They are all binary, with a single missing value each, in the same record.               > str(peimp.in)     'data.frame':   146 obs. of  164 variables:      $ id                        : chr  "195" "128" "218" "1106" ...          $ distress0                 : num  5 6 0 2 1 0 5 1 1 5 ...      $ pcs0                      : num  47.1 46.9 60.8 44 61.5 ...          $ mcs0                      : num  33.4 43.8 52.4 38.6 53.1 ...      $ pwb0                      : num  24 21 28 17 27 NA 8 16 24 8 ...          $ swb0                      : num  24 25 28 24 28 NA 17 28 28 24 ...      $ ewb0                      : num  18 14 20 15 20 ...          $ fwb0                      : num  14 21 28 18.7 26 ...      $ ccs0                      : num  14 22 25 25 26 NA 16 9 23 14 ...          $ pn.sf70                   : int  2 2 0 4 0 3 1 3 3 4 ...      $ pn.sf80                   : int  2 1 0 2 0 2 3 3 1 3 ...          $ pn.fcgp40                 : int  1 1 0 4 0 NA 1 3 3 3 ...      $ pn.aqol0                  : int  1 1 0 1 0 1 0 1 1 2 ...          $ pf.nbs0                   : num  51.8 55.6 55.6 40.3 57.5 ...      $ rp.nbs0                   : num  21.2 45.9 54.9 34.7 57.2 ...          $ bp.nbs0                   : num  42.6 46.7 62 34.6 62 ...      $ gh.nbs0                   : num  54.6 28.5 66.5 61.7 62.7 ...          $ vt.nbs0                   : num  49.6 49.6 70.4 40.7 58.5 ...      $ sf.nbs0                   : num  37.3 47.3 37.3 42.3 57.3 ...          $ re.nbs0                   : num  24.8 45.7 56.2 45.7 56.2 ...      $ mh.nbs0                   : num  37.8 43 53.5 28 50.9 ...          $ distress1                 : int  8 7 7 5 0 5 NA 8 0 2 ...      $ pcs1                      : num  36.4 23.3 43.9 27.5 33.9 ...          $ mcs1                      : num  32.1 32.2 45.2 32.5 51 ...      $ pwb1                      : num  17 3 13 13 25 19 14 0 24 21 ...          $ swb1                      : num  24 24.5 27 26.8 28 ...      $ ewb1                      : num  17 13 4 18 18 22 20 17 24 22 ...          $ fwb1                      : num  15 3.5 18.7 6 15 ...      $ ccs1                      : num  13 16 21 11 24 ...          $ pn.sf71                   : int  2 4 4 3 2 3 5 4 1 3 ...      $ pn.sf81                   : int  3 4 3 4 1 3 3 4 2 3 ...          $ pn.fcgp41                 : int  1 4 3 1 1 2 2 4 0 3 ...      $ pn.aqol1                  : int  1 1 2 1 1 0 1 2 0 1 ...          $ pf.nbs1                   : num  40.3 19.3 38.4 19.3 42.2 ...      $ rp.nbs1                   : num  28 21.2 57.2 25.7 21.2 ...          $ bp.nbs1                   : num  38.6 26.5 30.5 30.1 46.7 ...      $ gh.nbs1                   : num  34.2 28.5 59.4 46 42.7 ...          $ vt.nbs1                   : num  31.8 38.7 46.7 25.9 37.7 ...      $ sf.nbs1                   : num  27.3 22.2 22.2 17.2 47.3 ...          $ re.nbs1                   : num  31.8 21.4 56.2 31.8 49.2 ...      $ mh.nbs1                   : num  37.8 32.6 44.3 35.2 48.2 ...          $ distress2                 : int  7 7 NA 8 0 0 10 1 0 0 ...      $ pcs2                      : num  24 29.9 NA 37.7 51.7 ...          $ mcs2                      : num  35.6 32.5 NA 29.8 59.9 ...      $ pwb2                      : num  12 9 NA 11 28 22 17 20 26 24 ...          $ swb2                      : num  24 23.3 NA 8 28 ...      $ ewb2                      : int  22 13 NA 12 19 24 10 18 22 23 ...          $ fwb2                      : num  8 13 NA 2 28 10.5 7 4 22 18 ...      $ ccs2                      : num  10 10 NA 12 26 21 14 12 27 23 ...          $ pn.sf72                   : int  5 4 NA 2 0 3 3 3 1 1 ...      $ pn.sf82                   : int  4 3 NA 2 0 NA 4 3 1 0 ...          $ pn.fcgp42                 : int  4 4 NA 1 0 1 3 2 0 1 ...      $ pn.aqol2                  : int  3 2 NA 1 0 1 1 1 0 0 ...          $ pf.nbs2                   : num  21.2 28.8 NA 30.8 49.9 ...      $ rp.nbs2                   : num  21.2 30.2 NA 30.2 54.9 ...          $ bp.nbs2                   : num  21.7 30.5 NA 42.6 62 ...      $ gh.nbs2                   : num  38 26.1 NA 38 46 ...          $ vt.nbs2                   : num  43.7 37.7 NA 31.8 61.5 ...      $ sf.nbs2                   : num  22.2 27.3 NA 32.3 57.3 ...          $ re.nbs2                   : num  17.9 24.8 NA 31.8 56.2 ...      $ mh.nbs2                   : num  40.4 35.2 NA 27.3 58.7 ...          $ seegp2                    : Factor w/ 2 levels "No","Yes": NA 2 1 2 1 1 2 2 2 2 ...      $ admed2                    : Factor w/ 2 levels "No","Yes": NA 1 2 1 1 1 1 1 1 1 ...          $ admhsp2                   : Factor w/ 2 levels "No","Yes": NA 2 2 2 1 2 1 2 1 1 ...      $ hlp.paid2                 : Factor w/ 2 levels "No","Yes": NA 1 1 1 1 1 2 1 1 1 ...          $ hlp.unpaid2               : Factor w/ 2 levels "No","Yes": NA 1 1 1 1 1 1 1 1 1 ...      $ hlp.ff2                   : Factor w/ 2 levels "No","Yes": NA 1 1 1 1 1 2 2 2 1 ...          $ inhosp2                   : Factor w/ 2 levels "_","Unwell/in hospital": 1 1 1 1 1 1 1 1 1 1 ...      $ wthdrw2                   : Factor w/ 2 levels "_","Withdrawn": 1 1 1 1 1 1 1 1 1 1 ...          $ dth2                      : Factor w/ 2 levels "_","Deceased": 1 1 1 1 1 1 1 1 1 1 ...      $ distress3                 : num  3 4 0 5 0 4 10 0 0 1 ...          $ pcs3                      : num  42.1 24.9 46.9 32.3 55.3 ...      $ mcs3                      : num  43.1 44.2 69.2 30.3 58.5 ...          $ pwb3                      : num  17 10 28 10 28 ...      $ swb3                      : num  18 21 24 16 28 ...          $ ewb3                      : int  20 15 24 13 21 21 18 19 23 23 ...      $ fwb3                      : num  14 11 28 7 25 NA 7 11 28 16 ...          $ ccs3                      : num  20 13 28 10 25 ...      $ pn.sf73                   : int  3 4 0 3 0 2 4 1 1 3 ...          $ pn.sf83                   : int  2 4 0 3 0 0 3 0 0 2 ...      $ pn.fcgp43                 : int  2 4 0 3 0 1 4 1 0 2 ...          $ pn.aqol3                  : int  1 2 0 1 0 1 1 0 0 1 ...      $ pf.nbs3                   : num  46.1 26.9 26.9 32.7 55.6 ...          $ rp.nbs3                   : num  34.7 23.5 57.2 30.2 54.9 ...      $ bp.nbs3                   : num  38.2 26.5 62 34.2 62 ...          $ gh.nbs3                   : num  49.9 30.8 66.5 26.1 50.8 ...      $ vt.nbs3                   : num  37.7 46.7 70.4 31.8 61.5 ...          $ sf.nbs3                   : num  47.3 42.3 57.3 27.3 57.3 ...      $ re.nbs3                   : num  35.3 28.3 56.2 31.8 56.2 ...          $ mh.nbs3                   : num  48.2 40.4 64 29.9 58.7 ...      $ seegp3                    : Factor w/ 2 levels "No","Yes": 1 NA 1 2 2 2 1 2 2 2 ...          $ admed3                    : Factor w/ 2 levels "No","Yes": 1 NA 2 2 1 1 1 1 1 1 ...      $ admhsp3                   : Factor w/ 2 levels "No","Yes": 2 NA 1 2 1 2 2 1 1 1 ...          $ hlp.paid3                 : Factor w/ 2 levels "No","Yes": 1 NA 2 2 1 1 2 2 1 1 ...      $ hlp.unpaid3               : Factor w/ 2 levels "No","Yes": 1 NA 1 1 1 1 1 1 1 1 ...          $ hlp.ff3                   : Factor w/ 2 levels "No","Yes": 1 NA 2 1 1 1 2 2 1 1 ...      $ inhosp3                   : Factor w/ 2 levels "_","Unwell/in hospital": 1 1 1 1 1 2 1 1 1 1 ...          $ wthdrw3                   : Factor w/ 2 levels "_","Withdrawn": 1 1 1 1 1 1 1 1 1 1 ...      $ dth3                      : Factor w/ 2 levels "_","Deceased": 1 1 1 1 1 1 1 1 1 1 ...           [list output truncated]          - attr(*, "datalabel")= chr "Dataset for imputation"          - attr(*, "time.stamp")= chr "23 Jul 2012 09:19"          - attr(*, "formats")= chr  "%5s" "%9.0g" "%9.0g" "%9.0g" ...          - attr(*, "types")= int  5 254 254 254 254 254 254 254 254 251 ...          - attr(*, "val.labels")= chr  "" "" "" "" ...          - attr(*, "var.labels")= chr  "patient id" "0 distress" "0 pcs" "0 mcs" ...          - attr(*, "version")= int 12          - attr(*, "label.table")=List of 18           ..$ wthdrw  : Named num  0 1       .. ..- attr(*, "names")= chr  "_" "Withdrawn"       ..$ dth     : Named num  0 1           .. ..- attr(*, "names")= chr  "_" "Deceased"           ..$ inhosp  : Named num  0 1       .. ..- attr(*, "names")= chr  "_" "Unwell/in hospital"       ..$ noyes   : Named num  0 1           .. ..- attr(*, "names")= chr  "No" "Yes"           ..$ hosp    : Named num  0 1       .. ..- attr(*, "names")= chr  "Sydney" "Other"       ..$ TYPCANCE: Named num  1 2 3 4           .. ..- attr(*, "names")= chr  "primary rectal" "recurrent rectal" "primary other" "recurrent other"           ..$ TYPEEXEN: Named num  1 2 9       .. ..- attr(*, "names")= chr  "curative exenteration" "palliative exenteration" "not applicable"       ..$ GENDER  : Named num  1 2           .. ..- attr(*, "names")= chr  "male" "female"           ..$ MARITALS: Named num  1 2 3 4       .. ..- attr(*, "names")= chr  "single" "married / living with partner" "divorced" "widowed"       ..$ EMPLOYME: Named num  1 2 3 4 5           .. ..- attr(*, "names")= chr  "full time" "part time" "retired" "unemployed" ...           ..$ LABB    : Named num  1 2       .. ..- attr(*, "names")= chr  "yes" "no"       ..$ LABC    : Named num  1 2           .. ..- attr(*, "names")= chr  "yes" "no"           ..$ CLEAROPM: Named num  1 2       .. ..- attr(*, "names")= chr  "r0 / r1" "r2"       ..$ RECADPRE: Named num  1 2           .. ..- attr(*, "names")= chr  "recurrent" "advanced primary"           ..$ RECADNEO: Named num  1 2 3       .. ..- attr(*, "names")= chr  "yes, short course" "yes, long course" "no"       ..$ V132_A  : Named num  1 2           .. ..- attr(*, "names")= chr  "yes" "no"           ..$ RECADEXE: Named num  1 2 8       .. ..- attr(*, "names")= chr  "curative" "palliative" "no surgery"       ..$ OPASASCO: Named num  1 2 3 4 5       .. ..- attr(*, "names")= chr  "healthy patient" "mild systemic disease - no functional limitation" "severe systemic disease - definite functional limitation" "severe systemic disease - constant threat to life" ...


title: If model is significant but some coefficients are not

text: Suppose you have a linear model where the model is significant but some coefficients are not. How does one interpret the model when some coefficients are not significant?


title: soft question - Why and when is classification useful

text: From What I understand about classification is that it is used to distinguish unlabeled data points in a collection. So if we have data which is labeled(age, height, weight,...) then we don't need classification, its used only in cases where data might be multidimensional but not classified, i.e without variables(age, height, weight,...) Do I have the right understanding? If not, what exactly is classification? and why and when is it useful? Any help is much appreciated


title: How to back-calculate change from baseline from a p-value for a paired t test

text: I have a result from a study in a research paper which gives a mean/SD for 2 time points (pre & post treatment), and a p-value based on a paired-samples t-test. What I want to calculate is the SD of the change from pre to post. Is there any way I can back-calculate this using the p-value and mean/SD for each time point?   Thanks in advance- Julie


title: Measuring error in the lag when cross-correlating data

text: Cross-correlation can be used to measure a 'lag' (or shift or offset) between two data sets (e.g., data streams). Is there a standard way to measure the errors in the 'lag'? Or better yet, is there a _good_ way to measure those errors? (related: Directly compare subpixel shifts between two spectra -- and get believable errors!) Some more details: I'm performing 2D cross-correlations between images in order to achieve sub-pixel image registration. I would like to know how accurate the measured offsets are (i.e., ARE they sub-pixel?). These are images of astrophysical sources, but not stars - all of the emission is extended. For images with stars, you can measure the centroids of the stars to determine offsets between images: that is not possible for this data. The Gaussian approach suggested by Nestor is one that I've attempted, but I'm not entirely clear on some points, particularly how to account for individual pixel errors. I suspected that other fields might use this same technique but call it by another name, which is why I asked this question here.


title: Exploratory factor analysis using pooled longitudinal data

text: I am working on an analysis using exploratory factor analysis (EFA) with the common factor model. This question concerns a methodological issue. Any insights from someone who knows about EFA would be appreciated. We measured the wellbeing of people using a **15 question** survey. The same questionnaire was given the same **100 subjects** at **3 different times** (there were large gaps in time between when the surveys were administered, therefore we are assuming there was no time dependence [such as carry-over effect]). The goal is to use factor analysis to determine a fixed set of common factors and then calculate the factor scores and to compare the scores of subjects _between time points_. Below is my strategy for the analysis. Does this sound reasonable? Analysis Strategy: Combine all the data as it if were from a single location. That is to say, create a data set with 100 observations from time 1 and append the 100 observations from time point 2. Then append the 100 observations from time point 3. This would give a combined data set with a total of 300 observations. Using the above data, fit an EFA model using all 15 variables and 300 observations and obtain factors. After computing the factor scores, again treat the data as if it were from 3 different time points (i.e. 100 subjects measured at three time points). Analyze the factors scores with i.e. ANOVA to test if factor scores of a factor change between the time points.


title: Least mean square algorithm – conceptual doubt

text: I have used the least mean square (LMS) algorithm to estimate a signal in the presence of high chaotic and random noise. MSE values in db at each SNR for the coefficients is positive even though I converted to Db scale. What is the reason that LMS may not perform? For an MA(4) model, estimating the coefficients where regressors are composed of chaotic noise: $$u_t = e_t + a_1e_{t-1} + ...a_2e_{t-4} + AWGN $$ $e_t$ is obtained from a logistic map, so it is chaotic in nature. These are the values of MSE between the estimated and true coefficients for different SNR               SNR  MSE     0    3.32008926893368e+95     5    7.74457951051661e+50     20   8620.99943544499     etc      What is the reason for such high MSE values? The following is the code. The MA model is $u(t) = s(t) + 0.4s(t-1) + 0.2s(t-4)$ obtained from http://www.mathworks.com/help/econ/regression-model- with-arima-errors-specifications.html Here the $s$ is obtained from the logistic map as mentioned above. Then to $u(t)$ I have added additive white Gaussian noise. I must be doing something incorrect in the code. Please help rectify. Thank you EDIT               %function [e,w]=lms(mu,M,u,d);     % Call:     % [e,w]=lms(mu,M,u,d);     %     % Input arguments:     % mu = step size, dim 1x1     % M = filter length, dim 1x1     % u = input signal, dim Nx1     % d = desired signal, dim Nx1     %     % Output arguments:     % e = estimation error, dim Nx1     % w = final filter coefficients, dim Mx1     %inital values: 0     clc     clear all          A= 0.4;     B = 0.0;     C = 0.0;     D =0.2;               M = 4;  %model order     %number of samples of the input signal     N=256;     %Make sure that u and d are column vectors            trials = 5;                              for itr = 1:trials                   x(1) = rand();     % Pass it through the MA model.          for jj =1 : N              x(jj+1) = 4*x(jj)*(1-x(jj));      end     s =x;     y(1) = 0.0;     y(2) = 0.0;     y(3) = 0.0;     y(4) = 0.0;               for ii =5 : N         y(ii) = s(i) + 0.4*s(ii-1) + 0.2*s(ii-4);     end          d=y(:);  %desired signal     %u =y;       u = awgn(y,0.1,'measured','linear');  %linear scale     %u = awgn(y,40,'measured');  %DB scale     u=u(:);     w=zeros(M,1);          for n=M:N     uvec=u(n:-1:n-M+1);     e_samples(n)=d(n)-w'*uvec;     % Start with big mu for speeding the convergence then slow down to reach the correct weights      if n < 20             mu=0.32;         else             mu=0.15;      end     w=w+mu*uvec*conj(e_samples(n));          end     e_samples = 0;          e_coeff(itr,:)= (((w(1)- A).^2) + ((w(2) - B).^2) + ((w(3)- C).^2) + ((w(4)- D).^2))/4  ;     W1(itr) = w(1);     W2(itr) = w(2);     W3(itr) = w(3);     W4(itr) = w(4);               end     mse_coeff = mean(e_coeff);               MeanSqErr = [mse_coeff  mean(W1) mean(W2) mean(W3) mean(W4)]


title: Test to compare an effect on a variable on a daily basis

text: I would like to know the effect of an external factor on the number of visits on a website during a period of 4 weeks (only search visits from Google). For this, I measured the number of visits every day through Google research. For example I have:                     Date Visits     01/01/2014    456     02/01/2014    423     03/01/2014    399     ...      I did this for 4 months, including the 4 weeks when the external factor ran. I would like to prove if there is any difference in the mean number of visits during these 4 weeks. My first idea is to run a t-test between two similar periods of 4 weeks, one will be the 4 weeks with the external factor, and the other will be a series of 4 weeks without. I have difficulty selecting the test: should it be a t-test paired comparing each day of the 4 weeks? I don't really know if I can consider my samples as independent as there is a weekly seasonality (on the weekends the visits are lower than on the week)


title: Is it necessary to plot histogram of dependent variable before running simple linear regression?

text: I was working on an assignment. The data set was really simple, only consisting one independent variable $y$ and dependent variable $x$. Someone suggested me plot a histogram of $y$ before running simple linear regression. He told me that by plotting frequency at y-axis and value of $y$ on x-axis, and then examining whether it looks like normal distribution would help, but I don't see why that is relevant. This might be a newbie question, but any help would be great! **Update:** From what I understand, simple linear regression assumes that the $e_i$ follows normal distribution. Since $y_i=\beta_0 +\beta_1x_i +e_i$ by our assumption, then $y_i$ also follows normal distribution. But I assume that is not what a histogram would tell us...


title: joint p.d.f. of stationary time series variables

text: if a stationary time series verifies that each variable depends only on the variable before it, and the joint p.d.f. of xi and xi-1 is f(xi-1,xi), which is the joint p.d.f. of xi,xi+1,xi+2, and of xi,xi+1,xi+2,xi+3?


title: Probability proofs using ordered samples

text: Given an ordered i.i.d sample $X_{(1)}, \dots, X_{(n)}$ from a continuous distribution $F(x)$. How can it be shown that: (1) $\text{Pr}(X_{(k)} \leq x) = \text{P}r(N(x) \geq k)$ where $N(x)$ is the number of sample values less than $x$; furthermore, $N(x) \sim \text{Bin}(n, F(x))$ and (2) $p(x) = \frac{n!}{(k-1)!(n-k)!} (F(x))^{k-1}(1-F(x))^{n-k}f(x)$ where $f(x)$ is the density of $F(x)$. Beginning with the right-hand side of the equality given in (1), given $N(x) \sim \text{Bin}(n, F(x))$, it would seem that $\text{Pr}(N(x) \geq k) = 1 - F_{\text{Bin}(n, F(x))}(k-1)$ where $F_{\text{Bin}(n, F(x))}(k)$ is the CDF for the binomial distribution evaluated at $k$ with parameters $n$ and $F(x)$. From here, because (1) is an equality the left-hand side and right-hand side should be identical. Given that the term on the left-hand side in (1) is the CDF of $X_{(k)}$, the PDF ($p(x)$) could be found by solving: $\frac{d}{dx}(1 - F_{\text{Bin}(n, F(x))}(x))$ Hence, $p(x)$ given in (2) should simply be $f_{\text{Bin}(n, F(x))}$(x). Obviously, my attempt above isn't congruent with the solution. I was hoping somebody could help with this problem.


title: Ensemble of models with different feature spaces

text: **BACKGROUND** I have data in which the dependent variable is binary with a highly-skewed distribution: <1% records are 1 (doers), >99% records are 0 (non-doers). I'm using logistic regression to predict the **probability** that new records are doers. To handle this rare-event situation, I made multiple samples of non-doers that are size-matched to the number of doers (e.g., sample 1 has the 100 doers and 100 non-doers, sample 2 has the 100 doers and a different set of 100 non- doers, etc.). **QUESTION** **If I fit a logistic regression to each sample, how do I make an ensemble model to assign _probabilities_ to new records?** Samples have different observations and perform their own feature selection, so they have different feature spaces, which precludes averaging feature coefficients. **Do you have any suggestions for how I can build an ensemble model to take into account the models from all of my samples to compute _probabilities_?**


title: Skewed, humpy, left tail doesn't approach 0

text: Self-learning statistics by fiddling around with some modeling. Having trouble finding distributions for this data however.   * Defined for $x >0$   * $P(0) > 0$   * $P(X) \to 0$ as $X \to \infty$, (very quicky)   * One hump I haven't found any distributions that fit the description so far. Anyone got any ideas? Given the assumption of basically asymptotic sample size, how unadvisable would it be to just curve fit some polynomial of degree 10 to the data and call that the Taylor series of my weird unknown PDF?


title: Interaction in stepwise regression analysis

text: I did a stepwise regrssion analysis to predict energy expenditure using the variables, height, weight, age, gender and energy intake. The final model contains the variables gender and weight. Now does this final model take into account gender by weight interaction? Or do I have to construct a new equation from this final model, that will accomodate interaction?


title: Entropy estimation for a symbol sequence

text: I am looking for an R-implementation of the Lempel-Ziv data compression algorithm, to estimate the _source_ entropy of a time-series consisting of a sequence of symbols. Rather than simply measuring the entropy of the (time-aggregated) symbol- frequency distribution, the Lempel-Ziv algorithm can be used to estimate the entropy of the temporal ordering within the sequence itself. The sequences I am studying are derived from animal trajectories; each $x,y$ location is assigned an arbitrary label (A,B,...,Z). Represented in this way, the trajectory of a moving animal looks like B,G,G,S,Y,G,Y,H,D,S,A,B,B,G,G,G,H etc. Hence, I am interested in measuring the $\textit{predictability}$ of the sequence. The procedure for estimating the source entropy is (briefly) described in the Supplementary Information of Song et al (2010), and they describe the Lempel- Ziv estimator as: $$ S^{est} = \left( \frac {1}{n} \sum_{i=1} \Lambda_i \right)^{-1} ln~n $$ where $\Lambda_i$ is the length of the shortest substring starting at position $i$ which $\textit{doesn't}$ previously appear from position 1 to $i$-1. So, does anyone have any suggestions for how to compute the above in R? Incidentally, I should mention that I have multiple trajectories between which I would like to make entropy comparisons, however, some of the trajectories contain gaps that correspond to times when the animal was not observed. REFERENCE: Song, Qu, Blumm & Barabasi (2010) Limits of Predictability in Human Mobility. $\textit{Science}$. 19. http://www.sciencemag.org/content/327/5968/1018.short _**UPDATE:_** The function below is what I arrived at. However, rather than chugging through all the possible sub-sequences after the i'th index, and then finding the shortest novel sub-sequence, it would be faster to use a while loop. Similarly I am updating the dictionary of newly-acquired unique substrings using rbind, which is not efficient. OTOH the sequences I'm using are short (<300), so this isn't a problem, but there is certainly room for improvement.               PREDICTABILITY <- function(Sequence, Max_String)        {            output     <- matrix(NA,nrow=length(Sequence), ncol=1)       dictionary <- NULL       ## cycle through each entry, i, in the Sequence       for (i in 1:(length(Sequence)-Max_String) )         {         ## Compile list of increasingly-long substrings, starting at position i; i+0,i+1,i+2,...,i+Max_String         codons <- matrix(NA, nrow=Max_String, ncol=1)         for (STRL in 0:Max_String)           {           codons[STRL,] <- paste(Sequence [i:(i+STRL-1)], collapse="")           }         ## Find which of these codons have NOT been seen before         new <- codons[!codons %in% dictionary]         ## check for no new codons         ifelse ((length(new)>0),           record <- min(nchar(new)),    ## if we have new codons, find the shortest among them           record <- NA )                ## if none are new (because we aren't searching far enough ahead), assign NA...          ## find the shortest of these unseen codons         output[i,] <- record              ## Finally, add the unseen codons to the dictionary         dictionary <- c(dictionary, new)         }##i       ## Calculate source entropy (i.e. predictability) from formula in Song et al (2010)       n <- length(output[!is.na(output)])         S <- (1/mean(output, na.rm=TRUE)) * log(n)  ## Source entropy ?natural log means the units are nats?        return(S)}


title: Sybil detection metric

text: I have a set of user data and I want to build some kind of metric to evaluate the probability of the user being a sybil (a "fake" account). **But** I have a very limited set of users who are sybils with 100% certainty. How do I use machine learning here? Also, as for now, I've built a heuristic metric based on that data and need to evaluate it somehow. **To sum up:** I have a small fraction of data that is labeled and only negative class. And need to build a metric to evaluate users. On top of that I need evaluate the "goodness" of that metric? How do I approach this problem? **ps** It would be good if I could scale this process for big datasets.


title: What does it mean to explain variance?

text: In particular, I am wondering why we have this concept Multiple R (which I can understand as the correlation between observed and predicted scores in multiple regression), and then a separate concept R-squared which is just the square or R. I've been informed that R-squared is the percentage variation explained and R is not, but I don't understand the distinction that is being made between correlation and explained variation.


title: Proportional data with multiple responses

text: I am trying to analyze a data set in which the response variable is a proportion that can range from 0-1 (inclusive). This response was measured from the same individuals under four different experimental conditions. I am interested in determining the effect of these conditions on the response. I have read up on some possible means of analyzing these data via either tobit or beta regression (I realize that the latter can't be used when 0's and 1's are in the dataset, but I have found a transformation in the literature that fixes this problem). However, I have not found a way to implement such a model in R to account for the repeated measurements. I've gone through several packages that can perform such procedures but none of them seem to have this functionality (or I am completely missing it). I'm quite stuck here and would appreciate any suggestions as to how this can be done. As a background, the response variable is the proportion of times a female grasshopper responded to an acoustic signal (it is important to note that these are, and must be, averaged across several stimuli so that it is impossible to quantify "successes" and "failures" per se; thus they are true proportions and I can't directly use binomial functions). So these values can range anywhere from zero to one. This experiment was performed at 4 different levels of background noise, and I am interested in whether the noise affected the likelihood of female response. I ran these trials for 21 different females, all of whom were exposed to all signals in all noise levels.


title: What are the best criteria to select the model for Lasso regression?

text: I have two different formulations of the Lasso regression for the same problem. For each formulation, I selected the best model based on cross validation error. But Now, I want to compare two models that come from two different formulations of Lasso regression. The cross validation error of theses two models are close to each other, but the number of parameters is different. Now, the question for me is that, which test statistic should I use to say that, one model is significantly better than other model with taking into the number of parameters of the models.


title: Power Calculation - Calculation of Absolute Difference

text: I've been asked to undertake a power calculation because the regular statistician is not available. I am calculating power for a trial where the endpoint will be a percentage of patients not achieving remission at a given timepoint after receiving a drug or placebo. The ratio of case:controls is 1:1 I have two columns of a table to fill out (absolute difference and sample size), and have calculated required sample size at given response rates:               Power    %control   %cases   Abs Diff     Sample Size      0.8        35         5                       24      0.8        45        15                       34      0.8        60        30                       42      0.8        35        15                       71      0.8        45        25                       88      0.8        60        40                       97      I calculated the sample size using the function for comparing two proportions, 'pwr.2p.test', in the R package pwr               P0 <- c(0.35,0.45,0.60,0.35,0.45,0.60)     P1 <- c(0.05,0.15,0.30,0.15,0.25,0.40)     lapply(1:length(P0), function(i) pwr.2p.test(h=ES.h(P0[i], P1[i]), n=NULL, power=0.8))      I am struggling to identify what the column 'absolute difference' means in this context. Could anyone help?


title: Cov(a,bc) when Cov(a,b)=0, cov(a,c) not zero

text: Is there an expression for $\text{cov}(a,bc)$ in terms of $a$, $b$ and $c$ separately. More concretely, if $\text{cov}(a,b)=0$ and $\text{cov}(a,c)>0$, is $\text{cov}(a,bc)=0$?


title: When teaching statistics, use "normal" or "Gaussian"?

text: I use mostly "Gaussian distribution" in my book, but someone just suggested I switch to "normal distribution". Any consensus on which term to use for beginners? Of course the two terms are synonyms, so this is not a question about substance, but purely a matter of which term is more commonly used. And of course I use both terms. But which should be used mostly?


title: Intercept from standardized coefficients in logistic regression

text: I have fit a logistic regression model with original y and standardized x variables. Slope coefficients can be easily converted back to their original scale by $\beta^*_j/\sigma_{x_j}$ where $\beta^*_j$ is the output from the fitted model. However, I am not sure how to get the original intercept estimate from the output. Any help is appreciated. Thanks


title: Variance inflation factor for generalized additive models

text: In the usual VIF calculation for a linear regression, each independent/explanatory variable $X_j$ is treated as the dependent variable in an ordinary least squares regression. i.e. $$ X_j = \beta_0 + \sum_{i=1, i \neq j}^n \beta_i X_i $$ The $R^2$ values are stored for each of the $n$ regressions and VIF is determined by $$ VIF_j = \frac{1}{1-R^2_j} $$ for a particular explanatory variable. Suppose my generalize additive model takes the form, $$ Y=\beta_0+ \sum_{i=1}^n \beta_iX_i + \sum_{j=1}^m s_j(X_i) . $$ Is there an equivalent VIF calculation for this type of model? Is there a way I can control for the smooth terms $s_j$ to test for multicollinearity? Thank you. http://en.wikipedia.org/wiki/Generalized_additive_model http://en.wikipedia.org/wiki/Variance_inflation_factor


title: Scale of variables and the consequences on the solution

text: Let us say I have data set of distinct $x_i$. A Gaussian is fitted to it with maximum likelihood, obtaining some $\mu$ and some $\sigma^2$. I will also obtain a likelihood $\mathcal{L}$. Now, I copy the data set where I set $\tilde{x}_i = 2 x_i$. I estimate another Gaussian, obtaining $\tilde{\mu}, \tilde{\sigma}^2, \tilde{L}$. Since I am using ML I will get $2 \mu = \tilde{\mu}$, $4 \sigma^2 = \tilde{\sigma^2}$. Contrary to my intuition, $\mathcal{L} \neq \tilde{\mathcal{L}}$. This is maybe not so bad for a univariate Gaussian. But what if I perform model selection for factor analysis, Gaussian mixture models or more complicated models? Will the order of models ranked by their respective likelihoods be the same?


title: Why compare sum of squares with ANOVA (and not mean)?

text: I got the following result of a ANOVA analysis of five independent variables with 5,8,2,2,6 levels and 5 replications for each combination:                                  Df Sum Sq Mean Sq F value Pr(>F)         IV1                 4 129805   32451  243.35 <2e-16 ***     IV2                 7  67227    9604   72.02 <2e-16 ***     IV3                 1  64253   64253  481.83 <2e-16 ***     IV4                 1 396445  396445 2972.94 <2e-16 ***     IV5                 5  91672   18334  137.49 <2e-16 ***     Residuals        4781 637553     133                        ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1       I now want to quantifiy which parameter has how much effect on the variation of the results. In a statistics book and suggest by an answer to my previous question I found the suggestion to simply calculate the ratio of the Sum of Square of the total Sum of Square. In my case as the Total sum of squares = 1386955 this would lead to 9% influence of the IV1 and 45% influence due to errors. **Question: Why is the Sum of Squares used here and not the mean of the sum of squares?** The mean is calculated by dividing by the number of degrees of freedom. This seems much more logical to me: With more degrees of freedom, the sum of squares automatically gets higher so this must be compenstated somehow. The results would be totaly different: IV1 would have 6% influence and the error would go down to below 1%. Side question: Why is the Df for the residuals such a strange number? I thought the df of the residuals is a*b*c*d*...*(n-1), in my case 5*8*2*2*6*4=3840


title: Automatically detect "bimodal" curves in a large number of short time series

text: So, I have a huge number of amplitude modulated waves and I am attempting to determine which modulation parameter modulates each time series. To do this, I folded each time series over at the repetition frequency of the modulator wave, to produce an averaged modulation index for each time series. ![ "Unimodal-like" left, "Bimodal-like" right](http://i.stack.imgur.com/eNlC9.jpg) Here is a (not particularly good, but convenient to produce) example. The left curve, if smoothed, produces a single preferred time point. The right produces two (a weak preference at 7, and a strong preference at 17). Technically the same is somewhat true of the one on the left but, for the sake of argument, lets ignore that for now. I'd like to do something analogous to kernel density estimation for each time series, and then automatically determine the peaks from this shape. I know these aren't really distributions, but what I would like to do is automatically determine: A) the number of peaks per time series, B) the locations of these peaks. B sounds easy to me if I can do A. I suppose this means smoothing the curves and defining the local maxima via some sort of window, but thought some of you statistics experts might know of a proper method for this sort of thing. Thanks in advance, Joseph


title: How should we interpret the variable created by Principal Component Analysis?

text: I tried to model $\text{Saving} = a + b_1*\text{Income} + b_2*\text{Wealth}$ but found that $\text{Income}$ and $\text{Wealth}$ were highly correlated. I applied PCA to get a new variable $\text{New}$ based on the first component. Next, I modeled $\text{Saving} = c + b_3*\text{New}$. I then had two questions as follows:   1. How should I interpret the model: $\text{Saving} = c + b_3*\text{New}$?   2. If it would not make any sense to interpret $\text{Saving} = c + b_3*\text{New}$, should I replace $\text{New}$ with the corresponding weights (loadings) generated by PCA? The weights regarding $\text{New}$ were $1.5*\text{Income} + 2*\text{Wealth}$. (I mean should I do $\text{Saving} = c + b_3*(1.5*\text{Income} + 2*\text{Wealth})$ instead?)


title: Using logistic regression to fit a weird dataset

text: I have the variables $Y$, $x_{m1}$, $x_{m2}$, $\ldots$, $x_{mn}$, $x_{p1}$, $x_{p2}$, $\ldots$, $x_{pn}$, $x_{q1}$, $x_{q2}$, $\ldots$, $x_{qn}$. $Y$ is the class label (0 or 1). $x_{mj}$, $x_{pj}$ and $x_{qj}$ can be regarded as the same kind of property, but different combination of three properties can predict different $Y$. Actually, $x_{m1}$, $x_{m2}$, $\ldots$, $x_{mn}$, $x_{p1}$, $x_{p2}$, $\ldots$, $x_{pn}$ and $x_{q1}$, $x_{q2}$, $\ldots$, $x_{qn}$ came from three subclasses. The subclasses' label determine the final $Y$. I have such kind of a dataset, I want to do inference with logistic regression and then use the new model to predict new $y$ given $x$. The simplest method is to deem every $x$ as independent and do the logistic regression, I don' think this method is effective, using interaction may be another choice. Can anyone tell me how to solve this problem or have some other suggestions by using other models.


title: Is 0 a valid value in a Likert scale?

text: I have carried out my pilot study on language learning motivation using a 6 point Likert scale but from 0 (strongly disagree) to 5 (very much agree). I noticed a colleague in his survey used 1 to 6. Will my computed variables (sum and mean) be the same as if I had used 1 to 6? Is it normally recommended to not use a 0 for some reason? I am new to SPSS but have managed to do most of what I need to do but I'm now worrying my values are 'distorted'. I don't understand how SPSS adds a 0 into an equation. Thanks


title: Effect of mediating variable independent of effect of possible covariate

text: Here is a rather involved question. It can be summed up as follows: **how can I test the effect of a potential mediating variable, independent of the effect of a possible covariate with which it is confounded?** Here is a simplified description of my study and then the analysis problem I'm facing in more detail: > I recently conducted a study investigating the effects of different > instruction methods on participants' learning outcomes. Participants were > exposed to one of three TRAINING conditions. I then measured two types of > learning outcomes. The first measurement, SOLUTION, is a hand-coded binary > factor describing how well participants could verbalize the solution method > they were supposed to have learned during training. SOLUTION is mainly of > interest as a predictor of the other measurement, TRANSFER. TRANSFER is a > metric variable defined as performance on POSTTEST minus performance on > PRETEST. I was originally testing the following hypotheses.   1. TRAINING condition A would promote TRANSFER better than condition B.    2. TRAINING condition A would also promote better SOLUTION quality than B.    3. SOLUTION would mediate the effects of TRAINING on TRANSFER. For   1. I did an ANOVA with TRAINING as a between-subjects factor and TRANSFER as the DV. Unfortunately the effect was not significant, so (1) was not confirmed.    2. was confirmed, via a simple Chi-square analysis.  So my main question relates to how to test (3). Because (1) wasn't confirmed, I restate (3) as follows: SOLUTION has an effect on TRANSFER. In particular, better SOLUTION leads to higher TRANSFER. Now I already know that the data trend indeed is in that direction. To test whether this effect is statistically reliable, I would LIKE to simply add SOLUTION as a factor to the ANOVA I used to test (1). However, there is a complication. SOLUTION is correlated with PRETEST score, which in turn is correlated with TRANSFER. So, apparent effects of SOLUTION might actually be effects of PRETEST, right? Thus, do I need to add PRETEST into the model as a covariate? If so, I'd then be doing an ANCOVA, and I've been advised that doing so violates the assumptions of ANCOVA, namely that the covariate should NOT be correlated with the other factors, while in this case it IS correlated with one of the factors, i.e. SOLUTION. Any way around this? Alternate analysis I should be doing? Or is it not really a problem to begin with? If it matters, PRETEST is _positively_ correlated with SOLUTION, but _negatively_ correlated with TRANSFER. So maybe I'd be justified in leaving PRETEST out of the analysis because, although it is confounded with SOLUTION, their effects on transfer point in opposite directions, so leaving it out would only weaken my effect. Is this reasoning valid? I'm also wondering more generally whether it's OK for me to add SOLUTION as a factor to my original ANOVA at all, given that SOLUTION was not randomly assigned and is itself correlated with the other factor, TRAINING, as mentioned above. Is this a problem and if so, what should I do about it?


title: Conditionalising increases expected information?

text: I saw a talk that mentioned this theorem. I am now trying to find details of the theorem, but I'm not having any luck. I hope that if I can describe roughly what the theorem said, someone can recognise it and point me to some literature. So the theorem basically said something like the expected self-information of conditionalising on the value of some random variable is always higher than the prior's self-information. Does that ring any bells?


title: Pocket algorithm for training perceptrons

text: When you read about perceptron variants at Wikipedia there is explained an algorithm: Pocket Algorithm It is said that: > solves the stability problem of perceptron learning by keeping the best > solution seen so far "in its pocket" However, there is not much explanation about the algorithm, and I would like to see some pseudocode for it, as well as an explanation of how to implement by hand.


title: Regression with correlation structure

text: I have a theoretical question about regression models. Let's say I measured multiple responses from $n$ subjects and these responses are correlated with each other. For example, let's say I measured heart rate and body temperature from $n$ individuals with the following categorical factors: sex (male/female) and age (young/adult). It's quite possible that heart rate and body temperature are correlated to some extent in any individual. Ignoring the correlation between heart rate and body temperature I'll model each response independently as: $y_{heart\\_rate} = X\beta_{heart\\_rate} + \epsilon_{heart\\_rate}$ $y_{body\\_temperature} = X\beta_{body\\_temperature} + \epsilon_{body\\_temperature}$ where $X$ is identical for the two models and is of dimensions $n \times 2$, and: $\epsilon_{heart\\_rate} \sim N(0,\sigma^2_{heart\\_rate})$ $\epsilon_{body\\_temperature} \sim N(0,\sigma^2_{body\\_temperature})$ If I do want to take the correlation into account I would define $X^*$ as a block diagonal matrix of $X$'s. I.e., $X^* = \left( \begin{array}{cc} X & 0 \\\ 0 & X \\\ \end{array} \right)$ and $\epsilon^* \sim N_2(0,\left( \begin{array}{cc} \sigma^2_{heart\\_rate} & \rho \\\ \rho & \sigma^2_{body\\_temperature} \\\ \end{array} \right))$ For simplicity let's assume that $\rho$ (the correlation between heart rate and body temperature) is given. And my model will be: $ \left( \begin{array}{c} y_{heart\\_rate}\\\ y_{body\\_temperature}\\\ \end{array} \right) = X^* \left( \begin{array}{c} \beta_{heart\\_rate}\\\ \beta_{body\\_temperature}\\\ \end{array} \right) + \epsilon^*$ My question is whether $\hat{\beta}_{heart\\_rate}$ and $\hat{\beta}_{body\\_temperature}$ and their standard errors will be different under the two different models.


title: How to choose the split in Random forest for categorical predictors (features)?

text: I understand how best split is chosen for random forest for numerical predictors (features). Numerical predictors are sorted then for every value Gini impurity or entropy is calculated and a threshold is chosen which gives the best split. But how best split is chosen for categorical predictor as there is no specific ordering?


title: Data set to probability distribution to maximum likelihood estimation of sigma 1 and sigma 2

text: I have following 3 two dimensional datasets. Case 1: (Two continuous random variables) A = 1.3, 2.7, 3.9, 4.7, 5.6, 6.3, 7.5, 8.9, 9.1, 10 B = 7.4, 15.3, 24.4, 25.4, 29.6, 32.1, 34.5, 35.7, 27.8, 39.1 Case 2: (One discrete and one continuous random variable) C = 1, 1, 1, 1, 1, 0, 0, 1, 0, 1 D = 80.7, 85.3, 83.4, 84.6, 87.3, 87.8, 99.5, 83.5, 79.7, 82.4 Case 3: (One continuous and one discrete random variable) E = 1.3, 2.7, 3.9, 4.7, 5.6, 6.3, 7.5, 8.9, 9.1, 10 F = 1, 1, 1, 1, 1, 0, 0, 1, 0, 1 First, I want to represent them as a probability distribution and then want to estimate maximum likelihood of sigma 1 and sigma 2 parameters. I appreciate your help..


title: Permutation-based p-value to evaluate a classifier

text: I read an article entitled "Permutation Tests for Studying Classifier Performance" which explains how to use "permutation-based p-value" to test the performance of a classifier. I can not understand how to calculate the p-value. The article reports this formula:                p = (|{D′ ∈ D: e(f,D′) ≤ e(f,D)}|+1)/(k+1)      (DEFINITION 1)      Where e(f,D) is some error of the classifier f on the dataset D or D' (D' is the dataset obtained with permutation). This is the article: http://jmlr.org/papers/volume11/ojala10a/ojala10a.pdf But I do not understand the formula. Why the error with D' is less than the error with D? and can I calculate the absolute value of an equation with "<=" ? And what is k? The article says: > The empirical p-value of Definition 1 represents the fraction of randomized > samples where the classifier behaved better in the random data than in the > original data. Intuitively, it measures how likely the observed accuracy > would be obtained by chance, only because the classifier identified in the > training phase a pattern that happened to be random. Therefore, if the > p-value is small enough usually under a certain threshold, for example, α = > 0.05 we can say that the value of the error in the original data is indeed > significantly small and in consequence, that the classifier is significant > under the given null hypothesis, that is, the null hypothesis is rejected.