title: How to choose between z and t in 2 sample (non-paired) tests?

text: I understand that in the one sample case, z is employed where population variance is known and t for where it is unknown. In the **two sample case** : we are trying to test $H_0: \mu_x=\mu_y$ against $H_1: \mu_x\not= \mu_y$, for finite variance; $\sigma^2_x$ and $\sigma^2_y$; while $n_x$ and $n_y$ are respective sample sizes What test would we use if:   1. (both variance are equal) AND either      * (i) both population variance unknown (i.e. we use sample variance $S_x=S_y$)      * (ii) both population variance known ($\sigma^2_x=\sigma^2_y$ )   2. (variance not equal) AND either     * (i) both population variance unknown (i.e. we use sample variance, $S^2_x\not=S^2_y$)      * (ii) both population variance known ($\sigma^2_x \not=\sigma^2_y$ )     * (iii) one population variance known other unknown (i.e. given $\sigma^2_x$ and $ S^2_y$) * * * My answers:   1\. (i) use 2 sample t-test with $n_x+n_y-2$ degrees of freedom and variance: $S_p*\frac{n_x+n_y}{n_xn_y}$ (where $S_p$ is the pooled variance )   1\. (ii) use 2 sample z-test with variance: $(\frac{\sigma^2_x}{n_x}+\frac{\sigma^2_y}{n_y})$   2\. (i) no idea :(   2\. (ii) same as 1. (ii)   2\. (iii) no idea :(


title: Time series deseasonalization

text: How do you know when deseasonalization is not necessary? That is, from what I understand, if you want to just look at the trend and irregular components of a time series, then you just need to remove the seasonal component. However, let's say we're working with an irregular time series that is not affected by seasonality at all (for the sake of this example) and we want to see a trend of specific section of the remaining components. Would it be beneficial to remove the seasonal component anyways? Or would it skewer the remaining components? Any help is appreciated, thanks.


title: Is it still considered time series if one uses additional signals

text: Apologies in advance, time series is not my strength. Say I want to predict f(T+1) using f(T-1, T-2, ..., T-N) -- for example using a multi-level preceptron. If I want to enhance this using some other signals say using weather data on times T-1, T2, ..., T-N Would this still be considered Time Series Analysis? Thanks.


title: Overdispersed Poisson regression

text: In Gelman and Hill 2007 (http://www.stat.columbia.edu/~gelman/arm/), they mention that adding an $\epsilon \sim N(0,\sigma^2)$ term to a Poisson regression can be used to account for overdispersion. On a number of different sites I've seen this mentioned, but all refer you back to Gelman and Hill for an explanation of this, and the text doesn't address conceptually why this makes sense (although I believe them). I understand that a two-parameter distribution (eg negative-binomial) may be used in place of a Poisson, and this makes sense to me: it has an extra parameter to capture variance. However, it is not at all obvious to me why simply adding an $\epsilon$ term should capture overdispersion, as the resulting Poisson will still have mean equal to its variance (as always). Can anyone clarify what's going on here and why this should account for overdispersion, as an alternative to using a negative binomial? And is there any simple way of interpreting how much overdispersion this $\epsilon$ should capture? Thanks!


title: IRT/Rasch modeling with very large N

text: I want to fit a 1-parameter IRT model on a questionaire with 15 questions and about six million people. Considering the large N, standard errors aren't essential. It looks like the IRT world is sort of dizzying, and I was wondering if there were any tips as for what the proper software approach would be.


title: confidence interval and sampling distribution for the mean

text: Ok, I'm having trouble with standard error, sampling distributions and confidence intervals. The question I am trying to answer is: How many holes do I have to drill in the ground to get good idea of the fracture density (# of fractures / hole). I simulated a computer model with theoretical fractures. If I sample the model with 1 hole, the number of intersections could be from zero to about 5... Something in that range. It varies a lot, so I make multiple samples. Let's say I make 100 samples of 1 hole per sample. Then I can easily build a 95% confidence interval for the mean. Now, I increase to 2 holes per sample and perform 100 samples. In this case, when I build my 95% confidence interval am I using n = 100, or n = 200?


title: How to use out of sample prediction for zero utilization records in two-part health care cost model?

text: I'm considering a two-part, member-level model for predicting Year 2 inpatient health care costs based on health care utilization, demographics, and diagnosis data from Year 1. A logistic regression to predict Year 2 inpatient utilization (yes/no) + a loglinear OLS to predict Year 2 log cost (with perhaps a Duan Smear factor adjustment) given Year 2 costs > $0 seems appropriate. I'll train my model on 2011-2012 data, then score 2013 records to make 2014 predictions. However, I've encountered a major limitation with using a two-part predictive model: **How can I make severity predictions with the OLS component of the two-part model when I don't know a priori which members will have non-zero health care costs in Year 2 (2014)?** I've read a few forum posts recommending **out of sample prediction** for the zero utilization records in a two-part model - how is this done?


title: Odds of births of four generations

text: What are the odds of the following: Great Grandmother was born `January 14 (1/14)` My Grandfather was born `April 11 (4/11)` I was born on `January 14 (1/14)` My son was born on `April 11 (4/11)` FYI: My Great GrandMother lived to be 100 years old. She passed just a few years ago.


title: Asymptotically unbiased estimator using MLE

text: I am learning Maximum likelihood estimators for a inference class. And this is a problem I came across. Let $X_1,X_2,X_3,\ldots, X_n$ be a random sample with p.m.f $$p(X)=\theta(1-\theta)^x; x=0,1,2,\ldots\quad \mathrm{and}\quad 0<\theta<1$$.   As maximum Likelihood estimator I obtained $\hat{\theta}=\frac{1}{\bar{x}+1}$. Now the question asks to obtain a unbiased estimator for $\theta$ when $n$ is large and considering it's distribution for large sample approximate 95% confidence interval for $\theta$. What I am unable to do is to find a unbiased estimator. Then if $\hat\theta $ is unbiased for $\theta$ using the property that $\hat{\theta}\sim \mathrm{N}\left(\theta,{1\over I_x(\theta)}\right)$. I can construct a confidence interval. But how can I find a unbiased estimator for $\theta$?   I thought since **MLE's are not in general unbiased but are asymptotically unbiased to compute $\mathrm{E}\left(\frac{1}{\bar{x}+1}\right)$**. Here $\sum X_i$ follows a **negative binomial distribution** right? Then $\mathrm{E}\left(\frac{1}{\bar{x}+1}\right) =\mathrm{E}\left(\frac{n}{T+n}\right)$ where T=$\sum X_i$. Then I come up with $\sum n\left(\frac{1}{t+n}\right)\binom{n+t-1}{t} \theta^n(1-\theta)^t $. I am unable to show that this is asymptotically unbiased. Please help me to find a unbiased estimator for $\theta$ when $n$ is large.


title: Reliability Measure for Unequal Subjects between Raters?

text: I am comparing two medical imaging analysis techniques, and I would like to assess the inter- and intra-rater reliability for using each method. I am planning the study design, and would like help with finding an appropriate statistical analysis method. A description of the analysis techniques follows: For both methods 1 and 2, analysis requires filtering of the data to remove artifacts caused by the data collection process. Method 1 requires the rater to manually remove erroneous waveforms before analyzing each wave (i.e. subject). Method 2 requires the rater to choose filter settings before it automatically removes errors and analyzes each wave. The rater is blinded to the data set when performing the analysis with each method so the number of subjects varies between raters and methods. For a set of 10 waveforms, for example, the rater may remove waves 1, 3, and 5 using method 1 but remove waves 1, 3, 6, and 8 using method 2. I expect a different number of subjects for each method within a given rater and between raters. It is possible to note the waveform removed for each method, but informing the rater of this information may bias the analysis. A total of 4 raters will analyze a dataset containing 4000 waveforms (i.e. subjects) with each method on 2 separate occasions. Does anyone have a suggestion for a statistical analysis method to assess the reliability between these two methods? I would calculate ICC values, but I do not think unequal subjects between raters are possible.


title: Which method to find which features explain the target?

text: I have many categorical and numerical features and I sub-selected a few of them by a Random Forest. Now I'd like to determine which features probabilistically(/causally?) explain the target. Something that I can tell to business people which isn't a black box like a Random Forest. I can plot the class density versus some variable A, but sometimes it feels like there could be an alternative indirect correlation relation using features B, C, ... which better explains why the class density shows a pattern versus variable A. Maybe like when features explain away others. How can I determine which of these (marginal) class density is most useful to show to business people and ideally even be causal? I heard of graphical models - do they help? (Features can have a positive or negative correlation or sometimes two features together are needed since they cannot be correlated individually to the target)


title: Using copulas to sample a probability distribution

text: Say I have two random variables, X and Y. Their joint probability density function is a uniform distribution inside the triangle with vertices at (0,0), (0,1) and (1,2). The area is 1 so the joint pdf is 1. Calculating the marginal cumulative distribution function I get P_X = x^2 and P_Y = y - (y/2)^2. The question is how do I hook up two random number generators U(0,1) and V(0,1) to generate a sample of points distributed according to the above by using copula theory? Note that we have the full information and do not need to assume anything about the relationship of the marginals. I have inverted the cdfs but the problem is that the marginals do not encode the triangular restriction so some of the points generated will fall outside the specified triangle. Thanks very much for any help.


title: I've found two equations for regression slope, but they give me two different answers. What am I missing?

text: According to my notes, in a statistical model where $$Y_i=\beta_0 + \beta_1(x) + u$$ (where $u$ is the error term) the predicted slope is $$\hat{\beta}_1 = \frac{\operatorname{Cov}(X, y)}{\operatorname{Var}(X)}$$ According to this Khan academy video the slope equals $$\frac{\overline{X} \overline{Y} - \overline{XY}}{\overline{X}^2 - \overline{X^2}}$$ When I calculate the two given a basic data set the numbers are WAY off. Are these equations related, or am I way off here?


title: What are possible theories/models to describe social influence based on frequency?

text: I am working on social influence based on frequency, meaning if a person repeatedly appears at a location many times then he/she is more influential at that location. In other words, he/she can influence other people at that location with a higher chance of success compared to someone who is not there often. For example, if a professor at a university has worked there over many years, therefore he has been seen at the school many times, thus he should be influential at that school. I have been looking for a mathematical/statistic model that can, or at least, can partially describe the influence quantitatively based on the frequency of visiting a particular location. So my questions are: 1) Do you have any suggestions of what in math/stats that may be related to this? I am just asking for your opinions/suggestions in general, which do not have to be correct, but can still be immensely helpful. Note that I assume friendships have already be considered in this social influence problem. Here, I just pay attention to frequency of visits. The frequency may be related to social status, social power, which in turn can affect social influence. 2) Regardless of question 1, do you think the assumption "the higher the frequency, the more influential" is intuitively correct? Any response/agreement/disagreement/opinions would be much appreciated. Thanks


title: Methods to combine ( e1071 svm ) models in R to generate a more complete, accurate model

text: I am using the function `svm` from the package `e1071` in `R` to generate a support vector machine model. I have a **very large** data set, and for the moment, while in an exploratory mode, want to simply read in small slices of the data that can be modeled on my single machine.      After obtaining the model results, I would like to read in more data and generate another SVM model, etc. until all of the data is modeled, which will generate about 50 different models. I would then like to merge all of these models / results together in some fashion in order to get a grand model which would approximate what I would achieve if I could fit all the data in at once.      I know there are multiple ways to do this, theoretically. But, sticking with SVM modeling, what are my options in `R`? (I want to stick with the `e1071` package because it has some things about it that are not in the others that I saw.)      If it matters, my data is **_not_** genomic, and it is highly weighted (meaning that I am modeling T/F, and usually find an F, but the feature set data is not sparse).      Below is a "pseudocode" snapshot of what I'm proposing / hoping, in case that helps clear things up:                    require(e1071)     modelSet = NULL     for (i in c(1:50) ) { ## imagine there are 50 files       dF <- read.csv(paste("file",i,"csv", sep=".") )       modelSet[[i]] <- svm(myOutput ~ ., data= dF, probability= TRUE)     }     ## Now I would like to find a way to merge all 50 of the "modelSet" models together     ## to make 1 composite model derived from all 50 data files.            Thanks!   Mike


title: How does including a random effect, change the parameter estimate for a group level covariate?

text: I'll try and explain where I am getting stuck. I wish to model delay to treatment. I observe patients nested in hospitals. I have a mixture of patient and hospital level covariates. I suspect that there are unmeasured covariates that induce correlation between patients being treated in the same hospital. This would violate the basic assumptions about independent and identically distributed variables, and to get round this I add a random or fixed effect to the model. Without this, then the unmeasured covariates will be operating through the measured hospital level covariates and induce bias. If estimate a fixed effects model then I would have an individual mean 'delay' for each hospital, but am I correct in thinking that these would be 'collinear' with any measured site level variable, and therefore not identifiable. Does this also apply to a random effects model? And if it does then is the aforementioned bias completely partioned to the 'random effect' or is some still operating through the measured site level covariates? I have had a look at this great answer, and read the Bafumi and Gelman reference that was suggested, but I would very much appreciate any guidance.


title: Standard method for calculating contribution of individual variable to outcome

text: I'm looking at the change in vote shares for an election between two periods. I'd like to say something about the contribution of one variable in particular towards the election swing I document. I'm wondering if there is a standard approach in the literature. I'd like to be able to say something like X accounts for 5 points of the 20 point swing. Is there a standard way to calculate a variable's contribution? Predict the outcome with X set equal to zero and compare the difference? Standardize the coefficients and take the value of the beta coefficient as the percent contribution?


title: testing for the difference in slopes if independent data is correlated

text: So I am looking at whether the effect of rainfall on crop yield for 8 different States is similar across these States i.e. I am testing whether the slopes between rainfall and yield is significantly different from each other among the eight States. For each State, I have yield and rain data for 48 years. I ran the following code in r one with interaction between rainfall and State and other without interaction and then compared the model.                   mod1 <- aov(yield ~ rainfall*state)         mod2 <- aov(yield ~ rainfall + state)         anova(mod1,mod2)      If model shows that removing the effect of interaction significantly affect the fit of the model, then I can conclude that slopes of yield vs rainfall are different for different States. So far this is ok and please let me know if I am wrong in interpreting this. My problem is this: my rainfall data (my independent variable) collected over time for each State will not be independent of each other. In this case, will be results will be wrong? What other set of analysis I need to do if my results are wrong? Thanks for the advice.


title: Finding significant peaks -- evaluation of different methods

text: I have a small doubt. My real data looks like this Y values are random values of integers from 0 to 2000. X values run like 1,2,3,4,5,.. to 2 million. Now, my task is to identify significant peaks and remove background noise. To achieve this , I used 2 methods : I created a sliding window of 50 elements. I divide y[50] by mean of y[1]..to y[49].Now, the modified values are plotted and the significant peaks are recorded. I created a sliding window of 50 elements. I divide y[50] by y[49] and y[50] by y[48] and so on till y[50] by y[1], and take the average.Now, the modified values are plotted and the significant peaks are recorded. Which is the best method to use? You can also suggest other methods. Background info: The 2 million numbers are the nucleotide positions of a an entire genome. For eg if the nucleotide code is AUGCAUC .. and so on 1st position correspond s to A , 2nd to U and so on. For each of this position,you have a coverage value which is important to predict other things. First I calculate these values, and see if the peaks correspond to the results. After visual inspection I have to empirically choose the cut off value of the peak. Say for eg 100. I am more interested in the method. The first method does not seem to work out well. If the mean of previous 49 elements is low then the final value shoots up. The second one seems to be fine. But, I need a good opinion from a mathematician.


title: What is a least angle regression?

text: Conceptually, I don't understand what is a least angle regression and why it solves LASSO http://www.cc.gatech.edu/~isbell/reading/papers/lasso_simple.html.pdf We know that LASSO is $$\min_x||Ax - y||_2^2 + \lambda|\beta|_1$$ From my understanding LARS is exactly like line search from CGD, where you take the variable that gives you the most desire result and iterate through. Can someone give a walk through on how LARS solve LASSO? Thanks in advance


title: building model decision trees Bayesian network

text: I am analyzing marketing data and Y response are as follows. There are 2 different binary Y columns. Each row of the data has data about a customer and what he bought               a customer bought product A - Yes or No     a customer bought product B - Yes or No      In order to build a decision tree, I am planning to merge these columns into a single column (NOMINAL values). How should I create that column? Is below method okay? I am a bit hesitant because number 3 below will appear a very few times. Please let me know your inputs.                  1) if a customer bought product A then lets say that column has value 1        2) if a customer bought product B then lets say that column has value 2        3) if a customer bought product A and B both then lets say that column has value 3        4) if a customer didnt buy either A or B then lets say that column has value 4


title: Leave-one-out cross-validation for negative binomial

text: I've just started working with the GLM negative binomial and I was wondering if it's possible to carry out a leave-one-out cross-validation for negative binomial GLM as I'm having troubles doing it with R using {boot}. My variables are standardized!Can it be a problem? Any help will be very very appreciated! Sorry if this question is a bit prosy..Thanks in advance! Mau


title: Is there a way to customize my likelihood function for logit models using speedglm/biglm/glm packages?

text: My goal is to fit a custom logistic regression/survival analysis function using the `optim`/`maxBFGS` functions in R and literally defining the functions by hand and then extracting the betas. I was always under the impression that for the packages `speedglm`, `biglm`, and `glm`, the likelihood functions for logit models or whatever distribution were hardlocked. However, I was wondering if I was mistaken or if it was possible to specify my own likelihood functions. The reason being that `optim`/`maxBFGS` is a LOT slower to run than `speedglm`.


title: How to compare 4 independent means when there is no dependent variable?

text: There is one other question similar to mine (Compare independent numeric vectors (means) in R), but has not been answered (1 year has passed), so I'm repeating the question here. I have obtained exam scores (0-100) of 200 students in 4 subjects. Two of the subjects are taught using one teaching method, while the other two are taught using another method. The subjects are unrelated to each other (a good score on one does not facilitate a good score in another). Scores in one of the subjects are almost exclusively in the 50s. The other three subjects have consistently higher scores on average (mid to high 60s). I wish to see if there is a difference in scores between the two teaching methods. Initially I thought I could use ANOVA, but then realized that there is no dependent variable here. What are my options? What is the best way to approach this problem? Thanks in advance.


title: Why not use the T-distribution to estimate the mean when the sample is large?

text: Basic statistics courses often suggest using a normal distribution to estimate the mean of a population parameter when the sample size _n_ is large (typically over 30 or 50). Student's T-distribution is used for smaller sample sizes to account for the uncertainty in the standard deviation of the sample. When sample size is large, the sample standard deviation gives good information on the population standard deviation, allowing for a normal- distribution estimate. I get that. But why use an estimate when you can get your confidence interval exactly? Regardless of sample size, what's the point of using the normal distribution if it's just an estimate of something you can get exactly with the T-distribution?


title: piloting a multinomial logistic regression model

text: I have 11 variables in my data set.   farmers Group(1,2,3,4 this is my dependent variable) Independent variables   Total holding ,Crop area , barn capacity.....and barn capacity extent match and YPH.   These last 2 variables are calculated using the existing explanatory variables in the data set(e.g. YPH=volume/crop area). So obviously those have large correlation between them. Anyway I plotted multinormial logistic regression for my data set but only 3 variables are significant                                      Coef    S.E.   Wald Z Pr(>|Z|)     y>=2                   -1.3878 1.6039 -0.87  0.3869       y>=3                   -2.6493 1.6222 -1.63  0.1024       y>=4                   -3.8260 1.6332 -2.34  0.0192       Total.Holding..ha.      0.0032 0.0940  0.03  0.9731       Crop.Area...ha.         0.4709 0.4091  1.15  0.2497       Ban.Capacity            0.3288 0.3808  0.86  0.3879       Capacity.Extent.Match. -0.0061 0.0035 -1.74  0.0822       Volume...kg.           -0.0002 0.0001 -1.12  0.2629       YPH                     0.0013 0.0005  2.66  0.0077       X1F.                    0.0710 0.0272  2.61  0.0092       R.Grades..              0.0552 0.0502  1.10  0.2714       FIDS..                 -0.0191 0.0171 -1.11  0.2649       Net.Return....         -0.0025 0.0013 -1.99  0.0461       This should not be acceptable by me. Please help me. 1) Is this procedure correct? 2) Can I use this model when I have correlated explanatory variables? 3) Here I had a variable called farmers mark and I categorize it and created a variable called farmers group, then I used for logit model.   Normally it was a continuous variable. So, is there any other model that can I use for this data set(rather than simple regression model)? code I used in R               library(rms)     View(polo9variables)      attach(polo9variables)                Y<-cbind(Group1.no.recovery.)          X<-cbind(Total.Holding..ha.,Crop.Area...ha.,Ban.Capacity,Capacity.Extent.Match,     Volume...kg.,YPH,X1F.,R.Grades..,FIDS..,Net.Return....)           table(Y)           XVar<-cbind("Total.Holding..ha.","Crop.Area...ha.","Ban.Capacity","Capacity.Extent.Match",     "Volume...kg.","YPH","X1F.","R.Grades..","FIDS..","Net.Return....")           ddist<- datadist(XVar)           options(datadist='ddist')          ologit<-lrm(Y~X,data=polo9variables)          print(ologit)


title: How to obtain estimates for all levels in a mixed effects model that uses effect (deviation) coding?

text: I am running a binomial mixed effects logistic regression in R using `glmer` for a sociolinguistics project. I was asked to used deviation (effect) coding. From what I gather, in deviation coding the last level in a factor is assigned -1, because this is the level that is never compared to the other levels within that variable. Is it possible to obtain the `Estimate` (`Exp(B)` value) for the last level as well by using function `relevel`? I need to report the estimates for all the levels. For example, my model has the independent variable called **Orthography** with four levels (`s`, `sh`, `s1`, `sh1`). The dependent variable is **produced sibilant**. In deviation coding the fourth level (`sh1`) will not be compared to the other three levels, and estimates will be available for the first three (`s`, `sh`, `s1`). The intercept is the mean of the means of all four levels (`s + sh + s1 + sh1 / 4`). I am interested in obtaining the estimate for the last level (`sh1`) as well. Does anyone know how to get that? Do I have to rerun the model by changing levels? If so, does anyone know how to do that? I have been unsuccessful with using function `relevel` to do this. I have other terms in my model as well:   * following segment, which has two levels (`vowel`, `consonant`),    * position of sibilant in word (`initial`, `medial`, `final`),    * grammatical function (`noun`, `verb`, `adjectives`), and    * language of instruction (`English`, `Gujarati`). This is the code for my model:               model.final_si = glmer(prod_sib ~ orthography + foll_segment + word_position +                                        grammatical_func + language_instruction +                                        (1|participant) + (1|item),                             family="binomial",data=data)


title: Linear modelling question and answer book

text: This book gives 500 question and answers on introductory statistics. Are there any books out there which give questions and answers specifically on linear models (simple linear regression, multiple regression, logistic regression etc)?


title: Need help to calculate and maximize cost-effective solution

text: Suppose we have a task that can be a success or a failure. We pay a fee and get a success for any attempt at the task. We can pay more in order to get a higher success rate. Given sufficient resources, we can guarantee success. For my example this shows the percentage success (i.e. probability of success on any attempt) followed by the price we pay in some monitary units.               Success rate (P) - Cost (C)     20%  -  10     40%  -  20      60%  -  40     80%  -  80     100% - 160      The example shows that if we pay 10 units we will have a 20% success rate. If we fail, the next time we will pay the same price. We are permitted to retry anytime until we are successful. I would like to find the most cost-effective solution in terms of expected cost. Is it possible to do that?


title: statistical comparison of power spectrums over all frequencies

text: So I have N=147 control time-series signals and N=134 treatment (THC) signals. I want to compare if there are any significant differences in frequency power in the two groups. The top part of the figure below shows the median frequency power for each frequency of the two groups. Thus for each frequency, I run a 2-sample t-test to get the mean difference and confidence bounds for each frequency. The mean difference (black line) and confidence bounds (grey shading) are shown in the bottom of the figure. Obviously, whenever the CB's dont cross 0, it is a significant difference. Here is my issue: **The two groups have different overall power** (where power=E(x(t)^2]). Thus if I ran the t-test on the naive data, it would show that the group with more power has higher power across all frequencies. To avoid this, I normalized each signal by its power, so that they all had a total power of 1. This revealed a significant difference from 5-10Hz, which is exactly what the theory predicts, but also showed a negative difference from 12-20 Hz. I suspect this is an artifact of the normalization method, since all the power must now be equal, and the control group had more power in the 5-10Hz range, it meant that the normalization lowered the power in the 12-20 Hz range to compensate, resulting in the observed significant difference. So in conclusion, my questions are   1. how can I avoid this artifact of normalization?   2. Is there any theory/literature on the method I am using? Thanks so much guys! ![img](http://i.stack.imgur.com/aDGsp.png)


title: Prediction interval for robust regression with MM-estimator

text: In their book "Robust Statistics", Maronna et al. consider the following model for robust regression: $y_i = \beta x_i + u_i$, where $u_i$ are independent of the $x_i$, and are i.i.d, with finite variance. They go on to provide a robust estimate $\hat{\beta}$ of $\beta$ which is asymptotically normal and give the covariance matrix for $\hat{\beta}$. My question is that without any knowledge about the distribution of $u_i$, is it possible to provide a prediction interval for $y$ (without using bootstrap)? I'm asking this because               library(MASS)     robustModel = rlm(formula = myFormula, data = myData, method = "MM")     predict.rlm(robustModel, newdata = myNewData, interval = "prediction")       in R generates a prediction interval. For reference, this is the code for predict.rlm:               predict.rlm <- function (object, newdata = NULL, scale = NULL, ...)     {     ## problems with using predict.lm are the scale and     ## the QR decomp which has been done on down-weighted values.     object$qr <- qr(sqrt(object$weights) * object$x)         predict.lm(object, newdata = newdata, scale = object$s, ...)     }      It seems to me that the prediction interval that is obtained this way is for normally distributed $u_i$. Is that correct? What am I missing here?


title: Is Fisher's LSD as bad as they say it is?

text: When we perform experiments (on small sample sizes (usually the sample size per treatment group is about 7~8)) on two groups, we use a t-test to test for difference. However, when we perform an ANOVA (obviously for more than two groups), we use something along the lines of Bonferroni (LSD/# of pairwise comparisons) or Tukey's as a post hoc, and as a student, I have been warned off from using Fisher's Least Significant Difference (LSD). Now the thing is, LSD is similar to pairwise t-test (am I right?), and so the only thing it doesn't account for is that we're doing multiple comparisons. How important is that when dealing with say 6 groups, if the ANOVA is itself significant? Or in other words, is there any scientific/statistical reason for using a Fisher's LSD?


title: Regarding switch from pure mathematics to machine learning: opinions and information

text: I'm a first year postdoc in pure mathematics (geometry/topology with strong background in analysis) with undergraduate statistics and probability background (with also measure theoretic probability). I also have some programming knowledge in ForTran, C and Matlab, but I never used them in my pure mathematics career in my graduate school. In my next job, I'm considering doing a postdoc in machine learning. The reason behind this switch is: I'm satisfied with my pure knowledge so far and have been wanting to see some real-life applications of mathematics, and also keep my options open in both industry and academia. My questions are: 1) How hard is the switch going to be? I guess I've all the required mathematics background, but will it be hard to pick up the necessary computer science skills, even if I work in more theory-oriented problems? 2) Is there a website/email-list where I can get notifications on jobs in machine learning? I'm looking for jobs in Europe mostly. Thanks and much appreciation :)


title: Funny statistics exam answers

text: Correcting exams is probably the most boring task of a teacher. But it might be amusing for us to collect funny statistics exam answers. One entry per answer.


title: Neural Net on results from Partition Bootstrap Forest - validity?

text: I am seeing some promising results with this method for my 'nutrition' dataset but am wary about modelling on results of a model. **Problenms with this data** I cannot seem to reduce my many correlated variables in a way that is acceptable. -PCA factorises the variables and indeed does successfully remove correlation but i loose the original meaning of the X predictors which is important. I can almost make sense of the PCA components but not enough to a satisfactory level, i.e factor one appears to be protein and factor two fats/oils. Also my matrix is not positive definitive and i still do not quite understand how to solve this, I suspect the cause is due to heavy correlation. -I am wary of removing variables manually as some people seem to think it is acceptable others not. Also some of the correlations are very difficult to justify removal and there are so many interlinks. **Cleaning the dataset issues** -Outliers are a huge potential problem in that there are quite a few 'chunks' of outliers, perhaps 5% of the data in some variables, beyond the usual definition of an outlier but this is quite possible a true result due to the nature of the data. So in my mind it is not justified to just remove them -The outliers that are problematic are in non normal X predictor variables which further hinders justification, it is possible they are errors but equally possible they are not **Bootstrap to Neural** Partition seems to be the best initial method for bypassing all these problems. Specifically bootstrap forest. It manages to pull out variables that make sense in real terms AND variables that are consistent when using different randomly generated 3 way validation columns - K fold always seems to fit too far and the validation agrees with this overfit. Is this possibly due to correlated colinear variables? The bootstrap method does fit too far in the training set but validation and test do agree. I cannot seem to get the overfit down here. So I took the variables that the bootstrap has produced and put them into a ANN. This ANN appears to have less over fit in the training and also a better fit in the test set. This is also true across multiple validation layouts. Is this model on model approach acceptable? I should add that neural network really doesn't like the dataset before variables are removed (by partition), it makes very poor choices which are absolutely not predictors, I really don't understand why it picks the variables it does in this situation. They are variables that i would typically remove


title: Standard errors in weighted least squares on aggregated data

text: I am interested - mostly just for my own knowledge, and not for any real problem - in the use of weighted least squares to estimate a model on individual-level data and aggregated versions of those same data. Here's some simple R code to show what I'm talking about. Basically, I estimate a regression model `lm1` on some data (`old`), then I aggregate the data (to create `new`) so that I only have one observation for each combination of $X1$, $X2$, and $Y$, and weight based on the number of observations at each combination of $X1$, $X2$, and $Y$. The resulting coefficients in a regression (`lm2`) on `new` are the same as those in `lm1`, but the standard errors are larger. But they are larger by a fixed amount (i.e., the SEs from `lm1` are proportionate to one another). How do I figure out what that ratio is? I've searched around for calculations of SEs for weighted least squares, but that has proved unhelpful.               set.seed(1)     n <- 1000     old <- data.frame(x1=sample(1:3,n,TRUE), x2=sample(1:3,n,TRUE))     old$y <- old$x1 + old$x2 + sample(-3:3,n,TRUE)         old$i <- with(old, interaction(x1,x2,y, drop=TRUE))     levels(old$i) <- seq_along(levels(old$i))     s <- split(old, old$i)         new <- do.call(rbind, lapply(s, `[`, 1, 1:4))         new$w <- sapply(s, nrow)          lm1 <- lm(y ~ x1 + x2, data=old)     lm2 <- lm(y ~ x1 + x2, data=new, weights=w)          # coefficients equal?     all.equal(coef(lm1),coef(lm2))     ## [1] TRUE          # ratio of SEs     summary(lm1)$coef[,2]/summary(lm2)$coef[,2]     ## (Intercept)          x1          x2      ##   0.2453172   0.2453172   0.2453172      Why are the SEs from the aggregated model proportionate to those from the individual-level model? And how do I know what their ratio is?


title: Is building a multiclass classifier better than several binary ones?

text: I need to classify URLs into categories. Say I have 15 categories that I'm planning to zero down every URL to. Is a 15-way classifier better? Where I have 15 labels and generate features for each data point. Or building 15 binary classifiers, say: Movie or Non-Movie, and use the numbers I get from these classifications to build a ranker, to pick the best category, going to be better?


title: Best way to go about correcting for multiple comparisons (or should one even) for a small dataset?

text: Disclaimer: I have a very minimal stats background and am just getting into this stuff. Corrections to any terminology I use or an improved rephrasing would be super appreciated. I have a simple science/perception experiment that I've run, testing a single group of subjects on 2-dimensions of conditions: speed (3 possible values) and "encoding type" (3 possible values). On a given trial, a subject is presented with a stimulus of a random "encoding type" presented at a random speed. I have a small number of subjects (~10). For each subject, I've calculated their performance (an observation) for each speed+encoding type combination. I'm interested in only seeing if there is a significant difference in performance of of "encoding type" restricted to a given speed (e.g. is encoding type 1 at speed 1 different from encoding type 2 at speed 1). I don't care of comparisons across speeds (e.g. I don't care if encoding type at speed 1 is different from encoding type 1 at speed 2). I'm doing 9 comparisons total (3 for each speed). So far, I've only run some paired ttests (would this even be the right thing to do?). I am getting some significant results without correction (p values ~0.03).   * Is a paired ttest the best test to use?   * Should I be correcting at all?   * Is there something other than bonferroni/sidak I should try? I've attempted to use bonferroni correction, where I took 0.05/9 = 0.0055. This throws out a number of results. Should I be dividing by 9 (the total number of tests I've performed) or by 3, the total number of tests constrained to the speed condition?


title: Which distributions on [0,1] other than the beta distribution form nice compounds with the binomial distribution?

text: For which distributions x, other than beta, is the x-binomial distribution nice? The beta and binomial distributions are famously conjugate but I am curious if other non-conjugate distributions will give comparably simple compound pmfs. By nice I mean that the pmf is nice to compute without resorting to numerical integration; I'm not referring to the easyness of sampling from the compound distribution.


title: Expected log value of noncentral exponential distribution

text: Suppose $X$ is non-central exponentially distributed with location $k$ and rate $\lambda$. Then, what is $E(\log(X))$. I know that for $k=0$, the answer is $-\log(\lambda) - \gamma$ where $\gamma$ is the Euler-Mascheroni constant. What about when $k > 0$?


title: Gensim Topic Modeling

text: I want to build a standard topic classifier. I was told gensim is the way to do it. I have difficulty training the gensim system. How do we provide a training data in a fast way. Some forum suggested Wikipedia training, and I used it. But the code WikiCorpus('....en.bz2') is running since last 6 days with no result. Is there a quick solution to the problem?


title: What test is appropriate for binary outcome with repeated measures and binary, ordinal and continuous independent variables?

text: I'm trying to figure out what is the most appropriate statistic test for use with my data and was hoping for some advice. The primary data consists of a binary independent variable (patient test - positive or negative; there are some missing values which are being excluded for now) for which we want to determine if it is predictive for a given clinical outcome. The outcome is ordinal (negative or grade 1 - 4), however we are interested in any positive outcome so this could be simplified to a binary outcome (negative or positive), although if the test is only predictive for some grades this would be good to know. In the simplified case, I could use a chi-square test. I'm less certain how to deal with the ordinal outcome. Additionally, each patient has multiple tests and outcomes measured at different times, but the number of repeats and the times are not the same for all patients. There are also a number of other data pertaining to each patient, some of which are binary (eg. sex, one of 2 research centers, certain other test results), ordinal (different test results) and continuous (age, yet some more test results which - in some cases complicated by a continuous result or N/A). If the outcome were continuous I think I'd be looking at using a mixed effects linear model so I could determine if there were any effects or interactions with the various secondary independent variables. Is there something equivalent I can do for ordinal and/or binary outcome? Thanks for your advice! Ric * * * **Edit1** (in response to answer below by Macro): I did read your linked answer for more discussion re: random effects vs GEE. However, I'm honestly still unclear on when you would want an individual vs a population based coefficient. To address my data specifically, we want to determine how well, if at all, the hypothesized predictor can predict the outcome. If we find it is a good predictor then we would want to be able to say when testing new patients in the future that given they are positive for test X, they are Z times more likely to become positive for outcome Y -- this will affect the course of treatment. So would I want a random effects or a GEE model to answer this question? What would be a similar scenario where I would want the other model? As a preliminary experiment I did a Fisher's Exact test on X vs Y (ignoring the fact that many observations are repeated measures) and found evidence for a possible trend but not a statistically significant effect (p ~ 0.1). I repeated this after splitting the data into 2 sets based on research centre. Centre 1: p=1, odds ratio ~ 1. Centre 2: p=0.02, odds ratio ~ 3. Other independent factors may influence whether X predicts Y or may correlate with Centre. Is logistic regression still an appropriate method to explorer these effects on the predictive effect of X? If it just tells me that X is not predictive but some other factor is, while interesting, it would not answer our primary question. If it tells us that X is predictive, but only for Centre 2, and especially for older patients at early times since treatment, that is the sort of thing we are looking for. Thanks! * * * **Edit2** (interpreting results?): So I've been playing with this the last few days but still having trouble making sense of it. I've found some helpful texts on logistic regression using glm() but nothing on lme4 that I was able to fully grok. Perhaps if you know of a good tutorial on the subject... In any case, to fill in some details on the data set here's a table summarizing the variable structures:               variable  type      levels   missing values     id        factor      76           0     inst      factor       2           0     d1        Date        NA           0     d2        integer     NA           0     d3        Date        NA           0     d4        integer     NA           0     a         integer     NA           2     b         integer     NA           6     c         factor       2           2     d         factor       2           6     e         factor       2           7     f         factor       2           2     g         logical     NA           6     h         integer     NA          12     i         numeric     NA          43     j         factor       2           0     x         factor       2           0     y         factor       2           0      The primary hypothesis is that x can predict y. y is the response variable. I initially removed any observations in which x or y were missing. 192 observations remained for 76 patients. There are missing values, however, for several of the variables: especially h and i, which based on some preliminary tests look like they might be important. There are 50 records in which some variable has a missing value, leaving 142 records with no missing values in any variable. As far as I understand, glmer() has no way of dealing with missing values and by default will omit the entire record if it has any missing values in the variables being used. I think this may be problematic due to the relatively small sample and effect sizes so I'm not sure what the best way of dealing with this would be. "id" is the patient id is what I am trying to use as the random effect with all other variables being fixed effects. I've found putting too many variables in the formula will cause glmer to fail to converge. Below are the commands I used for a few simpler models and results: command:               lr1   <- glmer(formula = y ~ x + 1|id, data = s.dat, family = binomial);      result:               Generalized linear mixed model fit by the Laplace approximation      Formula: y ~ x + 1 | id         Data: s.dat         AIC BIC logLik deviance      238.9 252 -115.4    230.9     Random effects:      Groups Name        Variance Std.Dev. Corr        id     (Intercept) 0.97379  0.98681                     xpositive   0.01454  0.12058  1.000      Number of obs: 197, groups: id, 76          Fixed effects:                 Estimate Std. Error z value Pr(>|z|)         (Intercept)  -1.1343     0.2121  -5.349 8.85e-08 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1       command:               lr2   <- glmer(formula = y ~ x * inst + 1|id, data = s.dat, family = binomial);      result:               Generalized linear mixed model fit by the Laplace approximation      Formula: y ~ x * inst + 1 | id         Data: s.dat         AIC   BIC logLik deviance      252.7 288.8 -115.4    230.7     Random effects:      Groups Name             Variance Std.Dev. Corr                       id     (Intercept)      1.016604 1.00827                                    xpositive        0.072397 0.26907   1.000                            instI2           1.626269 1.27525  -0.683 -0.683                     xpositive:instI2 0.135144 0.36762  -0.817 -0.817  0.137      Number of obs: 197, groups: id, 76          Fixed effects:                 Estimate Std. Error z value Pr(>|z|)         (Intercept)  -1.1224     0.2099  -5.347 8.96e-08 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1       command:               lr3   <- glmer(formula = y ~ x * inst + h + i + 1|id, data = s.dat, family = binomial);      result:               Warning message:     In mer_finalize(ans) : iteration limit reached without convergence (9)      command:               lr4   <- glmer(formula = y ~ x + h + i + 1|id, data = s.dat, family = binomial);      result:               Generalized linear mixed model fit by the Laplace approximation      Formula: y ~ x + h + i + 1 | id         Data: s.dat         AIC   BIC logLik deviance      184.6 217.1 -81.31    162.6     Random effects:      Groups Name        Variance   Std.Dev.   Corr                       id     (Intercept) 1.1663e+01 3.41508554                                   xpositive   1.2627e+00 1.12369065 -0.284                            h           7.9415e-02 0.28180626 -0.998  0.229                     i           2.1064e-07 0.00045896  0.107  0.923 -0.164      Number of obs: 142, groups: id, 57          Fixed effects:                 Estimate Std. Error z value Pr(>|z|)         (Intercept)  -1.3359     0.2792  -4.785 1.71e-06 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1       I don't really understand what this is telling me, or why all my variables are showing up as random effects instead of fixed effects. Any enlightenment will be greatly appreciated! Thanks!


title: How do test whether two multivariate distributions are sampled from the same underlying population?

text: Say you are given two multivariate data sets, say an old one and a new one, and that they are supposed to have been generated by the same process(that you have no model for) but perhaps, somewhere along the line of collecting/creating the data, something went awry. You wouldn't want to use the new data as, say, a validation set for the old data or to add to the old data. You can do a bunch of 1-d stats (per variable), e.g Wilcoxon rank sum, and try some multiple test correction but I'm not sure that's optimal (to capture the intricacies of multivariate data let alone multi-test issues). One way is to use a classifier and see if you can discriminate between the two datasets (given an optimal classifier that's optimal). That does seem to work but still a) perhpas there's a better way b) It's not really designed to tell you why it's different (if nothing else it will use the best predictors and possibly miss other good predictors that were subsumbed by the better ones)


title: Complex Dependent Variable

text: I have been running a regression analysis in SPSS where my dependent variable is a composite. By this I mean that it is built up of 7 different items. How can I make it one dependent variable? Because in SPSS you have to give only one dependent variable but I have 7 items for one dependent variable. Is there a way to do this in SPSS?


title: Confidence interval for binomial distribution

text: I know how to plot confidence intervals (CI) to a plot based on their values, but I'm not sure if there's a proper way to plot CI on a plot based on fraction above/below; e.g.               #dummy data     d1=c(20,30,50,80,70,40,4,7,9,11,14)     d2=c(22,32,51,90,100,30,14,71,19,12,1)     data=list(d1,d2)     #calculating the fraction less than the threshold 40     fr<-lapply(1:length(data), function(i)       unlist(sum(data[[i]]<=40)/length(data[[i]])))     plot(c(1,2), unlist(fr), type='l')      How would one calculate and plot CI for such a binomial distribution? thanks.


title: Survey Analysis

text: I have survey data from approximately 500 companies, and on one question they were asked to rate their priorities from 1 to 8. It's a forced ranking so you can't answer with two 8s. The issue that I'm now facing is to make some logical groupings or of the companies analysis. What are typical strategies for this sort of analysis? I've tried factor analysis which didn't seem to yield good results based on goodness of fit measures. The lack of fit could certainly be a function of my unfamiliarity with factor analysis. I've also looked at clustering algorithms like kmeans and while it provided reasonable groupings information criteria like BIC suggest that 1 cluster is the best fit. Thoughts?


title: Predictive accuracy and correlation of inputs

text: I am reading Judgement under Uncerntainty and in p.65 it is stated that in the normal linear model, correlation of input variables decreases the predictive accuracy (in contrast to the human perception which is exactly the opposite). I can't put a specific quote because the meaning is scattered over the whole page. I am having trouble understanding this statement. Can anyone help?


title: How do you estimate the predicted probability of an integer value from a negative binomial regression equation?

text: I'm trying to estimate the predicted probabilities of an observation being a particular integer, $y$, after a negative binomial regression model. Long's _Regression models for categorical and limited dependent variables_ gives this predicted probability as (pg.237): $$ \hat{\text{Pr}}(y \mid x) = \frac{ \Gamma(y + \hat{a}^{-1}) }{ y!\Gamma(\hat{a}^{-1}) } \left( \frac{\hat{a}^{-1}}{\hat{a}^{-1}+\hat{\mu}}\right)^{\hat{a}^{-1}} \left( \frac{\hat{\mu}}{\hat{a}^{-1}+\hat{\mu}} \right)^y $$ Where $\hat{\mu}$ is the predicted mean of the variable, $\hat{a}$ is the dispersion estimate, and $\Gamma$ is the Gamma function. Now, my question is the statistical software I use takes both a _shape_ and a _scale_ parameter for the $\Gamma$ distribution, so I am confused as to how to actually estimate the predicted probabilities for any particular integer $y$. In the above equation, what does Long expect me to supply as the shape and the scale for the $\Gamma$ function?


title: What is the amount of certainty in this yes/no test?

text: A child must learn a poem by heart. The poem has 200 lines. To test the child, the teacher asks the child to complete ten lines of the poem given the first half of the sentence. If the child gets, say, seven lines correct, the teacher assumes that the child knows 140 lines of the poem correctly. What is the uncertainty in the teacher's guess of the child's knowledge? In general, given $n$ things to learn, which may be either right or wrong, what uncertainty is associated with a test which tests $m$ of them?


title: Why X-process is called a process?

text: I have recently learnt about kernels in machine learning. And I have been introduced to many different processes e.g. Gaussian process, Wiener process. Now my question is why a set of functions has been named as a process? For example this is the definition of Gaussian process: Let $\mu: X \to R$ be any function, $k: X \times X\to R$ be a Mercer kernel. A Gaussian process $p(f) = GP(f;μ,k)$ is a probability distribution over the function $f : X \to R$, such that every finite restriction to function values $f_X :=[f_{x_1},...,f_{x_N}]$ is a Gaussian distribution $p(f_X)=N(f_X;μ_X,k_{XX})$. Which parts exactly associates to a process??


title: Card deck reshuffling

text: Suppose one has a deck of cards. If one reshuffles it using a discrete probability model (based on combinatorics) one will still have the same probability of getting a particular card on each draw. This seems to be a little bit naïve for real cards, since reshuffling is somewhat non-random. Suppose you know the ordering of the cards in the deck before reshuffling. Is it possible to tackle non-randomness of reshuffling to get conditional (on original ordering) probabilities of getting a particular card from the deck?


title: Units for silhouette measure

text: While drawing graphs that compare the silhouette measure of different clustering algorithms, what unit should I specify for the silhouette width?


title: How can I estimate that a customer wont buy again?

text: I have a history of purchases from different clients and I'm trying to calculate a threshold in days after which the probability of the client making another purchase is below 50% (asumming that this means that it's more improbable for him to purchase again). My current approach was to create an histogram with the number of days since last purchase, but I don't even know what to do with that data. I'm honestly at a complete loss. I don't even know the correct terms for searching the web for this. Any pointers are appreciated.


title: How to handle the error of glmnet package for non-positive lambda?

text: I'm using `glmnet` package to learn regression models,it works fine, but for some models, I face an error and my script stops running. Here is my effort:               # goal : build a model for every row of T     # T: matrix     #MI: matrix feature matrix          model <- vector("list", nrow(T))     for(i in 1: nrow(T)){       x <- MI       model[[i]] <-cv.glmnet(x, t(T[i,]), upper=0.0, parallel = TRUE)       print(i)   # keep track of progress of script     }      Now I run it               source("myscript.R")     # [1]1     # [1] 2     # [1] 3      .      .      .      .          # [1] 498     # [1] 499     # [1] 500     # [1] 501     # [1] 502     # Error in if (any(lambda < 0)) stop("lambdas should be non-negative") :      #   missing value where TRUE/FALSE needed      Would someone help me, to add some conditions to the script, which if I faced with this error, just skip that model and go to build other models.


title: In CFA model, why observed measures operate as dependent variables?

text: I am curious that why in CFA model, observed measures operate as dependent variables. In my understanding, observed measures served as data point that estimating other parameters, such as latent variables and errors. Thanks,


title: statistics t test sample size

text: in the formula for t test, t = (mean1 - mean2)/sqrt[(var1 + var2)/N]   1. is N the total number of subjects in both groups being compared or the number of subjects in each group or condition?   2. df = N-2 is N the total number of subjects in both groups being compared or the number of subjects in each group or condition?


title: Gibbs Sampling Detecting Change point in time series

text: I was reading through this one page paper on using Gibbs sampling for detecting a change point in a time series like data. While I understand the part where the $\lambda$ and $\phi$ are chosen from a gamma distribution I do not understand how the author generated the posterior distribution for $k$. My thinking is, if the first $k$ points came from a Poisson distribution then the probability of $P(k|\lambda,\phi,Y_i)$ $$ P(k|\lambda,\phi,Y_i) = \Pi_{i=1}^{k}\frac{e^{-\lambda}\lambda^{Y_i}}{Y_i!} \times \Pi_{i=k+1}^{N}\frac{e^{-\phi}\phi^{Y_i}}{Y_i!} $$ The author shows an expression of the form $\frac{A}{B}$ where $B$ is the normalizing term. This part I understand. However it is unclear to me how he found the numerator term which appears to be a log of $$ P(k|\lambda,\phi,Y_i) = \frac{\Pi_{i=1}^{k}\frac{e^{-\lambda}\lambda^{Y_i}}{Y_i!}}{\Pi_{i=k+1}^{N}\frac{e^{-\phi}\phi^{Y_i}}{Y_i!}} $$ How does he arrive at that? Any help appreciated. Thanks in advance


title: What distribution results in adding two Pareto distributions

text: I'm wondering what distribution results in adding two (or more) type-one Pareto distributions of the form $x^{-\alpha}$. Experimentally, it looks like a two-mode power-law, asymptotic to the difference of alphas.


title: Calculation of variance in boundary cases of constrained likelihood

text: I am computing the MLE values of some parameters of a likelihood function, subject to non-negativity constraints. Given these constraints, some of the MLE values for my parameters are exactly in the boundary (i.e. they are equal to zero). Is it still valid to estimate the variance of the MLEs using the Hessian of the likelihood function, even in these boundary cases?


title: Endogenous treatment effects: exogenous covariates

text: I am interested in estimating an endogenous treatment effects model of the following form: \begin{eqnarray} Y_i = \alpha + \beta_x X_i + \beta_{z1} Z_{1i} + e_i \\\ X_i = a + \beta_{z2} Z_{2i} + v_i \end{eqnarray} where $Y$ is a continuous variable, $X$ is a binary variable that is endogenous in the first equation and $Z_1$ and $Z_2$ are exogenous variables. Furthermore, I have reason to believe that $Z_3$, a third exogenous variable, is a determinant of both $X$ and $Y$. My question is the following. Should $Z_3$ appear in the second equation that explains $X$ only, or should it appear in the first equation that explains $Y$ as well?


title: Strange intercept using lmp() in lmPerm library in R

text: Please consider this data set:               y <- c(2, 4, 6)     x <- c(1, 2, 3)      Now calculate a linear model using lmp():               library(lmPerm)      lmp(y ~ x)      What I get, but don't understand, is this:               Coefficients:     (Intercept)            x                 4            2        My question: Why on earth is the intercept given as 4? It should be 0. lm() seems to agree:               lm(y ~ x)          Coefficients:     (Intercept)            x                 0            2        This is just an example with a small data set. I get similarly disturbing results on my actual data set, even if I tell lmp() to perform 100000 permutations and to stop at Ca=0.0001.


title: Data analysis tool with semantic capabilities?

text: I am currently working with quantitative data in the form of user-provided keywords which means the same idea can be expressed differently (e.g. "teachers, teaching, university, students" are related to the same "Education" subject area). I am currently using excel which is totally inadequate as Excel is incapable of understanding such things. I tried working around this problem by tagging keywords based on reverse dictionaries lookups which works OK, however Excel just dies when these lookups are applied to large spreadsheets (290K rows), therefore I need another tool. **Do you guys know about a data analysis tool with semantics capabilities?** _I looked at tools such as Tableau, Pentaho, and Jaspersoft but couldn't see that semantic layer._


title: In a moving average model, how do I calculate the errors?

text: I am programming an ARIMA model. The MA part of the model uses past forecast errors. How do I calculate these errors? Should I use any other model?


title: What is a "discriminant function" and how to interpret it?

text: I want to test a model using discriminant function analysis. My question, as the title states, is very basic: _What is a discriminant function?_ That is, how can I interpret the different discriminant functions? I'll briefly explain this question. I have several continuous variables, and a categorical variable (with $N = 4$ groups) that I want to predict. Following a discriminant function analysis in SPSS, the output states $N-1$ discriminant functions, that is $3$ discriminant functions. My question is how I can interpret these functions, so that I can understand how my dependent variables predict the categorical variable.


title: Obtaining error bars for a difference

text: I have a pair of data sets, and I want to know whether the medians of these two sets differ significantly. I've calculated the median of each set and subtracted them from each other. How can I determine appropriate error bars for that point? I was thinking I would take the median absolute deviation (MAD) of both and add them together. But should that be the entire length of the error bar? Or half of it? Or should I divide it by $\sqrt n$? Any advice is appreciated, thanks.


title: Two simple questions regarding GLM

text: I'm currently doing a modelling project. However, I haven't taken a bunch of statistics classes, so I have to teach myself generalized linear models. I'm reading _Generalized Linear Models for Insurance Data_ (Heller and de Jong, 2008, CUP), and I have two questions: 1\. On page 64, it says: > Given a response $y$, the generalized linear model is > $f(y)=c(y,\phi)\exp{\frac{y\theta - a(\theta)}{\phi}}$. The equation for > $f(y)$ specifies that the distribution of the response is in the exponential > family. Is that the equation for the distribution of $E[y_i|x_i]$ or some other thing? If it's the distribution for $y$ corresponding to a fixed $x_i$, is it possible that even if the plot of $y$ against $x$ looks like a straight line, I should still use GLM instead of simple regression? **update** : I guess I should clarify myself a little bit. Currently I have a dataset and my dependent variable is $y$. I made a histogram of $y$ (with frequency on y-axix) and it looks like a gamma curve fits well. Does that essentially imply that I should choose $f(y)$ to be gamma? I kinda doubt it because I suppose $y_i|_{X=x_i}$ and $Y$ are essentially two different things. I hope I'm not confusing you guys. 2\. The book suggests that when assuming response $y$ follows a gamma distribution, it is a common practise to use a logarithmic link function. I don't quite understand the reason behind that. Any suggestion would be great. Thanks!


title: Survival - comparing Kaplan-Meier curve to handful of points

text: I am investigating the effect of genomic dysmethylation on cancer survival time, with data from multiple different cancers with very different survival curves. Normally, I would split the cases into high and low dysmethylation, and do a log-rank test on the Kaplan-Meier curves. This isn't possible because for each cancer I have only 5 or so cases where methylation was measured. However, for every cancer, I do have 100+ cases where survival time was measured, even if methylation wasn't. Can I do a statistical test where I measure the survival times of the cases where methylation was measured against this "background" survival curve? Is it possible for me to combine the different groups, so I get a single answer across all the different cancers? Thanks, Steph


title: What is the False Discovery Rate curve?

text: I have heard of the " **false discovery rate curve** " (e.g here) but never seen an example. If I recall correctly from a conversation with a colleague, the `y-axis` in the FDR curve measures the FDR itself, defined as $FDR = \frac{FP}{TP+FP}$ (i.e. $1 - \text{Precision}$), but what goes into the x-axis? What is a normal shape for this curve, and when is it a good idea to use it?


title: Why do my ANOVA tables keep returning $\chi^{2}$ values of 1?

text: I'm using ANOVA to test for differences between different values of the same factor for a mixed effects model which I produced. My model is: `m2 <- lmer (ovsize ~ d.sheetratio + (1|nid), REML=FALSE)`. Following this, I subsetted the data so that each separate model described just data values for just one value of ratio.               me <- lmer (eovsize ~ eratio + (1|nide), REML=FALSE)     md <- lmer (dovsize ~ dratio + (1|nidd), REML=FALSE)     mc <- lmer (covsize ~ cratio + (1|nidc), REML=FALSE)     mb <- lmer (bovsize ~ bratio + (1|nidb), REML=FALSE)     ma <- lmer (aovsize ~ aratio + (1|nida), REML=FALSE)      At this point I get an error message which states that 'fixed-effect model matrix is rank deficient so dropping 1 column / coefficient'. Then when I do the ANOVA to compare each (e.g. `me` vs. `ma`) I get an output which looks like this:               Data:      Models:     mb: bovsize ~ bratio + (1 | nidb)     ma: aovsize ~ aratio + (1 | nida)        Df     AIC     BIC logLik deviance Chisq Chi Df Pr(>Chisq)     mb  3 -323.05 -316.06 164.53  -329.05                             ma  3 -235.34 -229.37 120.67  -241.34     0      0          1      How can $\chi^{2}$ be 1 when each model is different? Is there an error in the method that I've used? Also, is the error message I received important?


title: Empirical logit transformation on percentage data

text: I have already used the logit transform on my outcome variables (which are displayed in percentages). However, this obviously gives me -INF values and since my data includes a lot of zeros in some instances, this makes it hard to analyse. I have now tried an empirical logit transform, adding the smallest non-zero promotion to the numerator and denominator of my variables to remove the -INF values (as suggested in http://www.esajournals.org/doi/abs/10.1890/10-0340.1). However, now my data are very non-normal again. I have tried experimenting with error terms to add to the logit transform but since have had no luck. Is there any way I can find a value to add to my transformation to ensure normality?


title: Linear regression and ordinal data

text: I am running a multiple regression looking at whether TV viewing predicts waist circumference (WC). When I ran through the tests with my tutor we placed WC as the dependant and TV as independent, then ran it again with some of the potential confounders. However, the TV variable is ordinal (1- never, 2- 0 to 59 mins, 3- 1 to 2 hours, 4- 2 to 3 hours, 5- 3 to 4 hours, 6- 4 to 5 hours, and 7- 5+ hours). Should I be recoding these as dummy variables? I also plan to run regressions with similarly collected snack food consumption ordinal data. I'm hesitant to do it differently to my tutor (who is away for 2 weeks), and think it may be too many dummy variables, but is it still ok to interpret the results as I would a continuous dataset (in terms of mins/hours)? Many thanks, Nick.


title: Initial Parameters for VARMA models?

text: I want to estimate parameters of VARMA model using maximum likelihood estimation using real data. The problem I face with is that I don't know how to set the initial values for the parameters. I tried to choose randomly but the objective function is always undefined at initial point. Please, is there any theory about this? I googled and found out that using this program http://www.mathworks.com/help/econ/vgxvarx.html **it is not necessary to have initial values for any parameters to be estimated.**. But I don't understnd the theory behind this. Thank you.


title: Is it possible to assign class probability to a random forest prediction?

text: When making predictions with a random forest model, is it possible to associate the probability of a test case belonging to a class? For example, for a given test case, can we say that the probability of that test case belonging to the `setosa` class is 90%, `versicolor` class is 7% and `virginica`class is 3%? Is there any other way to somehow quantize the association?               library(caret)     data(iris);      inTrain <- createDataPartition(y=iris$Species, p=0.7, list=FALSE)     training <- iris[inTrain,]; testing <- iris[-inTrain,]     modFit <- train(Species~ .,data=training,method="rf",prox=TRUE)     pred <- predict(modFit,testing);                > pred      [1] setosa     setosa     setosa     setosa     setosa     setosa     setosa     setosa     setosa     setosa     setosa     setosa     setosa         [14] setosa     setosa     versicolor versicolor versicolor versicolor versicolor versicolor versicolor virginica  versicolor versicolor versicolor     [27] versicolor versicolor versicolor versicolor virginica  virginica  virginica  virginica  virginica  virginica  versicolor versicolor virginica      [40] virginica  virginica  virginica  virginica  virginica  virginica


title: How to make IPF code faster and more concise in R

text: Iterative proportional fitting is a way of adjusting 'internal cells' in a multidimensional matrix to optimise fit. It is also known as 'raking' and can be seen as a subset of 'entropy maximisation'. The purpose for which I use IPF is to allocated individuals to zones. My code iteratively adjusts a weight for each individual in each zone, depending on how representative the individual is of each zone. A practical introduction to the technique, with many links to other resources, can be found here. Now, my problem is that the code seems longer and slower than it needs to be. I have thought this for a while, but am only asking the question now because I've recently discovered 4 alternative R implementations of IPF:   1. This very recently published **ipfp** package implements IPF in C, so will be fast.   2. There is a longer-standing version of IPF within the long-standing categorical data analysis package, **cat**.   3. There is also a custom `ipf()` function that seems to be used by some, described on CRAN email list servers, dating back to 2008.   4. An implementation for the Alaska Department of Labor and Workforce Development  Options 2-4 are mentioned in a previous Cross Validated post that is very useful and which triggered this question. The top answer to this question proposed yet another method: `glm()`. Now, I can see that there are large potential benefits of using one or more of these options in my code to speed things up and to make things run faster. However, I can't see how to implement these new functions in the framework of my code. The problem, I suspect, is that I'm re-weighting individual-level data, whereas the other implementations are re-weighting aggregate data. So, to put the question more concisely: which option is best suited to my application? Any guidance (preferably based on the reproducible example) on implementation greatly appreciated. **Reproducible example** This reproducible example can be found here. Cloning the repo https://github.com/Robinlovelace/IPF-performance-testing enables running additional iterations.               # Minimal IPF example     start_time <- Sys.time()     num.its <- 10 # iterations          # Read-in data (manually to start, will use scripts in future)     c.names <- c("id", "age", "sex")     ind <- c(       1, "50+", "m",                     2, "50+", "m",                      3, "16-49", "m",                      4, "50+", "f",                      5, "16-49", "f")     ind <- matrix(ind, nrow = 5, byrow = T) # Convert long data into matrix, by row     ind <- data.frame(ind) # Convert this into a dataframe     names(ind) <- c.names # Add correct column names          # Read in the data in constraints     category.labels <- c("16-49", "50+",  "m", "f") # Age, and sex constraints     all.msim <- c(        8,        4,     6,    6       ) # each row represents an area, each column a category     all.msim <- matrix(all.msim, ncol = 4, byrow = T) # Convert long data into matrix, by row     all.msim <- data.frame(all.msim) # Convert this into a dataframe     names(all.msim) <- category.labels # Add labels     con1 <- all.msim[,1:2] ; con2 <- all.msim[,3:4]          ind.cat <- data.frame(cbind(model.matrix(~ind$sex - 1)[,c(2,1)], model.matrix(~ind$age - 1)))     names(ind.cat) <- category.labels          # create weights in 3D matrix (individuals, areas, iteration)     weights <- array(dim=c(nrow(ind),nrow(all.msim),num.cons+1))      weights[,,num.cons+1][] <- 1 # sets initial weights to 1     ini.ws <- weights[,,num.cons+1]          # convert survey data into aggregates to compare with census (3D matix)     ind.agg <- array(dim=c(nrow(all.msim),ncol(all.msim),num.cons+1))     for (i in 1:nrow(all.msim)){       ind.agg[i,,1]   <- colSums(ind.cat) * weights[1,i,num.cons+1]}          # re-weighting for constraint 1 via IPF      for (j in 1:nrow(all.msim)){       for(i in 1:ncol(con1)){      weights[which(ind.cat[,i] == 1),j,1] <- con1[j,i] /ind.agg[j,i,1]}}     for (i in 1:nrow(all.msim)){ # convert con1 weights back into aggregates       ind.agg[i,,2]   <- colSums(ind.cat * weights[,i,num.cons+1] * weights[,i,1])}          # test results for first row     ind.agg[1,1:2,2] - all.msim[1,1:2] # should be zero          # second constraint     for (j in 1:nrow(all.msim)){       for(i in 1:ncol(con2) + ncol(con1)){       weights[which(ind.cat[,i] == 1),j,2] <- all.msim[j,i] /ind.agg[j,i,2]}}       for (i in 1:nrow(all.msim)){ # convert con2 back into aggregate     ind.agg[i,,3]   <- colSums(ind.cat * weights[,i,num.cons+1] * weights[,i,1] * weights[,i,2])}          # for multiple iterations     wf <- array(dim=c(dim(weights), num.its, 1)) # array to store weights its, wei     indf <- array(dim=c(dim(ind.agg), num.its, 1))     wf[,,,1,1] <- weights      indf[,,,1,1] <- ind.agg          # loop for multiple iterations (run e2.R repeatedly, saving each time)     for(it in 2:num.its){     source(file="e2.R")     wf[,,,it,1] <- weights     indf[,,,it,1] <- ind.agg     }          # Analysis - in general, see analyis files     a.v <- as.vector(as.matrix(all.msim)) # constraints in long form, for cor     g.v <- as.vector(as.matrix(indf[,,1,1,1]))     cor(a.v,g.v)     (total_time <- Sys.time() - start_time)


title: "Exponential" distribution with negative scale parameter?

text: The standard exponential density is only defined for $x>0$ and scale parameter $\lambda > 0$ where the density is given by $f(x) = 1/ \lambda * exp(-x/ \lambda)$. My question is, does a similar density where $\lambda < 0$ and $x < 0$ make any sense and have any interpretation? Similarly, what about the case of $\lambda = 0$? Thanks!


title: How to choose which demographics to include in randomized trial regression?

text: I'm rusty on my stats. I am wondering what the prerequisites are for including covariates in my regression. I'm not looking for a prediction model but just the effect of the "treat". The sample was randomized and the demographics appear pretty similar with a couple of exceptions. Do I include all or only those that are the same in both treatment groups? Do I need to first do a statistical test to compare the two groups? Few of the demographics appear to make a significant difference in the regression, but if I exclude them this hurts the significance level of the treat coefficient.


title: Type of model to use

text: Suppose we are interested in an outcome variable $Y$. It can take four values $1,2,3,4$. These are categories. What would a probit link function do as opposed to a proportional odds model? Would a probit link model the outcome as $1 \ \text{versus} \ 2$, $2 \ \text{versus} \ 3$ and $3 \ \text{versus} \ 4$? Whereas the proportional odds model is cumulative?


title: How can I separate each of 100 observations into groups as determined by the data?

text: I have 3 covariates for 100 observations. How can I separate each of my 100 observations into groups as determined by the data. I was thinking clustering. However, apparently, I need more than 3 dimensions to do hierarchical clustering. Would some other clustering method work? How about PCA? I've attached the data as R output below. `dput` output               structure(c(3.87, 0.672, 0.7392, 6.471, 0.12294, 1.0857, 16.701,      0.2754, 0.17328, 8.076, 0.12222, 1.1796, 8.625, 1.3998, 0.07233,      3.933, 0.017484, 0.9189, 4.134, 0.7338, 2.9517, 5.091, 0.017136,      0.6318, 6.672, 3.012, 0.08214, 15.834, 0.7968, 0.27768, 3.954,      0.02046, 0.705, 9.465, 0.15444, 1.2702, 15.012, 0.4263, 2.262,      21.438, 0.9291, 0.3399, 20.076, 1.023, 8.289, 5.601, 0.007992,      0.984, 12.396, 0.4869, 2.343, 11.697, 0.4296, 0.4932, 8.247,      0.27063, 0.3408, 3.273, 0.03954, 0.16446, 4.59, 0.0011037, 0.8937,      11.196, 0.17538, 0.9594, 14.688, 0.13527, 0.3672, 2.8554, 0.0027594,      0.5943, 0.26472, 0.0004233, 0.3315, 1.5633, 0.0363, 0.5232, 5.766,      0.005901, 0.342, 10.578, 2.4477, 0.28872, 10.803, 1.0185, 1.3935,      5.352, 1.1967, 0.5316, 4.8, 0.00672, 2.6418, 3.081, 0.15525,      0.6873, 13.899, 0.19149, 0.4674, 11.439, 1.6521, 0.3867, 4.005,      0.008328, 0.3675, 5.7, 0.027999, 0.3486, 13.035, 0.21639, 0.7293,      8.706, 1.0833, 0.198, 5.871, 0.5655, 2.0367, 4.218, 0.011547,      0.24234, 2.1603, 0.0011748, 0.4569, 5.385, 0.005091, 0.666, 9.651,      0.4392, 1.101, 11.178, 0.4179, 0.7005, 8.151, 0.00516, 0.27696,      6.864, 0.018753, 0.303, 3.792, 0.01449, 1.8345, 6.834, 0.03339,      0.28896, 5.073, 0.012951, 0.5013, 3.132, 0.008892, 0.3207, 1.1841,      5.292e-05, 0.006795, 9.432, 0.324, 0.5916, 8.55, 2.4642, 0.9576,      3.588, 0.006912, 1.089, 6.396, 0.04818, 1.4448, 20.604, 0.363,      0.7401, 11.712, 0.03897, 1.9491, 11.682, 1.149, 2.217, 3.882,      0.15963, 5.916, 6.702, 0.3174, 1.6392, 7.188, 0.03582, 0.2646,      8.853, 0.7761, 2.2446, 18.915, 0.3993, 0.002736, 9.699, 0.16638,      0.6855, 9.423, 0.011793, 0.7986, 14.667, 2.5146, 0.28512, 5.919,      0.06705, 0.25305, 8.184, 0.005262, 0.6492, 1.4604, 0.0138, 0.1872,      8.604, 0.3057, 0.8052, 8.142, 0.017808, 0.9564, 2.2824, 0.000243,      0.2565, 6.012, 0.16425, 0.3969, 12.633, 0.9408, 1.7154, 9.396,      0.21945, 2.73, 13.479, 0.16236, 14.433, 9.612, 0.24222, 1.3275,      12.486, 2.1543, 0.08652, 0.0003612, 5.364e-07, 0.3144, 9.942,      1.3674, 1.326, 2.4621, 0.00019425, 0.6684, 1.6341, 0.0006165,      0.5124, 11.796, 0.9798, 3.243, 11.73, 0.4716, 1.0248, 5.133,      0.04527, 0.3078, 11.886, 2.6718, 1.158, 5.421, 0.06027, 1.7655,      6.69, 0.00783, 5.907, 11.832, 0.9534, 3.228, 1.0323, 0.0016356,      0.861, 6.774, 1.1001, 1.1811, 8.856, 0.4185, 1.3521, 11.877,      0.2754, 2.5563, 0.0024852, 1.4796e-05, 0.6741, 6.774, 0.336,      2.5017, 1.6425, 2.478e-05, 0.09243, 8.973, 0.25473, 0.9942, 13.245,      0.3234, 0.6711, 10.35, 0.5148, 1.0578, 14.556, 0.774, 0.9225), .Dim = c(3L,      100L), .Dimnames = list(c("A", "B", "C"), c("000162434", "000151547",      "000133688", "000123954", "000184599", "000122987", "000117559",      "000121528", "000192459", "000196759", "000172539", "000155583",      "000185889", "000143968", "000128617", "000185423", "000158324",      "000114797", "000126134", "000185624", "000123385", "000188299",      "000195142", "000194666", "000113189", "000182457", "000173324",      "000162459", "000141996", "000155516", "000148231", "000176159",      "000135131", "000186287", "000187355", "000199513", "000125251",      "000116237", "000188675", "000147224", "000198156", "000119366",      "000132841", "000123791", "000138154", "000149758", "000157127",      "000167763", "000113718", "000128418", "000148221", "000139836",      "000194814", "000199972", "000168968", "000198853", "000128498",      "000149484", "000196219", "000184178", "000144155", "000114251",      "000114264", "000131697", "000154146", "000163257", "000112289",      "000114416", "000195761", "000128348", "000144337", "000167126",      "000159175", "000172296", "000182932", "000198134", "000127718",      "000166651", "000196877", "000174415", "000131167", "000165476",      "000195958", "000189229", "000119255", "000165984", "000119118",      "000164273", "000199986", "000136544", "000124271", "000191248",      "000126459", "000143728", "000182847", "000162785", "000193387",      "000119516", "000199516", "000145424")))


title: How to get the true mean forecast using the Arima package with a Box-Cox transformation

text: In the Arima package, using a Box-Cox transformation give wrong results when later applied to the forecast method. For example, consider this data:               library(forecast)     data<-c(2,3,2,3,2,3)      And for the sake of simplicity, consider an ARIMA(0,0,0) model. (The mean of this series is 2.5.) The mean forecast made without a Box-Cox transformation is correct:               forecast(Arima(data,order=c(0,0,0)))$mean      [1] 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5      However if we use a Box-Cox transformation, such as a log transformation with lambda=0, the "mean" forecast is wrong:               forecast(Arima(data,order=c(0,0,0), lambda=0))$mean     [1] 2.44949 2.44949 2.44949 2.44949 2.44949 2.44949 2.44949 2.44949 2.44949 2.44949      It seems that to produce the mean forecast of Y=exp(X), it does E(Y)=exp(E(X)). Is there a way to correct this? Is there a package with a correct implementation of forecasts with Box-Cox transformations?


title: interpreting the coefficient on a logged independent variable from a probit model

text: I am using a double model with log transformed independent variables and have calculated average partial effects. Now I am now not sure how to interpret the coefficients; particularly those from the first stage which is supposed to be capturing the probability of participation using a probit model.


title: Repeated measures mixed modeling. Why is Time (day, second etc.) a fixed effect?

text: Let's say I have an experiment where I repeatedly measure something over Time (1:10), say 10 days. Here I simulate some data...               set.seed(101)          N = 10 # number of repeats     n = 10 # number of individuals          # Data     Time = rep(1:N,n) # measured over time      Note that I set the fixed effect of time to 0.               # Fixed effects     X = matrix(c(Time), ncol = 1)     B = c(0) # Set fixed effect of Time's coefficient to 0           # Random effects     id = factor(c(0L, cumsum(diff(Time) < 0))) # make id for each block of Time     Z = model.matrix(~id+id:Time) # create design matrix          u = c(rnorm(n),rnorm(n,2,1)) # random slopes and intercepts     e = rnorm(N*n) # error               y = X%*%B + Z%*%u + e # Model        Now let's say I fit two models:               library(lme4)     fit1 = lmer(y ~ 1 + Time + (1 + Time|id))     fit2 = lmer(y ~ 1 + (1+ Time|id))      Comparing the two with anova()               anova(fit1,fit2)     Models:     fit2: y ~ 1 + (1 + Time | id)     fit1: y ~ 1 + Time + (1 + Time | id)      Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)         fit2  5 380.80 393.83 -185.40   370.80                                  fit1  6 329.41 345.05 -158.71   317.41 53.389      1  2.736e-13 ***     ---     Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1      Here we see that `fit1`, which includes `Time` as a fixed effect, performs much better, but why is this true when its coefficient is set to 0? I have a feeling it has something to do with this line.               Z = model.matrix(~id+id:Time) # create design matrix      but my intuition is failing me.


title: Nonlinear functions of other features as new features in SVM model with RBF kernel

text: Can some one give me some conceptual insight on the potential advantages of disadvantages of adding features that are (nonlinear) functions of existing features in training an SVM model with an RBF kernel? For example, if I include area and perimeter as features then technically all of the information in something like the hydraulic radius (area/perimeter) is already in the feature set. However I think the relative value of adding features that may be nonlinear functions of other features may lie in the ability of the classification boundary to be nonlinear in a way that can capture higher level functional dependencies that lie in the new feature. For example, adding the new feature would be more advantageous in linear SVM than with an RBF kernel. Am I thinking about this correctly? And can someone point out the relative advantages or disadvantages of doing this with an RBF kernel? Would this be heading into the area of feature duplication? Thanks!


title: How to forecast with quantile regressoin

text: If you have three quantile regression models with taus of 0.25, 0.5 and 0.75 and their coefficients how do you use these models to forecast a set of data not used to calculate the coefficients. In other words how do you do a cross validation quantile regression.


title: Basic Question about the Latent Dirichlet Allocation Generative Model

text: So Here is the LDA Generative Model ![enter image description here](http://i.stack.imgur.com/oW48y.png) The $\alpha$ and $\beta$ nodes represent the parameters for two Dirichlet distributions. The $\theta$ and $\phi$ nodes represent the parameters for two multinomial distributions. My question is about the $Z$ node (the topic assignment for the word). It represents a sampled topic from the topic distribution for a document from $\theta$. So I don't think it is actually a random variable, it is a sample. But then why is it a node in the generative model if it is not a random variable? This has me confused. Is $Z$ a random variable or not?


title: Include insignificant predictor with low variance?

text: In my regression, I have a predictor (measured in proportions) with very low variance (.002). It had no significant effect on the outcome due to a very high standard error of the regression coefficient. When reporting on the effect in a paper, is it right to say there was no effect of my predictor in my study? Even though I know the lack of effect is most likely due to the low variance, and I add that thought in the paper? It feels like I might mislead the non-attentive readers, or I am reporting useless stuff. Maybe it is best practice to discard this predictor from the model, and note in the paper that it was impossible to test this effect due to the data limitations? And as extra I wanted to ask if it is possible to get a lower variance when taking the log of the values? I observed this with my variable.


title: How would you validate a random walk model?

text: I have used a random walk model and Gibbs sampling (more specifically RJAGS) in order to obtain posterior of the state given the observations. In this case the state is the true proportion of the population that will vote one way or the other. I have used diagnostics to make sure that the Markov chains have reached the targeted distribution but I am wondering what else do I need to do in order to validate the model. Thanks


title: Significant digit 354 rule from Particle Data Group

text: The particle data group (PDG) has a particular rule to round the value of a measurement. You can find the complete description at this link (page 13). Summarizing: > ... if the three highest order digits of the error lie between 100 and 354, > we round to two signiﬁcant digits. If they lie between 355 and 949, we round > to one signiﬁcant digit. Finally, if they lie between 950 and 999, we round > up to 1000 and keep two signiﬁcant digits. In all cases, the central value > is given with a precision that matches that of the error do you have an idea where this `354` rule come from?


title: Simple scheme to sample from the Bernoulli distribution

text: I was implementing a simple scheme of Bernoulli distribution sampler. $ X \sim B(p) $. I have a function that generates a uniform random number $r \in (0,1)$. Then, I set $ X = 1 $ if $p > r $, and $X =0$ otherwise. Is this correct?


title: How to cluster with binary, nominal, ordinal and continuous variables?

text: I am using SAS to do clustering analysis on a huge dataset. Since my dataset contains various types of variables, I am confused about the appropriate method to perform the analysis. Here are my questions:   1. How to deal with nominal variables?   2. I standardize the data to be analyzed using `PROC STANDARD (mean=0 std=1)` before doing `PROC CLUSTER`. However, the result will be sill based on binary variables, meaning the clusters I got are heavily divided based on binary variables used, such as gender. How can I deal with this issue?   3. If I make some variables have a larger std, will these variables have larger effect on the cluster generated, i.e., make my clusters more dependent on these variables?


title: Mixed effect logistic regression in R: choosing random effects

text: I conducted an experiment which measured a binary response for each subject. The subjects were in 1 of 3 groups. There were two other fixed factors, each of which were continuums (cont1, cont2) ranging from 0 to 10. In other words, for each step in cont1, there was a corresponding 0-10 step continuum (cont2). Cont1 refers to formant frequency of a vowel and cont2 refers to vowels duration. Therefore, there was a total of 121 stimuli repeated 7 times. Each subject provided 847 responses. I would like to know how the three groups differ in there responses for the two continua. Groups 1 and 2 are control groups (I know more or less how they behave), but I am interested in group 3 and if they perform more like group 1 or group 2. I have used the following model in R using lme4:               full.mod <- glmer(response ~ cont1+cont2+group+(random effects), data=df, family=binomial)      I think that I need the interaction between the three fixed effects to answer my research question; however I am not sure.               full.mod.int <- glmer(response ~ cont1*cont2*group+(random effects), data=df, family=binomial)      Therefore, my first question is whether full.mod or full.mod.int is the best choice. My experience with ANOVA leads me to believe that I need the interactions, however, being that there are 11 levels for each continuum, the output of coefficients is going to be enormous, and this makes me think that something is not right. With regard to random effects, I know that I should include random intercepts for each subject (1|subject). However, it is not clear to me whether or not I need a random effect for items too. In this case, both cont1 and cont2 are items. I have a column in my data farm called stimuli which gives information like: cont1_0_cont2_0, cont1_0_cont2_1, etc. Assuming this is what I should do, my model now looks like this:               full.mod.int <- glmer(response ~ cont1*cont2*group+(1|subjects)+(1|stimuli), data=df, family=binomial)      Would there be any benefit to adding random slopes? If so what would they be (cont1,cont2?)?


title: Hypothesis testing for two data series

text: I am currently pursuing research in management, but I have a serious problem with selecting the right statistical method. I have quarterly data for a couple of financial ratios (for example Return on Equity) for two groups of companies: those with CEOs with MBAs and those without MBAs. I want to test whether the two groups are statistically different. The first test which came to my mind is the student's t-test. My only concern is that this is a data series, so only the results from each quarter would be compared, and hence a t-test is not suitable. What would be the most appropriate test for this situation?


title: Multinomial logistic regression with "constant variable"

text: I am encountering some problems with setting up a multinomial logistics regression (stepwise). My data/model is the following:   * Dependent variable: strategy1/strategy2/no preference   * "Indenpendent" variable: risk perception value (prob*impact in percentage)   * Moderator: risk profile (likert scale)   * +some other control variables (demogrpahics etc.) Since the risk perception value is given, there is no variance in the values. I basically have two scenarios, each with different risk peception values (high_prob*high_impact and low_prob*high_impact) and these values are the same for all respondents. This is because the choice/preference for a certain strategy is supposed to be dependent on the (given) risk perception). However, when I run the regression for ONE of the scenarios with only this variable it gives me the following output: No significance level for the varaible but only for the intercept. For the significance level of the variable it only shows a "." (dot). So, SPSS does not give me a significance level for the variable (probably because of no variance), only for the intercept. Do I just have to accept it or is there any way to tackle this problem? Is this output even useful? Would it be better to generate a random perceived value in the range of 0.9 to 0.99 for e.g. high probability to get some variance? And then cluster them in two categories, thereby obtaining two perceived values for each scenario? Thank you very much in advance.


title: How to remove individual means and deviations for repeated measures

text: I currently have a experiment with a within subject design, where I repeatedly measure reaction times (RT). Reaction times for each category of reaction were measures multiple times. However now the problem is, that the mean of all reaction times for each subject may differ by quite a lot, as some person may react much slower altogether. The same goes for the standard deviation of each person. Now I was thinking I could clean up the data by removing individual means and deviations. Removing the mean of each persons RT is quite simple by just substracting the mean for that person from each RT for that person. However this leaves the difference in personal deviation in the sample. I could reduce this by dividing by the personal deviation, thus getting z-scores based on individual man and deviation for all reaction times for each person. By calculating the means of each z-score for each category, I then could see if categories lie above the mean, or below it and how much they differ in terms of the individual standard-deviation from the mean. Is this procedure acceptable to remove these individual factors from the analysis? Or is there another procedure, which should be used to clean up this kind of data?


title: Confusion related to a derivation

text: I was reading this paper http://cs.ru.nl/~perry/publications/2011/ICANN2011/groot-icann2011.pdf and I am a bit confused how this was derived $p(f|Y) \propto p(f)*p(Y|f) \propto exp^{(-\frac{1}{2}f^TK^{-1}f-\frac{1}{2}(Y-f)^T\Sigma^{-1}(Y-f))}$ is Gaussian processN(m,V) $ \propto exp^{-\frac{1}{2}f^TV^{-1}f + fV^{-1}m}$ with mean and covariance given by $m=(K^{-1}+\Sigma^{-1})^{-1}\Sigma^{-1}Y$ $V=(K^{-1}+\Sigma^{-1})^{-1}$ Is this the right way to derive it $p(f|Y) \propto p(f)*p(Y|f) \propto exp^{(-\frac{1}{2}f^TK^{-1}f-\frac{1}{2}(Y-f)^T\Sigma^{-1}(Y-f))}$ $\propto exp^{(-\frac{1}{2}f^TK^{-1}f-\frac{1}{2}(f-Y)^T\Sigma^{-1}(f-Y))}$ product of two gaussian so the mean and covariance of the product of two gaussians ![enter image description here](http://i.stack.imgur.com/PVx2K.png) and I eventually get $m=(K^{-1}+\Sigma^{-1})^{-1}\Sigma^{-1}Y$ $V=(K^{-1}+\Sigma^{-1})^{-1}$ However, I still don't get how it was manipulated(the expression) to get $ \propto exp^{-\frac{1}{2}f^TV^{-1}f + fV^{-1}m}$


title: Confusion related to derivation of product of gaussian

text: I am confused. Given ![enter image description here](http://i.stack.imgur.com/WAUg6.png) how come this is derived?enter image description here ![enter image description here](http://i.stack.imgur.com/JQsVu.png)


title: qqPlot with Two Outliers

text: If I have a QQ plot with two extreme outliers (picture below) how should I interpret it? Do I drop the outliers? Can I treat it as normal?![enter image description here](http://i.stack.imgur.com/N3ONG.png)